title: Estimating the break point in a broken stick / piecewise linear model with random effects in R [code and output included]

text: Can someone please tell me how to have R estimate the break point in a piecewise linear model (as a fixed or random parameter), when I also need to estimate other random effects? I've included a toy example below that fits a hockey stick / broken stick regression with random slope variances and a random y-intercept variance for a break point of 4. I want to estimate the break point instead of specifying it. It could be a random effect (preferable) or a fixed effect.               library(lme4)     str(sleepstudy)          #Basis functions     bp = 4     b1 <- function(x, bp) ifelse(x < bp, bp - x, 0)     b2 <- function(x, bp) ifelse(x < bp, 0, x - bp)          #Mixed effects model with break point = 4     (mod <- lmer(Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject), data = sleepstudy))          #Plot with break point = 4     xyplot(             Reaction ~ Days | Subject, sleepstudy, aspect = "xy",             layout = c(6,3), type = c("g", "p", "r"),             xlab = "Days of sleep deprivation",             ylab = "Average reaction time (ms)",             panel = function(x,y) {             panel.points(x,y)             panel.lmline(x,y)             pred <- predict(lm(y ~ b1(x, bp) + b2(x, bp)), newdata = data.frame(x = 0:9))                 panel.lines(0:9, pred, lwd=1, lty=2, col="red")             }         )      Output:               Linear mixed model fit by REML      Formula: Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject)         Data: sleepstudy        AIC  BIC logLik deviance REMLdev      1751 1783 -865.6     1744    1731     Random effects:      Groups   Name         Variance Std.Dev. Corr                Subject  (Intercept)  1709.489 41.3460                               b1(Days, bp)   90.238  9.4994  -0.797                       b2(Days, bp)   59.348  7.7038   0.118 -0.008       Residual               563.030 23.7283                     Number of obs: 180, groups: Subject, 18          Fixed effects:                  Estimate Std. Error t value     (Intercept)   289.725     10.350  27.994     b1(Days, bp)   -8.781      2.721  -3.227     b2(Days, bp)   11.710      2.184   5.362          Correlation of Fixed Effects:                 (Intr) b1(D,b     b1(Days,bp) -0.761            b2(Days,bp) -0.054  0.181      ![Broken stick regression fit to each individual](http://i.stack.imgur.com/HnAfg.jpg)


title: Why does ${\rm cov}(Y_{ij}, Y_{ik}) = \sigma_u^2$ in the model $Y_{ij} = \mu + u_i + e_{ij}$?

text: Consider $Y_{ij} = \mu + u_i + e_{ij}$ where $\mu$ is a constant and $u_i ∼ N(0,\sigma_u^2 )$ is the subject $i$ random effect which is independent of the error term $e_{ij} ∼ N(0,σ_e^2)$. I can't figure out why $${\rm cov}(Y_{ij}, Y_{ik}) = \sigma_u^2$$ Can anyone please explain?


title: Random matrices with constraints on row and column length

text: I need to generate random non-square matrices with $R$ rows and $C$ columns, elements randomly distributed with mean = 0, and constrained such that the length (L2 norm) of each row is $1$ and the length of each column is $\sqrt{\frac{R}{C}}$. Equivalently, the sum of square values is 1 for each row and $\frac{R}{C}$ for each column. So far I have found one way to achieve this: simply initialize the matrix elements randomly (e.g. from a uniform, normal, or laplace distribution with zero mean and arbitrary variance), then alternately normalize rows and columns to ${\rm length} = 1$, ending with row normalization. This seems to converge to the desired result fairly quickly (e.g. for $R=40$ and $C=80$, variance of column length is typically ~ $~0.00001$ after $2$ iterations), but I'm not sure if I can depend on this fast convergence rate in general (for various matrix dimensions and initial element distributions). My question is this: is there a way to achieve the desired result (${\rm row \ lengths} = 1$, ${\rm column \ lengths} = \sqrt{\frac{R}{C}}$) directly without iterating between row/column normalization? E.g. something like the algorithm for normalizing a random vector (initialize elements randomly, measure sum of square values, then scale each element by a common scalar). If not, is there a simple characterization for the convergence rate (e.g. num iterations until error $< \epsilon$) of the iterative method described above?


title: Do expits reveal effect sizes with logistic regression?

text: I used a logistic regression model (glme4) to determine whether two different types of words (emotional vs. unemotional) are recognized more accurately (accurate=0, error=1) when shown in a particular color (green, red, blue). The logistic model shows a significant difference between red and blue with a coefficient of .064. This difference was not influenced by word type (a coefficient = .0042 for the corresponding statistical interaction). My question is whether the “expit” function               expit <- function(x) { 1/(1+exp(-x)) }      is sufficient for the specification of the effect size (change in probability) for each of the two coefficients. How can the critical coefficients of .064 and .0042 be expressed in percentage values (effect sizes)? Numerically, the mean error rates for red and blue are 12% and 15%, respectively. I used the lme4 function               lmer(Correct ~ WordType * Color + (1| Participants), data = df)      with `contrasts(df$WordType) <- contr.sdif(2))` and `contrasts(df$Color) <- contr.sdif(3))`. Results are given below:               Fixed effect    Estimate    SE       Z          intercept       -1.86       .144    -12.96     WdTpe2-1          .058      .055      1.05     Color21           .006      .028       .228     Color 32          .064      .017      3.79     WdType:Color21   -.008      .055      -.14      WdType:Color32    .0042     .033       .127


title: Confidence intervals for difference in time series

text: I have a stochastic model used to simulate time series of some process. I am interested in the effect of changing one parameter to a specific value and want to show the difference between the time series (say model A and model B) and some sort of simulation based confidence interval. I have been simply running a bunch of simulations from model A and a bunch from model B and then subtracting the medians at each time point to find the median difference throughout time. I used the same approach to find the 2.5 and 97.5 quantiles. This seems like a very conservative approach since I am not considering each time series jointly (e.g., each point is considered independent of all others at previous and future times). Is there a better way to do this?


title: Ordering in VAR models

text: In Enders' 'Applied Econometric Time Series', I repeatedly stumbled upon the notion of the "ordering of a VAR model" and I am not sure I understand the concept right. As far as I understand it, the concept relates to the problem of recovering the true $\epsilon$ of the underlying structural equations. Looking e.g. at the following VAR in standard form: $y_t = a_{10} + a_{11}y_{t-1} + a_{12}z_{t-1} + e_{1t}$ $z_t = a_{20} + a_{21}y_{t-1} + a_{22}z_{t-1} + e_{2t}$ The problem is that the two errors are actually composites of the underlying shocks $\epsilon_{yt}$ and $\epsilon_{zt}$ of the structural equations, which we do not know. Here, Enders mentions that "(...) if the correlation coefficient between the $e_{1t}$ and $e_{2t}$ is low, the ordering is not likely to be important", and vice versa. My problem is that I do not grasp the concept of ordering and how it would affect the identification problem of the disturbances in the structural model. Does ordering refer to the order of the independent variables in the VAR?


How much do self organising maps suffer from local minima problems? I assume that the original orientation of the grid can have major orientation on the resulting shape of the grid, but does this have large impacts on the fit (e.g. mean square difference from data to assigned grid points), or are the differences fairly small?


title: When analyzing dyadic data, can unpaired cases be used in any part of the analysis?

text: I am currently in the process of analyzing a data set comprised of manager- subordinate dyads. Data were collected cross-sectionally and the data set contains some of the same variables collected from both members of the dyad (e.g., age, relationship tenure, interpersonal trust) as well as some unique variables for each member. Unfortunately, I have a large number of participants whose partners did not fill out their surveys. I was wondering if these unpaired cases could be of any utility in my analysis? For example, when conducting an EFA or CFA on a scale, could I include all managers (including the unpaired individuals), or is this bad practice since the unpaired individuals will not be a part of my regression and/or path analyses? I have been reading Dyadic Data analysis by David A. Kenny et al., but I have been unable to find an answer to this question so far.


title: Making horizontal max-min-average chart in Excel

text: I'm trying to make a min, max, average chart with horizontal lines, where the vertical axis represents years in the series, something similar to this chart at ![enter image description here](http://i.stack.imgur.com/69lBP.gif) I've already created the data for the chart, i just don't know how to proceed


title: Sampler method to choose in Monte Carlo Markov chain estimation

text: When estimating the posterior via MCMC, are there guidelines on the best sampling method to use depending on the nature of the model? There are a variety of forms of MCMC - the Gibbs sampler, the Metropolis algorithm, and the Metropolis-Hastings algorithm for example. Gilks et al. (1996) has a nice discussion of the MCMC variations. My particular use-case involves estimating a posterior by mixing observed sample security returns with a prior equilibirum model of security returns. Return series are often heteroskedastic and have cross-sectional correlations and are non-stationary across market regimes.


title: How can I calculate a probability from a likelihood, e.g. in the Metropolis-Hastings algorithm?

text: This is a follow-up to my previous question, how can I compute a posterior density estimate from a prior and a likelihood I am having difficulty understanding how it is possible to calculate the probability of a parameter set given its likelihood. For example, in the answer to the previous question, the likelihood is used to weight the probability of sampling from the prior distributions. However, the answer suggested I use a sequential sampling algorithm, so I have chosen to try a Metropolis-Hastings approach. However, it is not clear to me how I can use likelihoods to determine the probability that I will accept a proposed parameter set. I can empirically estimate (using the R `fitdistr` function) that the likelihood is approximately $L\sim\textrm{log-}N(8, 0.7)$, I suppose that I could calculate the probability as `plnorm(L, 8, 0.7)` and the probability of acceptance as $P(L_\textrm{new})/P(L_\textrm{old})$. But is there a way to do this without assuming that the likelihood has some particular distribution? Or is it standard practice to make such an assumption?


title: Reporting the relationship between variables: correlation or regression?

text: I'm doing my dissertation in Accounting and I developed a questionnaire to test certain hypotheses on the relationship between independent variables and a dependent variable. All were measured on a 5-point Likert scale. I am using SPSS to analyse the data. The variables are all non-normal (I used Kolmogorov- Smirnov test) and so I will be using Spearman correlation to test the relationships. I was wondering whether this is enough? Should I use some form of regression as well? And if I do, do I still include the Spearman results? Do studies usually include both? I have not studied any statistics so this is quite new to me.


title: How to estimate the given function using Rao Blackwellization approach?

text: I have a function $X$ which is lognormal (0,1) and then another function $\log Y = 4 + 2 \log(X) + \epsilon$ where $\epsilon \sim \mathcal N(0,1)$ I want to estimate $E(Y|X)$ as a Rao Blackwellized estimator. I tried to attempt the problem like this:               Nsim <- 10000     x <- rlnorm(Nsim)     z <- rnorm(Nsim,4+2*log(x),1)      Standard montecarlo estimator would be               sum(epx(z))/Nsim      I am confused about Rao Blackwellized estimator. Will it be just               p <- 4+2 * log(x)     sum(p)/Nsim


title: What are the regularity conditions for Likelihood Ratio test

text: Could anyone please tell me what the regularity conditions are for the asymptotic distribution of Likelihood Ratio test? Everywhere I look, it is written 'Under the regularity conditions' or 'under the probabilistic regularities'. What are the conditions exactly? That the first and second log-likelihood derivatives exist and Information matrix not be zero? Or something else entirely?


title: Original source for the "play tennis" dataset

text: A famous toy example in machine learning, especially with learning decision trees, is the well known "play tennis" dataset. Is there an official source for the dataset which could be quoted in scientific work?


title: How could I know that 3 sets of data have the same/different distribution?

text: I have 3 sets of data and I would like to know that what test I could use to find out that whether those sets of data have the same distribution. Edit 1: If I use hypothesis tests to show there are no difference between means and variances of those data sets, does it mean that they have the same distribution? Thank you.


title: Designing Computer Science experiments (Hypotheses and errors?)

text: For my master thesis I want to set-up several experiments. One thing my professor has been complaining about is that most computer science experiments are lacking scientifically. To not fall into this trap I want to properly design my experiments. Unfortunately I could not find much literature on CS experiment design. (Except for D. Feitelson's "Experimental Computer Science: The Need for a Cultural Change" available here which is more of an appeal than a guide.) I remember that in statistics people often use a H_0 and H_1 hypothesis together with the Type I and Type II error. I wonder if this is applicable in the field of computer science. A quick example that I could think of: > H_0: Both methods create the same output in this scenario > > H_1: One method creates a worse output in this scenario (worse: according to > metric X) But then how would I calculate the Type I and Type II errors? In a speed- benchmark I could run the experiments multiple times and then possibly calculate how likely it is that the benchmark result is still wrong. However often I just want to compare the output of two deterministic methods. Measuring the quality according to a certain metric. This result doesn't change when running the experiment multiple times. What do I do then? Is the null hypothesis approach not a good fit for these kind of experiments?


title: Is it appropriate to suppress the intercept in a regression graph?

text: In the event of a non-significant interaction between X and M on Y, is it sound to suppress the intercept in the graph that illustrates the relationship between X and Y (X being significant). The moderator line in this case will pass through the origin and represent each level of the factor (M) at zero (i.e. male and female for gender). SPSS provides this option under scatterplot.


title: Classifying new observations into two bivariate categories

text: I have two bivariate distributions $A$ and $B$. For some new observation $x$, I would like to calculate the probability that it belongs to $A$ and not $B$. I can assume that $A$ and $B$ are bivariate normal distributions (but if there is a method that doesn't assume that I'd like to know), but I can't assume their covariance matrices are the same. So far, I've calculated the mahalanobis distance from $x$ to $A$ and from $x$ to $B$, but I'm not sure where to go from there. I should also point out that I'll be doing this for many subjects for whom the distance between $A$ and $B$ will not be the same. If there is some way to take this into account when evaluating $P(x\in A ~\wedge \notin B)$, that'd be nice too.


title: Heuristic to determine the number of latent classes

text: Do anyone know if in the field of topic modeling and latent class analysis, it exists some heuristic to determine a good number of latent classes?


title: Discretize frequency of words that follow zipfian distribution

text: How could I discretize the frequency of words found in a corpus that follow a zipfian distribution? Are there standard methods? It should create bins of exponential-increasing size. My goal is to convert this continuous real frequency into multiple bins to create binary features in Conditional Random Fields (CRFs). I would like to be able to easily define how many bins/divisions I want, to be more or less precise as desired. For clarification, I have generated a list ordered by `freq` that gives me the frequency and cumulative frequency of each word in a studied corpus:               #1 word1 freq1 cumFreq1     #2 word2 freq2 cumFreq2     #3 ...      Then when generating features in a training set, for every wordx I want to assign it the feature-bin it belongs to or a `NOT_IN_DICTIONARY` label if the wordx was never found in the studied corpus. I want to create bins-labels so that they group words by similar frequency.


title: Permutational MANOVA and Mahalanobis distances in R

text: In this paper > Anderson, M.J. (2001). A new method for non-parametric multivariate analysis > of variance. _Austral Ecology_ , **26** , 32–46 I found very useful MANOVA approach that can, to some extent, cope with non- normality. I was wondering if it is possible and plausible to include Mahalanobis distances as similarity method. In fact my question is about R vegan package function `adonis` that implements this type of MANOVA but does not allow Mahalanobis distances as similarity measure. It allows though through `designdist` function the creation of any distance measure to be used in adonis, but I totally can't deal with that syntax. Any experiences regarding this?


title: Rolling Window Forecasts in R

text: I have a question: how do I use rolling window forecasts in R: I have 2 datasets:   1. monthly data which I downloaded from Google.    2. monthly data I downloaded from the CBS (central bureau of statistics in Holland) I want to test whether I can build a valid forecasting model, based on say 6years of Google Data, by using rolling window forecasts. My first intuition was to fit an ARIMA model on the Google Data and test this model on the actual CBS data by calculating the RMSE. But this is not efficient since I need to do this for every month, and I have a lot of months to analyse. Is there an automated way in R to do this?


title: Best approach for model selection Bayesian or cross-validation?

text: When trying to select among various models or the number of features to include for, say prediction I can think of two approaches.   1. Split the data into training and test sets. Better still, use bootstrapping or k-fold cross-validation. Train on the training set each time and calculate the error over the test set. Plot test error vs. number of parameters. Usually, you get something like this:![enter image description here](http://i.stack.imgur.com/6PrVd.jpg)   2. Compute the likelihood of the model by integrating over the values of the parameters. i.e., compute $\int_\theta P(D|\theta)P(\theta)d \theta$, and plotting this against the number of parameters. We then get something like this:![enter image description here](http://i.stack.imgur.com/TwHCk.jpg) So my questions are:   1. Are these approaches suitable for solving this problem (deciding how many parameters to include in your model, or selecting among a number of models)?   2. Are they equivalent? Probably not. Will they give the same optimal model under certain assumptions or in practice?   3. Other than the usual philosophical difference of specifying prior knowledge in Bayesian models etc., what are the pros and cons of each approach? Which one would you chose? **Update:** I also found the related question on comparing AIC and BIC. It seems that my method 1 is asymptotically equivalent to AIC and method 2 is asymptotically related to BIC. But I also read there that BIC is equivalent to Leave-One-Out CV. That would mean that the training error minimum and Bayesian Likelihood maximum are equivalent where LOO CV is equivalent to K-fold CV. A perhaps very interesting paper "An asymptotic theory for linear model selection" by Jun Shao relates to these issues.


title: why is the denominator of the correlation coefficient the SD of X multiplied by SD of Y?

text: I don't quite understand what is going on in the correlation coefficient formula. In the numerator we have the covariance, and in the denominator we have the standard deviation of variable x multiplied by the standard deviation of variable y. So ultimately it is a ratio of covariance to the product of the two standard deviations. What is dividing by the product of the two standard deviations doing to help us determine the corrrelation ? I have tried to draw it out to help me visually understand it as I find this helps, but am stuck as to what I should be looking for.


title: SOM automated/objective clustering

text: So as I understand it SOM is primarily a visualization tool and clustering is a logical next step after you construct a SOM from data. Typically, the clustering is subjective in that after looking at your SOM you can 'see' N clusters and then from there you would go on to cluster (k means, hierarchical etc.) the SOM nodes with this N as your parameter. Is there a nonsubjective (maybe nonparametric?) way to cluster the SOM nodes (i.e. without first looking at the SOM to determine the number of distinct clusters)? Thanks


title: How to compute co-occurrence probability assuming independence from unigram probabilities?

text: I have a text corpus, from which I have computed the unigram probabilities of every word. Now, let $p_{i}$ & $p_{j}$ be the unigram probabilities of $i$ and $j$ respectively, what will be the probability with which a pair of words $i$ and $j$ will co-occur(within a maximum gap of 1 word) when independence between $i$ and $j$ is assumed? **EDIT:** Here $p_{i}$ is the probability of any given word being $i$. By co-occurrence probability, I mean the chance that word $i$ & $j$ co-occur in a random sequence of three words drawn from the corpus. By co-occurring I mean that $i$ and $j$ occurring anywhere in this 3 word random sequence within a maximum word distance of 1. (For example, $i$*$j$ and $ij$ are both valid co-occurrences of $i$ & $j$ but $i**j$ is not a valid co-occurrence since the distance between $i$ & $j$ is more than 1.) And finally, $ijij$ counts as 3 co-occurrences.


title: Comparing sample distributions. What should I normalize by for lower population sample distributions?

text: I'm sorry if I'm not describing this accurately. I have a distribution of lengths of a protein. For proteins that are near the average length, I have many more data points, considering there are many more proteins with this average. So I am more confident that the properties of those proteins with a large sample size are more native. However, when you approach the tails of the distribution, you sample size decreases, so the properties of those proteins with large or small lengths may have greater error since the sample size is lower. Can anyone help me with what I should normalize with in order to compare the lengths? I know this is a well known problem in statistics but I'm not sure what it's called.


title: Reference for dimension reduction techniques

text: This is a follow-up question to Is PCA appropriate for comparing subsets of panel data?. It turns out that, yes, PCA is appropriate. But there are also many other ways to reduce n-dimensional data to 1-3 dimensions that can be visualized, and I'm having trouble getting a global view of the available set of techniques for doing so. For instance, several are mentioned in this question and this question. The latter question has a link to Nonlinear Dimensionality Reduction on Wikipedia, which is relatively comprehensive, even encyclopedic, for a Wikipedia article on a technical subject. But it's still Wikipedia and I'd like to have everything in one place. _edit:_ there's also this question from today. **Is there a somewhat-comprehensive reference for dimension reduction?** For instance, I'd really like an explanation of why the "manifold-based" approaches (e.g. ISOMAP) are different or better than distance-based approaches (e.g. classical MDS). Ideally it would include some treatment of different distance metrics as well. Ultimately **I'm trying to get a handle on the kinds of questions each technique is designed to answer**. Obviously this could turn out to be a bigger subject than I realize. I'm aware that, e.g., self-organizing maps are a type of neural network, and neural networks are quite different "under the hood" from, say, multidimensional scaling. I'm hoping that a reference or set of references exists that doesn't require me to know the gory underlying details of each approach. Learning everything there is to know about statistics is a long process and I'm not quite there yet (although I'm making good progress, in part thanks to this site).


title: Patterns in residuals plot from linear regression: do they tell us what model to use?

text: I have carried out a linear regression. This is the form of the model: bounded (0-1) response variable ~ factor1 (2 levels) + factor2 (5 levels) + interaction between factor1:factor2 + factor3 (2 levels) Sample size is about 350. I have plotted the residuals against the predicted values from that model: ![enter image description here](http://i.stack.imgur.com/UJ1Db.jpg) As there is a distinct pattern in the residuals, it appears linear regression is not a suitable model. I have carried out several other linear regressions lately with different data, and I have repeatedly seen this pattern in the residuals. Does this pattern provide information on what model would be best to use? Should an interaction term be added? Is there a predictor missing? Should a non-linear model be used? Or instead, does this particular pattern in the residuals actually not provide any indicators as to what model to use?


title: how to explain the anova object for a multivariate linear regression?

text: I performed a multivariate linear regression such that:               fit<-lm(as.matrix(y)~mwtkg+mbmi+mage,data=x)      where $y$ is a $500 \times 26$ multivariate outcomes. Then, I am wondering how to explain the `anova(fit)`:               > anova(fit)     Analysis of Variance Table                       Df  Pillai approx F num Df den Df    Pr(>F)         (Intercept)   1 0.99959    63064     25    651 < 2.2e-16 ***     mwtkg         1 0.03506        1     25    651    0.5403         mbmi          1 0.20862        7     25    651 < 2.2e-16 ***     mage          1 0.09016        3     25    651 4.567e-05 ***     Residuals   675                                                  ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      What do the three Dfs, Pillai, and P values mean for the model?


title: Model diagnostics for multivariate linear regression

text: I am wondering whether there is a routine way to diagnose a multivariate linear model. For example:               y   <- rmvt(n=100, sigma=diag(3), df=3)  # 3 response vars as spherical multivariate t-dist     x1  <- rchisq(100, 2)                    # x1 distributed as chi-squared     x2  <- rt(100)                           # x2 distributed as t     fit <- lm(as.matrix(y)~x1+x2)            # fits the multivariate regression      How can you diagnose whether the model is fitted well?


title: Ok to cut continuous variable into irregular intervals?

text: I have a predictor variable (water height in metres) and a response variable (feeding rates of birds). The relationship between them looks similar that shown in the plot below.               tidal_effect <- iris[,3:4]     library(plyr); tidal_effect <- rename(tidal_effect, c("Petal.Length" = "feeding.rate", "Petal.Width" =  "water.height"))          tidal_effect_extra <- data.frame(water.height = sample(seq(2.5, 3, 0.1), 50, replace = T), feeding.rate = sample(seq(4, 4.5, 0.1), 50, replace = T))          tidal_effect_extra <- rbind(tidal_effect, tidal_effect_extra)          library(ggplot2)     ggplot(tidal_effect_extra, aes(water.height, feeding.rate)) + geom_point() + xlab("water.height (m)")      ![enter image description here](http://i.stack.imgur.com/w9OyE.jpg) I'm considering cutting water height into a three level factor. As can be seen in plot above, the obvious point at which to cut water height is about 0.7, 2.5 and 3, which means that the distance between the minimum/maximum values at each level of the factor are not equal I can cut the predictor variable like this:               tidal_effect_extra$water.height.cut <- cut(tidal_effect_extra$water.height, breaks = c(0, 0.7, 2.5, 3), labels = c("low", "mid", "high"))      And then perform this model:               lm(feeding.rate ~  water.height.cut, tidal_effect_extra)      My question: is there anything wrong with cutting water height into groups in which distances between the minimum/maximum values at each level of the factor are not equal?


title: Tone mapping a two dimensional distribution

text: I'd like to display some two dimensional distributions, but I'm not sure how to properly tone map them. My distributions have some narrow spikes, so linearly mapping the probabilities to color hides most of the information. I've tried two approaches, doing a gamma adjustment or cutting off high values at some threshold, in both cases I try to target a particular gini distribution for the displayed colors. The results aren't very good though. I feel that the gamma transforms don't really approximate the inverse CDF of my color distribution. I feel reluctant to use a purely local method, because I don't want to lose all absolute contrast difference. Are there other parametric transforms I could try?


title: In a correlation, compensating for one variable being a function of the other

text: I want to correlate a functional outcome measure with pain. However, the functional outcome measure includes pain as one of its 21 questions. Based on how the outcome measure is computed, pain constitutes 5% of its value. Is 5% small enough that a correlation is still valid? The resulting correlation coefficient was 0.71, suggesting that >50% of variability in the outcome measure can be explained by pain. Since 50% is much larger than 5%, can I conclude that the correlation between pain and functional outcome is well above the amount mathematically required by the outcome measure being calculated using pain?


title: How to model prices?

text: I asked this question on the matemathics stackexchange site and was recommended to ask here. I'm working on a hobby project and would need some help with the following problem. ### A bit of context Let's say there is a collection of items with a description of features and a price. Imagine a list of cars and prices. All cars have a list of features, e.g. engine size, color, horse power, model, year etc. For each make, something like this:               Ford:     V8, green, manual, 200hp, 2007, $200     V6, red, automatic, 140hp, 2010, $300     V6, blue, manual, 140hp, 2005, $100     ...      Going even further, the list of cars with prices is published with some time- interval which means we have access to historical price data. Might not always include exactly the same cars. ### Problem I would like to understand how to model prices for any car based on this base information, most importantly cars not in the initial list.               Ford, v6, red, automatic, 130hp, 2009      For the above car, it's almost the same as one in the list, just slightly different in horse power and year. To price this, what is needed? _What I'm looking for is something practical and simple, but I would also like to hear about more complex approaches how to model something like this._ ### What I've tried Here is what I've been experimenting with so far: 1) using historical data to lookup car X. If not found, no price. This is of course very limited and one can only use this in combination with some time decay to alter prices for known cars over time. 2) using a car feature weighting scheme together with a priced sample car. Basically that there is a base price and features just alter that with some factor. Based on this any car's price is derived. The first proved to be not enough and the second proved to not be always correct and I might not have had the best approach to using the weights. This also seems to be a bit heavy on maintaining weights, so that's why I thought maybe there is some way to use the historical data as statistics in some way to get weights or to get something else. I just don't know where to start. ### Other important aspects   * integrate into some software project I have. Either by using existing libraries or writing algorithm myself.   * fast recalculation when new historical data comes in. **Any suggestions how a problem like this could be approached?** All ideas are more than welcome. Thanks a lot in advance and looking forward to reading your suggestions!


title: Binary Logistic Regression: Direction of B's different in multiple than in bivariate cases

text: I am interested in the combined effects of 12 continuous/interval predictors on a dichotomous/binomial outcome and thus am using logistic regression. First I ran 12 separate bivariate logistic regressions. All had positive B's. Then when I ran the multiple logistic regression (enter method), some of these predictors have now become negatively associated with the outcome. I'm not sure what to make of this result. Is this unusual? What does this mean when this happens? Thanks so much for any guidance! Bonnie


title: Problem with unequal distribution of classes in sentiment classification

text: I am performing a binary sentiment classification (positive/negative) with RapidMiner. My problem is that I have about 400 positive and 1350 negative documents. I get pretty good accuracy but therefore my precision and recall for the positive class is around 60-70%. Are there any operators or methods which can help me to get a better recall/precision (f1 score)? I thought "Sample (Bootstrapping)" can help me. But this operator just sizes the whole dataset up or down?


title: Dantzig Selector, LASSO, LAD LASSO

text: I am wondering about this. When is it best to use Dantzig Selector (the infinity normed error measure plus the L1 regularizer) , the LASSO (the mean square error measure plus the L1 regularizer), and the LAD LASSO (the error measure containing the least absolute deviation plus the L1 regularizer)? To make things clear I will quote the minimization problems solved in each: LASSO $$ \min_{\beta} \left| \mathbb{y} - \mathbb{X} \boldsymbol{\beta}\right|_2 + \lambda \left|\boldsymbol{\beta}\right|_1 $$ LAD - LASSO $$ \min_{\beta} \left| \mathbb{y} - \mathbb{X} \boldsymbol{\beta}\right|_1 + \lambda \left|\boldsymbol{\beta}\right|_1 $$ Dantzig Selector $$ \min_{\beta} \left| \mathbb{X}^T( \mathbb{y} - \mathbb{X} \boldsymbol{\beta})\right|_{\infty} + \lambda \left|\boldsymbol{\beta}\right|_1 $$ Here $\mathbb{y}$ represents the observations, $\mathbb{X}$ represents the basis functions for each observation, $\boldsymbol{\beta}$ are the regression coefficients. The subscripts indicate the choice of the norms used. $\lambda$ is a positive coefficient determined by $n$ fold cross-validation typically. If there exists a criterion can this be determined from the dataset? Edit: I think a more reasonable question is the following: What are the features of the fit obtained from each one of these methods. For example I am looking for something like, the LAD-LASSO is less sensitive to outliers in contrast to the LASSO. http://hansheng.gsm.pku.edu.cn/pdf/2007/LAD-Lasso.pdf Also please provide references if you can.


title: "Leave one object out" cross validation

text: I have a training set with 140 instances and no separate cross validation set. The data set contains 7 measurements from each of 20 objects, hence the 140 instances. Each of 7 measurements have the same label. My problem is the choice of cross validation (CV) scheme. I wanted to do 20-fold CV. But this will almost always be biased since the random subsets will have different object measurements mixed. Instead I decided to "leave one object out" each time, i.e., at each run of 20-fold CV, a validation set contains one of the 7 instances and and the remaining 133 will be the training set. Do you think this is a fair CV scheme?


title: Comparing strength of trend and seasonality for various levels of aggregation

text: I have a monthly time-series for which I've fitted the time-series regression model including linear trend and 11 dummy variables for seasonality. Then, I aggregated the original time-series to quarterly level. I then again fit the same model for quarterly data with linear tend and three dummy variables. Question is: **How to compare the significance of estimated trends in both models as well as seasonality?** Recall that response data in both models are different, size of data is different, and the number of parameters as well.


title: Error using ltm R package

text: I'm using the R package ltm to create a 2-parameter logistic regression. The input matrix is sparse - many users have taken a small subset of the items in the item bank. For some of my data sets i'm running into this error:               Error in if (any(ind <- pr == 0)) pr[ind] <- sqrt(.Machine$double.eps) :        missing value where TRUE/FALSE needed      Not sure what the issue is. Doesn't repro on most of my data sets. Any thoughts?


title: Minimum tickets required for specified probability of winning lottery

text: In a lottery, 1/10 of the 50 000 000 tickets give a prize. What is the minimum amount of tickets one should buy to have at least a 50% chance to win? Would be very glad if you could explain your methodology when resolving this. Please consider this as homework if you will.


title: Restricted Boltzmann Machines and Markov Networks: relationship in inference?

text: I am wondering if there is some equivalence between retricted Boltzmann machines and pairwise Markov networks in terms of MAP inference. More specifically, let $y \in \\{0,1\\}^m$ be the output/visible layer and $h\in \\{0,1\\}^k$ be the hidden layer. Specifically $k < n$ (may be $\ll$.) The MAP estimate over RBMs can be written as \begin{equation} \text{arg}\max_{y \in \\{0,1\\}^n ,h \in \\{0,1\\}^k} y^T A h + u^t y + v^t h ~~~(1). \end{equation} for some given $A, y,$ and $v$. Now, a consider a pairwise Markov network defined only over $y$ where the inference is given by \begin{equation} \text{arg}\max_{y \in \\{0,1\\}^n} y^T B y + w^t y ~~~(2), \end{equation} for some given $B$ and $w$. My question is that, given an RBM with matrix $A$ and vectors $u$ and $v$, do there exist $B$ and $w$ such that (1) and (2) are maximized for the same $y$? If not generally, is this true under some conditions? Is the vice-versa (i.e. given $B$ and $w$, come up with equivalent RBM) true with $k < n$?


title: Transforming Percent Change: Lognormal Distribution?

text: I am dealing with a dependent variable relating to the price of an asset and the percent change in the price of that asset. The problem with this type of analysis, like when stocks are the unit of analysis, is that it follows what looks like a log-normal distribution nut has an interval of [-1, +inf). I have never dealt with transforming this type of DV for an OLS analysis. Everything I find in terms of literature is about stocks and why they follow a log-normal and not how to deal with it for any type of linear regression. I would also be open to another type of analysis if that is more 'normal' (no pun intended) and better. An example: If there is an asset X that is worth 10 when it is bought (t = 0) and tomorrow (t = 1) it is valued at 0.01 then this is approximately a 100% decrease in price, and it can basically go no lower. However, if instead it increased to a worth of 1,000 (at t = 1) this is a 9900% increase. Because my data fits a pattern like an asset I have a long right-tail and a distribution that looks log-normal. I was wondering how to deal with it.


title: Whether to leave the data unaltered in the face of outliers and non-normality when performing structural equation modelling?

text: I recently received this email from a graduate student, and I get similar questions often enough, that I thought I'd post it here: > I'm using factor analysis, multiple regression, and SEM and currently > checking statistical assumptions. I have found numerous univariate and > multivariate outliers. If I deleted them all, it would mean a large chunk > out of my sample size ($N \approx 350$). I also have problems with non- > normality, non-linearity, heteroscedasticity (Multiple regression), and > large standardised residual covariances (SEM). > > I have tried reducing the influence of the outliers (allocating them a value > one unit larger/smaller than the next most extreme non-outlier value), and > transformations (mostly the variables remained skewed and some outliers > remain). When I compare original results with altered data, there is little > effect. Given this, I am wondering whether it would be acceptable to leave > the data as it is? I'm inclined to, particularly because this data is from a > non-clinical population and I have used clinical measures.


title: Training a model on a broad patient population, then testing on a specific subset population

text: This is a question that has been bugging me because I see this practice done in many medical papers. Here's the scenario: you create a prognostic model from a general population of breast cancer patients with sensible exclusion criteria like no previous cancer. Your variables end up being age, HER2 receptor status, and estrogen receptor status. What I see in many papers is the authors validate the model, but then also test it on multiple groups of patients with very specific characteristics. Like testing performance for HER2+ and ER+ patients only, or only on young patients who are Stage III. They then say things like, "the model performed well on ER+/HER2- patients but not HER2+ patients overall." Something in my gut says that it is not great practice to train a model on a broad patient population then test it on multiple different subsets of patients (e.g. only Stage III patients, or only patients <50 yrs who are Black). I am not a statistician by training so I have no evidence for this. Anyone know if there are bias issues with this practice?


title: Can the outcome variable in a forward stepwise regression have negative numbers?

text: I'll be the first to admit I'm probably in over my head but...trying to perform a forward stepwise regression (using SPSS 20). Have 7 predictor variables (factors comprising items from a survey) and one outcome variable that comprises values from -43.07 to +25.45 where 0.0 is the "best place to be." ENG is the outcome (so trying to find the factors which best predict engagement) but within the engagement model I'm using, engagement is the result of a balance between challenge and capacity so the difference between challenge and capacity gives an engagement score "of sorts"; perfect alignment would be 0.0; more challenge than capacity = positive numbers; less challenge than capacity = negative numbers. Directionality is important within the model as movement away from engagement (0.0) towards positive values results in people being overwhelmed (more challenge than capacity); the other way is underutilized (less challenge than capacity). Hope this makes sense.


title: Different confidence interval results in R, why?

text: Using the following population:               sample <- c(41.5, 56.7, 54.2, 98.9, 56.7, 43.9, 35.8, 28.8)      I get a different result for the upper and lower confidence intervals when calculating them "manually" than what result from the standard library `t.test()` function.               s <- sd(sample) # standard deviation      se <- s/sqrt(NROW(sample)) # standard error           # Using t.test()     lower <- (t.test(sample))$conf.int[1] # yields 34.11755     upper <- (t.test(sample))$conf.int[2] # yields 70.00745          # Calculating manually     lower <- mean(sample)-(1.96*se) # yields 37.18821     upper <- mean(sample)+(1.96*se) # yields 66.93679      Can somebody explain what is going on here? **Update:** Thanks for the information everybody! This was really enlightening.


title: Using neural network trading in stock exchange, part 2

text: In several months, I have developed a framework for using neural networks (FANN library) in a chart trading software. The framework allows me to combine any inputs for the NN, choose a learning rate, momentum, activation function's variables, number of layers and number of neurons in them...almost anything FANN library offers. The following is the UI. (tl;dr version below) ![UI of the framework](http://i.stack.imgur.com/lcoSJ.png) Basically, I don't know what to do with it now. The internet has been so far my only source of information and I'd appreciate ANYTHING since I'm pretty much stuck right now and can't move anywhere (the project is a requirement for my master's degree). I have tried different setups and right now, I'm just aimlessly changing the inputs and the parameters of the NN. The problem is, everything I do and then start with teaching, the NN is just approximating itself to zero. I have tried incremental learning, batch learning, various values of learning rate and momentum. As an output with which I'm teaching the network, I'm using the slope of the linear regression curve for the following X bars (10 in the picture in the field Output Window). I've been using all sorts of things to the input. The picture shows two CCI oscilators - they are both having 10 as the input window which means there are total 20 inputs to the neural network, the difference between them is the length of the inner simple moving average (15 and 50) and they are both scaled by the maximum and minimum of the whole data set to the range -1, 1 since I'm using the bipolar sigmoid as the activation function. But really, I can insert the price itself, the volume or the day/hour, any average or any derivative such as the CCI as the input and choose how it should be scaled (by the min/max of the whole dataset or by the min/max of the input window). I have tried simple things like choosing only the close prices of the last 10 bars or incredibly complicated things like three CCIs with various lengths, two exponential moving averages, hour, day, last volume change.....there were over 50 inputs, I've used different hidden layer setups... ...but in the end, the neural network is just approximating to zero. I don't know what that means. I thought that the neural network is simply saying "Don't bother with that, there is no relation between the inputs you use and the linear regression curve slope." but then I think it cannot be true - I can see there really are double tops and double bottoms, I can see that there really are some patterns combining the two CCIs that are working repeatedly (patterns described in a book for intraday trading begginers). I'm using 3 months of trading data of the S&P 500 e-mini grouped by 1000 ticks. I was thinking that maybe 3 months is too much, I should probably try just one month and, well, the approximation to zero is still occuring but it is only taking a longer time. Only by using a really limited dataset (like just a few days), the neural network is actually trying to fit itself but I can clearly tell that it is being overtrained for such a short timespan. I have been studying the neural networks a lot, yet I feel I have made some false assumptions about how they should or should not be used. So far, I thought that using various inputs and the linear regression slope as the output will lead the NN to find a relation between the two. I'd like the NN to exaggerate even the slightiest relation found and I thought that might be done using the high momentum and high learning rate...NN is still approximating to zero :( Maybe NN is not supposed to work like that at all ... I've learned a lot about them yet I know very little to come to any conclusions. If you have anything you think could help me (apart from leaving the project), I'd be very grateful. **TL;DR** : My neural network is approximating to zero for different sorts of inputs with diferrent settings, do you know why? Also, any ideas, advices, interesting sources on the topic of neural network usage for trading? Or do you know, at least, where else to ask? **EDIT - problem found, see my answer below**


title: Distribution of 'unmixed' parts based on order of the mix

text: Suppose I have paired observations drawn i.i.d. as $X_i \sim \mathcal{N}\left(0,\sigma_x^2\right), Y_i \sim \mathcal{N}\left(0,\sigma_y^2\right),$ for $i=1,2,\ldots,n$. Let $Z_i = X_i + Y_i,$ and denote by $Z_{i_j}$ the $j$th largest observed value of $Z$. What is the (conditional) distribution of $X_{i_j}$? (or equivalently, that of $Y_{i_j}$) That is, what is the distribution of $X_i$ conditional on $Z_i$ being the $j$th largest of $n$ observed values of $Z$? I am guessing that as $\rho = \frac{\sigma_x}{\sigma_y} \to 0$, the distribution of $X_{i_j}$ converges to just the unconditional distribution of $X$, while as $\rho \to \infty$, the distribution of $X_{i_j}$ converges to the unconditional distribution of the $j$th order statistic of $X$. In the middle, though, I am uncertain.


title: Should I take the arithmetic mean of variates which vary on smaller time scales in a linear regression?

text: I am consulting with another student with a stats background, but am pretty sure he his saying something that contradicts what I have been told about my data. I have a daily response variable measured on 16 individuals in a GLMM where the random effect is the individual. Over 6000 days there are 6000 values. I have several co variates only one of which is measured on the same time scale, daily temperature. The other variates are either vary on a monthly scale, or monthly x individual scale. I have always worked on the presumption that GLMM is powerful enough to handle this, but this guy who is advising me says I have to average my response to a monthly basis, (same for the temperature co variate) so that they vary on same temporal scale. But I am sure that this is unnecessarily throwing away information that is important. When i search for problems relating to scale in stats I get very different types of answers to do with relative values of variates so even if you might shed some light on what to call my problem and where I could find a reference that resolves my query either way.


title: Using LaTeX expression in gnuplot

text: I am having difficulty to get $\bar u$ (to denote average of values $u_i$) as label for the x-axis in the output eps figure. Any help would be appreciated! Many thanks!


title: Bayes rule and base rate

text: My homework question: > An inspector suspects that the food in the factory she is inspecting has > been contaminated with a harmful chemical c. Such chemical contamination > occurs in 5% of factories producing this food. The inspector has a test A > for the chemical which registers positive with 100% certainty when the > chemical is present, but the test also registers positive in 10% of cases > where the chemical is not present. She decides to use this test to help her > decide whether there is contamination. > >   1. Assume that the prior probability of contamination is equal to the base > rate, and that the inspector’s test shows a positive result. Compute the > posterior probability of contamination. >   2. The inspector has another test B for chemical c which only registers > positive 50% of the time when c is present, but has the advantage of never > giving a false positive (i.e., if c is not present, the test will never say > it is). The results of the two tests, A and B, are independent given the > presence or absence of c. It turns out that when the inspector uses test B, > the results are negative. In addition, the inspector knows that the factory > is poorly maintained. The rate of contamination in factories with poor > maintenance is twice as high as the rate in factories overall. Compute the > posterior probability of contamination. > I know these should be rather basic, but I'm getting stuck. For #1 I've reached an answer, but I'm not sure it's correct: Using Bayes rule the probability should be $P(c|positive) = \frac{P(positive|c)P(c)}{P(positive)}$ Now I think that $P(positive)$ should be: $P(positive) = P(positive,c)+P(positive,\neg c) = P(positive|c)P(c) + P(positive|\neg c)P(\neg c)$ Thus: $P(c|positive) = \frac{P(positive|c)P(c)}{P(positive|c) P(c) + P(positive|\neg c) P(\neg c)}$ $P(c|positive) = \frac{1 * 0.05}{1 * 0.05 + 0.1 * 0.95} = 0.34$ **Is this correct?** In part #2 I thought since they are supposed to be independent this must hold: $ P(A_{positive} \cap B_{negative}|c) = P(A_{positive}|c)*P(B_{negative}|c) $ First I adjusted the base rate and recalculated #1: $P(c|A_{positive}) = \frac{1 * 0.1}{1 * 0.1 + 0.1 * 0.9} = 0.52$ The questions says: "the results are negative". So B should be: $P(c|B_{negative}) = \frac{P(B_{negative}|c)P(c)}{P(B_{negative}|c) P(c) + P(B_{negative}|\neg c) P(\neg c)} = \frac{0.5 * 0.1}{0.5 * 0.1 + 1 * 0.9} = 0.05$ $ P(c|A_{positive} \cap B_{negative}) = \frac{P(A_{positive} \cap B_{negative}|c) * P(c)}{P(A_{positive} \cap B_{negative})} = \frac{P(A_{positive}|c) * P(B_{negative}|c) * P(c)}{P(A_{positive} \cap B_{negative})} $ **What should I do next?**


title: Data mining: How should I go about finding the functional form?

text: I'm curious about repeatable procedures that can be used to discover the functional form of the function `y = f(A, B, C) + error_term` where my only input is a set of observations (`y`, `A`, `B` and `C`). Please note that the functional form of `f`is unknown. Consider the following dataset:               AA  BB  CC  DD  EE  FF     ==  ==  ==  ==  ==  ==     98  11  66  84  67  10500     71  44  48  12  47  7250     54  28  90  73  95  5463     34  95  15  45  75  2581     56  37  0   79  43  3221     68  79  1   65  9   4721     53  2   90  10  18  3095     38  75  41  97  40  4558     29  99  46  28  96  5336     22  63  27  43  4   2196     4   5   89  78  39  492     10  28  39  59  64  1178     11  59  56  25  5   3418     10  4   79  98  24  431     86  36  84  14  67  10526     80  46  29  96  7   7793     67  71  12  43  3   5411     14  63  2   9   52  368     99  62  56  81  26  13334     56  4   72  65  33  3495     51  40  62  11  52  5178     29  77  80  2   54  7001     42  32  4   17  72  1926     44  45  30  25  5   3360     6   3   65  16  87  288      In this example, assume that we know that `FF = f(AA, BB, CC, DD, EE) + error term`, but we're not sure about the functional form of `f(...)`. What procedure/what methods would you use to arrive at finding the functional form of `f(...)`? (Bonus point: What is your best guess at the definiton of `f` given the data above? :-) And yes, there is a "correct" answer that will yield an `R^2` in excess of 0.99.)


title: How much undersampling should be done?

text: My aim is to predict quarterly customer-default probabilities: I have data on ~ 2 million individuals, who default on average with a probability of ~ 0.3 percent. Therefore I am thinking about undersampling the majority class (non-defaults) to save computing-time (kernel methods can be quite costly, I know about correcting the probabability predictions) -- the other option being just taking a sub-sample of the data. What do you think would be a good ratio of defaults to non-defaults in my learning sample? Thanks for your help!


title: Can one use Cohen's Kappa for two judgements only?

text: I am using Cohen's Kappa to calculate the inter-agreement between two judges. It is calculated as: $ \frac{P(A) - P(E)}{1 - P(E)} $ where $P(A)$ is the proportion of agreement and $P(E)$ the probability of agreement by chance. Now for the following dataset, I get the expected results:               User A judgements:        - 1, true       - 2, false     User B judgements:        - 1, false       - 2, false     Proportion agreed: 0.5     Agreement by chance: 0.625     Kappa for User A and B: -0.3333333333333333      We can see that both judges have not agreed very well. However in the following case where both judges evaluate one criteria, kappa evaluates to zero:               User A judgements:        - 1, false     User B judgements:        - 1, false     Proportion agreed: 1.0     Agreement by chance: 1.0     Kappa for User A and B: 0      Now I can see that the agreement by chance is obviously 1, which leads to kappa being zero, but does this count as a reliable result? The problem is that I normally don't have more than two judgements per criteria, so these will all never evaluate to any kappa greater than 0, which I think is not very representative. Am I right with my calculations? Can I use a different method to calculate inter-agreement? Here we can see that kappa works fine for multiple judgements:               User A judgements:        - 1, false       - 2, true       - 3, false       - 4, false       - 5, true     User A judgements:        - 1, true       - 2, true       - 3, false       - 4, true       - 5, false     Proportion agreed: 0.4     Agreement by chance: 0.5     Kappa for User A and B: -0.19999999999999996


title: Regression with Pooled data in SPSS

text: ### Context: I have "pooled data" with time and cross section dimensions. It is unbalanced data without a full range of time observations for each cross section of observations. It is kind of like unbalanced panel data? I wish to proceed with the limited number of observations without losing any data. ### Questions:   * How can I get regression results (generalized least squares?) in SPSS with this kind of data?    * What should I pay attention to when evaluating the results?


title: How to compute CDF probability of normal distribution

text: > **Possible Duplicate:**   >  Evaluate definite interval of normal distribution Title was changed and question edited bellow. How is possible that a probability density function defined as following can return some values that are higher than 1 ?! Is there any error in my definition of this pdf ?               float gauss_pdf = exp( -pow(x-mean, 2.0) / (2.0 * pow(sigma, 2.0)) )   /   ( sigma * sqrt(2*M_PI) );      For example, for x = 0.0908182, mean = 0.0552096 and sigma = 0.0241953, it will return 5.583 **EDIT:** Well, what I need is the corresponding probability, thus I need to compute the CDF (i.e. P(X < x)) for that mean and sigma. How can I compute the CDF ?


title: ROC vs precision-and-recall curves

text: I understand the formal differences between them, what I want to know is when it is more relevant to use one vs. the other.   * Do they always provide complementary insight about the performance of a given classification/detection system?    * Is it reasonable to provide them both, say, in a paper? (when does it make sense to provide both of them instead of just one of them?).   * Are there any alternative (maybe more modern) descriptors that capture the relevant aspects of both ROC and precision recall for a classification system? I am interested in arguments for both binary and multi-class (e.g. as one-vs- all) cases.


title: How to describe a design with a mix of experimental conditions, predictor variables, and multiple outcome variables?

text: My research design is as follows: **I have these Between Subjects IVs:**   * Experiment Condition - 5 levels   * Facebook User status - 2 levels (yes/no) **My supervisor also wants me to see if there are significant effects of:**   * Gender - 2 levels   * Relationship Status - 2 levels   * Relationship Satisfaction - 2 levels **And these within subjects IVs:**   * Mood at Time 1   * Mood at Time 2 **I also have three DVs**   * Attraction level   * Frequency of Thought   * Mood State I am looking at the effect of a certain condition on Facebook users, vs non Facebook users. Within this I want to look at whether there are differences between the genders; whether there is a difference between people in a relationship vs. not in a relationship and the effect of relationship satisfaction. It's a big study, but I'm not sure about how to describe the design. ### Question   * How should I describe the design of this study? ### Initial thoughts I'm thinking it has to be a MANOVA design, however what I am finding confusing are the various IVs. Is it a 2 x 5 Factorial MANOVA? I'm just so confused with all the IVs flying around.


title: Representing missing data in an OpenOffice spreadsheet

text: I'm doing a time series analysis. I'm doing most of my analysis in R, where I can use "NA" to represent "not available" (e.g. a missing data point). But I'm doing some data preparation in OpenOffice; currently, I'm leaving cells blank for missing data. Is there a better way to "declare" that a cell is NA?


title: GLMM output interpretation (correct text)

text: I used the `lmer` function in the `lme4` package in order to assess the effects of 2 categorical fixed effects (1º Animal Group: rodents and ants; 2º Microhabitat: bare soil and under cover) on seed predation (a count dependent variable). I have 2 Sites, with 10 trees per site and 4 seed stations per tree. Site and Tree are my (philosophically) random factors, but given that I have only two level for Site, it must be treated as a fixed factor. I have questions about how to interpret the results:   1. I made a model selection criterion based on QAICc, but the best model (lower QAICc) does not result in any significant fixed effect and other models with higher QAIC (e.g. the Full Model) did find significant fixed factors. Does this make sense?    2. Given a fixed factor that is important to the model, how do I distinguish which level of fixed factor is influencing the response variable?  Finally, correlation between the fixed factors implies an incorrect estimation of the model?               FullModel=lmer(SeedPredation ~ AnimalGroup*Microhabitat*Site + (1|Site:Tree) +                                         (1|obs), data=datos,  family="poisson")           QAICc(FM)104.9896              enterGeneralized linear mixed model fit by the Laplace approximation      Formula: SP ~ AG * MH * Site + (1 | Site:Tree) + (1 | obs)         Data: datos         AIC   BIC logLik deviance      101.8 125.6  -40.9     81.8     Random effects:      Groups    Name        Variance Std.Dev.      obs       (Intercept) 0.20536  0.45317       Site:Tree (Intercept) 1.19762  1.09436      Number of obs: 80, groups: obs, 80; Site:Tree, 20          Fixed effects:                      Estimate Std. Error z value Pr(>|z|)       (Intercept)       0.01161    0.47608   0.024   0.9805       AGR             -18.97679 3130.76500  -0.006   0.9952       MHUC             -1.60704    0.63626  -2.526   0.0115 *     Site2            -0.91424    0.74506  -1.227   0.2198       AGR:MHUC         19.92369 3130.76508   0.006   0.9949       AGR:Site2         1.02241 4431.84919   0.000   0.9998       MHUC:Site2        1.80029    0.86235   2.088   0.0368 *     AGR:MHUC:Site2   -3.49042 4431.84933  -0.001   0.9994       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Correlation of Fixed Effects:                 (Intr) AGR    MHUC   Site2  AGR:MHUC AGR:S2 MHUC:S     AGR          0.000                                                 MHUC        -0.281  0.000                                          Site2       -0.639  0.000  0.180                                   AGR:MHUC     0.000 -1.000  0.000  0.000                            AGR:Site2    0.000 -0.706  0.000  0.000  0.706                     MHUC:Site2   0.208  0.000 -0.738 -0.419  0.000    0.000            AGR:MHUC:S2  0.000  0.706  0.000  0.000 -0.706   -1.000  0.000 code here          BestModel=lmer(SP ~ AG * MH + (1|Site:Tree) + (1|obs), data=datos,                      family = "poisson")           QAICc(M) 101.4419          Generalized linear mixed model fit by the Laplace approximation      Formula: SP ~ AG + AG:MH + (1 | Site:Tree) + (1 | obs)         Data: datos         AIC   BIC logLik deviance      100.3 114.6 -44.15     88.3     Random effects:      Groups    Name        Variance Std.Dev.      obs       (Intercept) 0.76027  0.87194       Site:Tree (Intercept) 1.14358  1.06938      Number of obs: 80, groups: obs, 80; Site:Tree, 20          Fixed effects:                  Estimate Std. Error z value Pr(>|z|)     (Intercept)   -0.5153     0.4061  -1.269    0.205     AGR          -18.7146  2603.4397  -0.007    0.994     AGA:MHUC      -0.7301     0.5045  -1.447    0.148     AGR:MHUC      17.7221  2603.4397   0.007    0.995


title: Is it possible to prove a null hypothesis?

text: As the question states - Is it possible to prove the null hypothesis? From my (limited) understanding of hypothesis, the answer is no but I can't come up with a rigorous explanation for it. Does the question have a definitive answer?


title: how to predict Yn value in this formula with Metropolis Hastings or Gibbs?

text: I have a model with this formula: $$ Y_n=aX_n^b + e_n $$ $$ X_n \in [0,2] \quad\quad a = 1.5 \quad\quad b = 0.5 \quad \quad e_n = N(μ = 0, σ^2 = 1) $$ I want to predict "$Y_n$" value with using "Metropolis Hastings" algorithm or "Gibbs" algorithm in these situations: (1) if N = 1 and $X_1$ = 1 then $Y_1$ = 0.5 (2) if N = 1 and $X_1$ = 2 then $Y_1$ = 1.2 (3) if N = 2 and $X_1$ = 1 and $X_2$ = 2 then $Y_1$ = 0.5 and $Y_2$ = 1.2 Why did you choose "Metropolis Hastings" instead of "Gibbs" and vice versa ? Your reasons ? Advantages ? Disadvantages ?


title: How to choose the test statistic for permutation test?

text: Is the purpose of permutation test to test the null that several groups of samples come from the same distribution? I found its steps are > The steps in a permutation-based computation of the significance level of a > test statistic are as follows: > > i) Choose a test statistic, eg. a t-score for a comparison of two groups, > > ii) Compute the test statistic for the gene of interest, > > iii) Permute the labels on samples at random, and re-compute the test > statistic for the rearranged labels; repeat for a large number (perhaps > 1,000) permutations, and finally, > > iv) Compute the fraction of cases in which the test statistics from iii) > exceed the real test statistic from ii). What kinds of test statistic should one choose in the first step? The example uses the t-score, which measures the difference between two groups. But it seems to me that any statistic will work, not necessarily measuring the difference between two groups. Is it correct? Thanks and regards!


title: One-way independent ANOVA where I expect two levels to be very similar and the other to be different. What's the right test?

text: I'm expecting Groups 1 & 2 to be very similar and Group 3 to be different. I'm using a very large data set and feel that power won't be an issue. Initially I felt that I'd then be expecting a significant ANOVA, and then in planned comparisons I'd expect 1 vs 2 to be non-significant while 1 vs 3 and 2 vs 3 should be significant. While that is actually how my data turned out I'm now thinking this was the wrong way to go about things. When comparing 1 vs 2 should I have used equivalence testing (e.g. http://www.cscu.cornell.edu/news/statnews/stnews85.pdf) instead? If so, how could I have gone about that using SPSS or some other commonly available software?


title: Generating random numbers based on partial correlation data

text: I need to generate random numbers based on already existing partial correlation data (not correlation or covariance data). Specifically, a 168*12 matrix based on a 12*12 partial correlation matrix. The idea is to simulate a data matrix that can be used for testing a few components of a project. Any help in this regard would be appreciated. I have looked around but have not found any threads that talk about doing this with partial correlation data. If someone has ideas about implementation in MATLAB, that would be a bonus! Thanks a lot in advance! Additions: Apologies for any ambiguity. -What I mean by partial correlation matrix is a matrix containing the partial correlations, calculated for any two pairs by partialling out effect of all other pairs. -The goal is: given a matrix of partial correlation values, is there a way I can generate a data set (168*12) that would have these partial correlation values? -If there is a method to convert partial correlation to correlation values, that would be appreciated as well. Thanks again!


title: How to find the effective bandwidth correctly using statistics?

text: I have done an experiment to find the effective n/w bandwidth. The data I got in kbps is 223, 221, 510, 220, 471, 229, 222, 221, 220, 221 How can I find the effective bandwidth? Averaging gives 275.8. But if I have done only first 4 rounds then the average is 293.5. How can I find out a more reasonable value as the effective bandwidth. Or is averaging the correct way of doing this?


title: Why do we calculate Information value?

text: I have the data with categorical variables and continuous variables, but is the need for finding information value in explanatory data analysis. Just give the reason for why we are calculating the information value for each variables at the beginning of the data analysis and what will be the cutoff point of INFORMATION VALUE for taking in care of the analysis


title: Test accuracy higher than training. How to interpret?

text: I've a dataset containing at most 150 examples (split into training & test), with many features (higher than 1000). I need to compare classifiers and feature selection methods which perform well on data. So, I'm using three classification methods (J48, NB, SVM) and 2 feature selection methods (CFS, WrapperSubset) with different search methods (Greedy, BestFirst). While comparing, I'm looking at training accuracy (5-fold cross-folding) and test accuracy. Here is one of the results of J48 and CFS-BestFirst: { "accuracyTraining" : 95.83, "accuracyTest" : 98.21 } Many results are like this, and on the SVM there are many results that indicate that test accuracy is much higher than training (training: 60%, test: 98%) How can I meaningfully interpret these kind of results? If it was lower, I would say it's overfitting. Is there something to be said about bias and variance in this case by looking all the results? What can I do to make this classification meaningful, such as re-selecting training and test sets or just using cross-validation on all data? I have 73 training & 58 test instances. Some answers didn't have this info when they were posted.


title: How do you prepare longitudinal data for survival analysis?

text: I'm trying to plot a Kaplan-Meier curve of my data with R. Currently, the data is in the following format:               patient_id;number_of_days;survival     1         ;100           ;T     1         ;200           ;F     1         ;300           ;F     2         ;50            ;F     ...      Survival in my context is not to be interpreted literally: it means progression of disease, which is only assessed at more or less equally spaced time-points. I'm having trouble determining what the `Surv` function expects as input? How am I supposed to prepare the data? I've tried feeding my data directly to `Surv`:               Surv(data$number_of_days,data$survival)      However this seems to count each patient multiple times (once per time-point). So my guess is that I must apply some kind of transformation to my data first. I'm fairly certain this problem of converting longitudinal data to "event data" must be fairly common, so how is this called? Is it already implemented somewhere? Thanks!


title: Independence of samples

text: Recently I started doing verification on daily temperature forecasts by different providers. I calculated the absolute forecast errors for each of the providers, and I am planning on performing ANOVA to check if the differences in mean absolute error between the providers is significant. So lets say that my sample consists of 110 maximum daily temperature forecast errors for 5 different providers. The forecasts errors are obtained as the difference between the actual temperature, obtained from the nearest weather reporting station and the forecasts by each of the providers. What puzzles me here is the independence of the sample, since I use same set of actual temperature values to obtain the errors for each of the forecasts. The more I read about independence of samples the more I get confused. Is the independence assumption violated in this case? Moreover, the distribution of the sample is nearly exponential, so I tried with transformation to approximate it to normal distribution. However, the normality significance tests indicate that the distribution is still non- normal, with a p-value of 0.01. How can I determine if the normality distribution assumption is not seriously violated, since I read on many threads that the test can be performed even if the normality assumption is moderately violated?


title: Conditions under which the power curve is not smooth or continuous in Hypothesis testing

text: Not sure if this is a valid question. I was wondering if there is are good answers to this question: > Under what conditions is the power curve not smooth or continuous in > Hypothesis testing? Thanks.


title: Interpreting the mathematical formula of a mixed effect model

text: I am a bit confused about the function of a parameter in setting up a linear mixed effect model (hierarchical/multilevel model). This is how I understand a (random intercept and slope) multilevel model is set up: Level 1: $$ y_{ij}=\beta_{0j}+\beta_{1j}x_{ij}+e_{ij} $$ (The $i$-th observation in group $j$ is estimated by an intercept for group $j$, a slope for group $j$ of the (here first and only) predictor and an error term for group $j$.) Level 2: `intercept per group:` $\beta_{0j}=\beta_0+u_{0j}$ (The intercept for group $j$ is estimated by a global intercept and a random variable for group $j$.) `slope per group:` $\beta_{1j}=\beta_1+u_{1j}$ (The slope for group $j$ is estimated by a global intercept for the first predictor and a random variable for the $j$-th group of the first predictor.) But I have encountered another way to write the intercept per group: `intercept for group:` $\beta_{0j}=\beta_0+\beta_{0ij}x_{0ij}+u_{0j}$ My interpretation of this would be: The intercept for group $j$ is estimated by a global intercept and a slope for the $i$-th observation in group $j$ multiplied with a predictor $x$ for the intercept and a random variable for group $j$. My first question is if I correctly interpret the subscripts of the slope and predictor of the group specific intercept. Here is what I would guess: $\beta_{0ij}$: $0:=$ just there for consistency; $i:=$ observation; $j:=$ group. (1) Why is the variable $x_{0ij}$ and the parameter/coefficient $\beta_{0ij}$ included? (2) Do they actually function as a fixed predictor in estimating the slope? (3) Are really three subscripts needed? (4) I read the claim that the $x_{0ij}$ and $\beta_{0ij}$ are included in order "make the model symmetric". What does this mean?


title: Average of percentile ranks

text: I have a number of objects (e.g. students) and their scores (possibly transformed into percentile ranks) across various events (e.g. tests). I would like to identify the students who consistently achieved high scores (e.g. top 10%) across the events. My first idea was, for each student, to average (or take the median of) his or her percentile ranks across the tests. However, would you recommend this approach? Are there any caveats that I should be aware of? I tried to do a quick web search and haven't seen this done by anyone.


title: When should normalization never be used?

text: Lately, there have been numerous questions about normalization What are some of the situations where you never ever ever should normalize your data, and what are the alternatives?


title: How do I calculate the SE for PPV, NPV, and DOR?

text: I am attempting to calculate the standard error (SE) for the positive predictive value (PPV), negative predictive value (NPV), and diagnostic odds ratio (DOR) that I have obtained using the rates of true positives, false positives, true negatives, and false negatives in a sample. I am able to get 95% CIs but not SE. Thank you!


title: Do autocorrelated residual patterns remain even in models with appropriate correlation structures, & how to select the best models?

text: ### Context _This question uses R, but is about general statistical issues._ I'm analysing the effects of mortality factors (% mortality due to disease and parasitism) on moth population growth rate over time, where larval populations were sampled from 12 sites once a year for 8 years. The population growth rate data displays a clear but irregular cyclical trend over time. The residuals from a simple generalised linear model (growth rate ~ %disease + %parasitism + year) displayed a similarly clear but irregular cyclical trend over time. Therefore, generalised least squares models of the same form were also fitted to the data with appropriate correlation structures to deal with the temporal autocorrelation, e.g. compound symmetry, autoregressive process order 1 and autoregressive moving average correlation structures. Models all contained the same fixed effects, were compared using AIC, and were fitted by REML (to allow comparison of different correlation structures by AIC). I'm using the R package nlme and the gls function. ### Question 1 The GLS models' residuals still display almost identical cyclical patterns when plotted against time. Will such patterns always remain, even in models that accurately account for the autocorrelation structure? I have simulated some simplified but similar data in R below my second question, which shows the issue _based on my current understanding of the methods needed to assess temporally autocorrelated patterns in model residuals_ , which I now know are wrong (see answer). ### Question 2 I have fitted GLS models with all possible plausible correlation structures to my data, but none are actually substantially better fitting than the GLM without any correlation structure: just one GLS model is marginally better (AIC score = 1.8 lower), whilst all the rest have higher AIC values. However, this is only the case when all models are fitted by REML, not ML where GLS models are clearly much better, but I understand from stats books you must only use REML to compare models with different correlation structures and the same fixed effects for reasons I won't detail here. Given the clearly temporally autocorrelated nature of the data, if no model is even moderately better than the simple GLM what is the most appropriate way to decide which model to use for inference, assuming I'm using an appropriate method (I eventually want to use AIC to compare different variable combinations)? ### Q1 'simulation' exploring residual patterns in models with and without appropriate correlation structures Generate simulated response variable with a cyclical effect of 'time', and a positive linear effect of 'x':               time <- 1:50     x <- sample(rep(1:25,each=2),50)     y <- rnorm(50,5,5) + (5 + 15*sin(2*pi*time/25)) + (x/1)      y should display a cyclical trend over 'time' with random variation:               plot(time,y)      And a positive linear relationship with 'x' with random variation:               plot(x,y)      Create a simple linear additive model of "y ~ time + x":               require(nlme)     m1 <- gls(y ~ time + x, method="REML")      The model displays clear cyclical patterns in the residuals when plotted against 'time', as would be expected:               plot(time, m1$residuals)      And what should be a nice, clear lack of any pattern or trend in the residuals when plotted against 'x':               plot(x, m1$residuals)      A simple model of "y ~ time + x" that includes an autoregressive correlation structure of order 1 should fit the data much better than the previous model because of the autocorrelation structure, when assessed using AIC:               m2 <- gls(y ~ time + x, correlation = corAR1(form=~time), method="REML")     AIC(m1,m2)      However, the model should still display almost identically 'temporally' autocorrelated residuals:               plot(time, m2$residuals)      Thank you very much for any advice.


title: Difference in implementation of binary splits in decision trees

text: I am curious about the practical implementation of a binary split in a decision tree - as it relates to levels of a categorical predictor $X{j}$. Specifically, I often will utilize some sort of sampling scheme (e.g. bagging, oversampling etc) when building a predictive model using a decision tree - in order to improve its predictive accuracy and stability. During these sampling routines, it is possible for a categorical variable to be presented to a tree fitting algorithm with less than the complete level set. Say a variable X takes on levels `{A,B,C,D,E}`. In a sample, maybe only levels `{A,B,C,D}` are present. Then, when the resulting tree is used for prediction, the full set may be present. Continuing from this example, say a tree splits on X and sends `{A,B}` to the left and `{C,D}` to the right. I would expect the logic of the binary split to then say, when faced with new data: "If X has value A or B, send to the left, otherwise, send this case to the right". What seems to happen in some implementations is "if X has value A or B, send to the left, if X has value C or D send to the right". When this case takes on value E, the algorithm breaks down. What is the "right" way for a binary split to be handled? It seems the much more robust way is implemented often, but not always (see Rpart below). **Here are a couple examples:** Rpart fails, the others are ok.               #test trees and missing values          summary(solder)     table(solder$PadType)          # create train and validation     set.seed(12345)     t_rows<-sample(1:nrow(solder),size=360, replace=FALSE)     train_solder<-solder[t_rows,]     val_solder<-solder[-t_rows,]          #look at PadType     table(train_solder$PadType)     table(val_solder$PadType)     #set a bunch to missing     levels(train_solder$PadType)[train_solder$PadType %in% c('L8','L9','W4','W9')] <- 'MISSING'               #Fit several trees, may have to play with the parameters to get them to split on the variable          ####RPART     mod_rpart<-rpart(Solder~PadType,data=train_solder)     predict(mod_rpart,val_solder)     #Error in model.frame.default(Terms, newdata, na.action = na.action, xlev = attr(object,  :      #factor 'PadType' has new level(s) D6, L6, L7, L8, L9, W4          ####TREE     mod_tree<-tree(Solder~PadType,data=train_solder,split="gini")     predict(mod_tree,val_solder) #works fine          ####ctree     mod_ctree<-ctree(Solder~PadType,data=train_solder,control = ctree_control(mincriterion = 0.05))     predict(mod_ctree,val_solder) #works fine


title: Applying Generalized Linear Model to a data with high collinearity

text: my supervisor has given me a response variable that is a count variable, along with 15 explanatory variables. Since the response is a count, I was thinking of doing either Poisson or Negative Binomial regression on the response. However, after computing correlation table and the VIFs, I realized that some of the explanatory variables are highly correlated. I asked my supervisor if I can get rid of some of the explanatory variables in order to avoid the multicollinearity problem, however he said no. How can I safely impose GLM to a data that has high collinearity (some of the correlation coefficients are over .90)?


title: how do I compare two logistic regression models?

text: I have two logistic regression models, using the same data set and same dependent binary variable but with different sample sizes due to different IV's. How would I go about comparing the two models aside from using a classification matrix?


title: Two-sample test for multivariate normal distributions under the assumption that means are the same

text: Let $\\{x_i\\}_{i=1}^n$ be a sample from a multivariate Gaussian distribution ${\cal N}(0, \Sigma_X)$ and $\\{y_i\\}_{i=1}^m$ be a sample from ${\cal N}(0, \Sigma_Y)$. Are there hypothesis tests for $\Sigma_X = \Sigma_Y$? Pointers to relevant literature would be very appreciated.


title: Project management for remote collaboration in prediction

text: Are there any tools for remote collaboration in prediction or machine learning settings? I am looking for a computing environment that includes appropriate source control, keeps track of how different datasets are being matched to different algorithms, and facilitates blending of various predictions. What tools would be useful for this?


title: Interpreting size and direction of fixed effects in a linear mixed effect model

text: I have some data which, after lots of searching, I concluded would probably benefit from using a linear mixed effects model. I think I have an interesting result here, but I am having a little trouble figuring out how to interpret all of the results. This is what I get from the summary() function in R:               > summary(nonzero.lmer)     Linear mixed model fit by REML      Formula: relative.sents.A ~ relative.sents.B + (1 | id.A) + (1 | abstract)         Data: nonzero          AIC    BIC logLik deviance REMLdev      -698.8 -683.9  354.4   -722.6  -708.8     Random effects:      Groups   Name        Variance   Std.Dev.       id.A     (Intercept) 1.0790e-04 0.0103877      abstract (Intercept) 3.0966e-05 0.0055647      Residual             2.9675e-04 0.0172263     Number of obs: 146, groups: id.A, 97; abstract, 52          Fixed effects:                      Estimate Std. Error t value     (Intercept)      0.017260   0.003046   5.667     relative.sents.B 0.428808   0.080050   5.357          Correlation of Fixed Effects:                 (Intr)     rltv.snts.B -0.742      My question involves the relationship between the dependent variable ("relative.sents.A") and "relative.sents.B" once the random factors are factored out. I gather that the t-value of 5.357 for relative.sents.B should be significant. But does this show what the direction of the effect is? I am thinking that because the coefficient for the slope is positive that this means that as relative.sents.B increases, so does my dependent variable. Is this correct? The book I've been using briefly mentions that the correlation reported here is not a normal correlation, but goes into no details. Normally, I'd look there to figure out the direction and magnitude of the effect. Is that wrong? If I'm wrong on both counts, then what is a good (hopefully reasonably straightforward) way to discover the direction and size of the effect?


title: How to change X-axis format in plot

text: I have fit a linear model and I now doing a lack-of-fit test. If the model has a good fit, the residual deviance should follow a $\chi^2$ distribution. The model I just ran has a very poor fit, and now I want to show this visually. I have used the following R code to plot the value of the residual deviance and the $\chi^2$ distribution it should come from (under the hypothesis of good fit).               df <- 384162     resdev <- 750749     x <- seq(from=350000, to=800000, by=1)     y <- dchisq(x, df=df)     plot(x, y, type='l')     abline(v=resdev, col='red')      Everything is working fine, except now the X-axis values are displaying in scientific notation. How to I control the formating of the axis values? I understand it may be crowded if I use normal notation, so can I format labels something like `350k` for 350,000?


title: What are correct values for precision and recall in edge cases?

text: Precision is defined as:               p = true positives / (true positives + false positives)      Is it correct that, as `true positives` and `false positives` approach 0, the precision approaches 1? Same question for recall:               r = true positives / (true positives + false negatives)      I am currently implementing a statistical test where I need to calculate these values, and sometimes it happens that the denominator is 0, and I am wondering which value to return for this case. P.S.: Excuse the inappropriate tag, I wanted to use `recall`, `precision` and `limit`, but I cannot create new Tags yet.


title: How to visualize what ANOVA does?

text: What way (ways?) is there to visually explain what is ANOVA? Any references, link(s) (R packages?) will be welcomed. Thanks.


title: Explaining PLS-DA to a layman

text: I recently learned about PLS-DA in a statistics class. I am able to perform PLS-DA mathematically, but I am having trouble really explaining it. I was wondering if anyone could help me with how to explain PLS-DA to laymen in non mathematical way?


title: If my counts are from a proportion (resampling) of a sample, does it influence the confidence interval?

text: The aim is to estimate the proportion (or counts) of a certain species among a collected sample. Due to the large number, it is not possible to count through the whole sample. So I have to select a proportion of the sample (roughly divide the sample by size or weight), for instance 1/3 of this sample. My question is: will this extra resampling influence the confidence interval of the estimated proportion for this species? especially if the species is rare (proportion is almost 0)? Is there any analytical interpretation? Thanks


title: Simple way to algorithmically identify a spike in recorded errors

text: We need an early warning system. I am dealing with a server that is known to have performance issues under load. Errors are recorded in a database along with a timestamp. There are some manual intervention steps that can be taken to decrease the server load, but only if someone is aware of the issue... Given a set of times that errors occurred, how can I identify the beginning of a spike in errors (in real time)? We can calculate periodically or on each error occurrence. We are unconcerned about occasional errors, but don't have a specific threshold. I could just notify someone any time we get, say, three errors in five minutes, but I'm sure there's a better way... I'd like to be able to adjust the sensitivity of the algorithm based on feedback from the sysadmins. For now, they'd like it to be fairly sensitive, even though we know we can expect some false positives. I am not a statistician, which I'm sure is obvious, and implementing this needs to be relatively simple with our existing tools: SQL Server and old- school ASP JScript. I'm not looking for an answer in code, but if it requires additional software, it probably won't work for us (though I welcome impractical but ideal solutions as a comment, for my own curiosity).


title: Some questions about confidence intervals, quantiles and cdf/pdf

text: My question is related to : Maximum likelihood estimation and the n-th order statistic estimation-and-the-n-th-order-statistic. I will try to answer the questions, and I would like to ask you guys to look at it, and criticize it, because I want all the answers to be correct obviously :) > Given $$ f_{\theta}(x)=e^{ (\theta -x)}$$for $x \ge \theta$, otherwise > $f_\theta = 0$ > > I showed that the MLE $\hat \theta_n $of $\theta= min(x_1,x_2, ... x_n) = > X_{(1)}$ I have to answer the following questions about this distribution: > a.) What is the cdf of $\hat \theta_n $ > > b.) What is the pdf of $\hat \theta_n $ and draw the distribution > > c.) We take a sample of size $n$ from this distribution and by chance we > took a very unusual value $\hat \theta_n = q_{1-\alpha}$ where > $q_{1-\alpha}$ is the $1-\alpha$ quantile of the distribution of > $\hat\theta_n$. Determine $q_{1-\alpha}$ as a function of of $\theta$ and > $\alpha$ > > d.) Again we take sample of size $n$ and now we got very lucky by obtaining > $\hat\theta_n=q_0$ where $q_0$ is the 0th quantile of the $\hat \theta_n$ > distribution. Determine $q_o$ > > e.)We observe $X_1=3, X_2=4, X_3=10$. With the answers of c.) and d.) we > should now be able to obtain a 95%-confidence interval for $\theta$ My try: a.) the cdf is $F_{\hat\theta_n}1-(1-(F_{\theta}(x))^n= 1-e^{n\theta - nx}$ b.) $f_{\hat\theta_n}= ne^{n\theta - nx} = \frac{d}{dx}( 1-e^{n\theta - nx})$. c.) $$1-\alpha =F(q_{1-\alpha})= 1-e^{n\theta -nq_{1-\alpha}} \implies \ln{(1-\alpha)}=n\theta-nq_{1-\alpha} \implies -nq_{1-\alpha}= \ln{(1-\alpha)}-n\theta \implies q_{1-\alpha}=-\frac{\ln{(1-\alpha)}}{n} +\theta$$ d.)$$0=F(q_0)= 1- e^{n\theta -nq_0} \implies 1=e^{n(\theta-q_0)} \implies q_0=\theta$$ e.) We need a $1-\alpha$-confidence interval for $\theta$. Let $c_l$ and $c_u$ denote the lower and upper boundaries of the confidence interval of $\theta$. I will solve $$P(c_l < \hat\theta_n < c_r)= 1-\alpha $$ by requiring $$P(\hat\theta_n \le c_l) =P(\hat\theta_n \ge c_r)= \frac{1}2{}\alpha $$ Now I find the following equations $$1-e^{n\theta-nc_l} = \frac{\alpha}{2}\implies \ln(1-\frac{1}{2}\alpha)=n\theta-nc_l \implies c_l = -\frac{\ln(1-\frac{1}{2}\alpha)}{n}+\theta $$ Similarly $$e^{n\theta-nc_u}= \frac{1}{2} \alpha \implies c_u=-\frac{\ln(\frac{1}{2}\alpha)}{n}+\theta $$ Now I have found that, whatever the value of $\theta$: $$P(-\frac{\ln(1-\frac{1}{2}\alpha)}{n}+\theta <\hat\theta_n <-\frac{\ln(\frac{1}{2}\alpha)}{n}+\theta ) = 1-\alpha $$ By re-arranging the inequalities, we see this is equivalent to $$ P(\frac{\ln(\frac{1}{2}\alpha)}{n}+\hat\theta_n <\theta < \frac{\ln(1-\frac{1}{2}\alpha)}{n}+\hat\theta_n ) = 1-\alpha$$ Therefore, a $100(1-\alpha) $% confidende interval for $\theta$ is given by: $$\left( \frac{\ln(\frac{1}{2}\alpha)}{n}+\hat\theta_n <\theta < \frac{\ln(1-\frac{1}{2}\alpha)}{n}+\hat\theta_n \right) $$ Looking at the observation, we have min(X1, X2, X3)= 3 and n = 3. So our 95% confidence interval should be: $$\left( \frac{\ln(\frac{0.05}{2})}{3}+3 <\theta < \frac{\ln(1-\frac{0.05}{2})}{3}+3 \right) $$This interval is [1.77, 2.99]


title: How to compare the results of a leave-one-out strategy?

text: I tried to reproduce the experiments described in this paper and wanted to compare the output of my system with the one described in the article. I am looking for a statistical comparison of the systems (I'm not really confident when it comes to hypothesis testing). To be more precise the dataset consists of images of certain categories of rooms (bathroom, living room, etc.) acquired in 6 different houses. In the experiment, we use 5 houses to train a classifiers and the 6th to test it (the author call it leave-one-out but strictly speaking I'm not sure it's exactly that). The results is the correct classification rate for each room in each "fold". Do you have any idea?


title: Given Two 95% Confidence Intervals

text: Suppose we have given a two $95 \%$ confidence intervals for $X_1$ and $X_2$. They are normally distributed. From this how would we get a $95 \%$ confidence interval for $X_{1}/X_{2}$?


title: A simple mathematics?

text: I am struggling on how to put below problem with simple mathematics. Let say I have 2 variables N & p. The 2 instances of N are N1 & N2 and 2 instances of p are p1 and p2. Now please consider below expression: abs(N1p1 - N2p2) (i.e. absolute diff. between 2 instances of Np) Now I need to answer this question: Out of above difference, how much (in percentage) is due to variance for N and how much due to p. Therefore I need to calculate: A1 = abs(N1 - N2) A2 = abs(p1 - p1) and since overall difference in abs(N1p1 - N2p2) is only due to difference between N1, N2 and p1, p2, I should be able to show A1 + A2 = 100% So my question is how to arrive A1 & A2 such that they sum up at 100%? Or is there any other way to answer the question: Out of the difference abs(N1p1 - N2p2), how much (in percentage) is due to variance for N and how much due to p Can someone here please help me out on this regards? Thanks for your help.


title: Bayesian Stats--Normal Distribution: known variance, unknown mean

text: I am given a situation where I have an unknown mean, mu, and a known variance. I have prior information about the mean that it follows a normal distribution and has a particular mean and variance that is given. I do not have sample info. Where do I go from here to find probabilities of my variable of interest?


title: Aligned Rank Transform for mixed effects ANOVA in SPSS

text: In a book, I've found SPSS syntax for performing an aligned rank transformation in a two way ANOVA. The procedure is: \- save residuals by performing a standard ANOVA \- use Aggregate to determine effects for group means (mij for interaction, ai as first factor, bj for second factor) \- eliminate interaction effect from residual to determine interaction (ran), eliminate single effects to determine main effects (ra) \- transform residuals into ranks \- run complete model on these ranks for interaction / main effects, respectively The syntax is:               Unianova  x by group drugs       /save=resid (rab)       /design=group drugs group*drugs.     Compute ra=rab.     Aggregate       /outfile=* mode=addvariables       /break=group drugs    /mij=mean(x).     Aggregate       /outfile=* mode=addvariables       /break=group    /ai=mean(x).     Aggregate       /outfile=* mode=addvariables       /break=drugs    /bj=mean(x).     Aggregate       /outfile=* mode=addvariables       /break=         /mm=mean(x).     Compute rab=rab-(mij - ai - bj + mm).     Compute ra=ra-(ai + bj - 2*mm).     Rank variables=ra rab (A)       /rank into rar rabr.     Unianova  rabr by group drugs       /design=group drugs group*drugs.     Unianova  rar by group drugs       /design=group drugs group*drugs.      Now I would like to adapt this to my design, which has an additional random factor "subjects". My initial ANOVA syntax would look like this:               Unianova x by condition repetition subject      /random=subject      /method=SSTYPE(3)      /intercept=INCLUDE      /save=RESID (rab)        /criteria=ALPHA(.05)      /design=condition repetition subject condition*repetition condition*subject repetition*subject.      However now I'm not sure **how to include the third random effect** in the subsequent aligned rank transform procedure. Any tips are greatly appreciated! Thanks! (Source: http://www.uni-koeln.de/~a0032/statistik/buch/nonpar-anova-buch.pdf, in German)