title: Can the scaling values in a linear discriminant analysis (LDA) be used to plot explanatory variables on the linear discriminants?

text: Using a biplot of values obtained through principal component analysis, it is possible to explore the explanatory variables that make up each principle component. **Is this also possible with Linear Discriminant Analysis?** Examples provided use the The data is "Edgar Anderson's Iris Data" (http://en.wikipedia.org/wiki/Iris_flower_data_set). Here is the _iris data_ :                 id  SLength   SWidth  PLength   PWidth species              1      5.1      3.5      1.4       .2 setosa         2      4.9      3.0      1.4       .2 setosa         3      4.7      3.2      1.3       .2 setosa         4      4.6      3.1      1.5       .2 setosa         5      5.0      3.6      1.4       .2 setosa         6      5.4      3.9      1.7       .4 setosa         7      4.6      3.4      1.4       .3 setosa         8      5.0      3.4      1.5       .2 setosa         9      4.4      2.9      1.4       .2 setosa        10      4.9      3.1      1.5       .1 setosa        11      5.4      3.7      1.5       .2 setosa        12      4.8      3.4      1.6       .2 setosa        13      4.8      3.0      1.4       .1 setosa        14      4.3      3.0      1.1       .1 setosa        15      5.8      4.0      1.2       .2 setosa        16      5.7      4.4      1.5       .4 setosa        17      5.4      3.9      1.3       .4 setosa        18      5.1      3.5      1.4       .3 setosa        19      5.7      3.8      1.7       .3 setosa        20      5.1      3.8      1.5       .3 setosa        21      5.4      3.4      1.7       .2 setosa        22      5.1      3.7      1.5       .4 setosa        23      4.6      3.6      1.0       .2 setosa        24      5.1      3.3      1.7       .5 setosa        25      4.8      3.4      1.9       .2 setosa        26      5.0      3.0      1.6       .2 setosa        27      5.0      3.4      1.6       .4 setosa        28      5.2      3.5      1.5       .2 setosa        29      5.2      3.4      1.4       .2 setosa        30      4.7      3.2      1.6       .2 setosa        31      4.8      3.1      1.6       .2 setosa        32      5.4      3.4      1.5       .4 setosa        33      5.2      4.1      1.5       .1 setosa        34      5.5      4.2      1.4       .2 setosa        35      4.9      3.1      1.5       .2 setosa        36      5.0      3.2      1.2       .2 setosa        37      5.5      3.5      1.3       .2 setosa        38      4.9      3.6      1.4       .1 setosa        39      4.4      3.0      1.3       .2 setosa        40      5.1      3.4      1.5       .2 setosa        41      5.0      3.5      1.3       .3 setosa        42      4.5      2.3      1.3       .3 setosa        43      4.4      3.2      1.3       .2 setosa        44      5.0      3.5      1.6       .6 setosa        45      5.1      3.8      1.9       .4 setosa        46      4.8      3.0      1.4       .3 setosa        47      5.1      3.8      1.6       .2 setosa        48      4.6      3.2      1.4       .2 setosa        49      5.3      3.7      1.5       .2 setosa        50      5.0      3.3      1.4       .2 setosa        51      7.0      3.2      4.7      1.4 versicolor        52      6.4      3.2      4.5      1.5 versicolor        53      6.9      3.1      4.9      1.5 versicolor        54      5.5      2.3      4.0      1.3 versicolor        55      6.5      2.8      4.6      1.5 versicolor        56      5.7      2.8      4.5      1.3 versicolor        57      6.3      3.3      4.7      1.6 versicolor        58      4.9      2.4      3.3      1.0 versicolor        59      6.6      2.9      4.6      1.3 versicolor        60      5.2      2.7      3.9      1.4 versicolor        61      5.0      2.0      3.5      1.0 versicolor        62      5.9      3.0      4.2      1.5 versicolor        63      6.0      2.2      4.0      1.0 versicolor        64      6.1      2.9      4.7      1.4 versicolor        65      5.6      2.9      3.6      1.3 versicolor        66      6.7      3.1      4.4      1.4 versicolor        67      5.6      3.0      4.5      1.5 versicolor        68      5.8      2.7      4.1      1.0 versicolor        69      6.2      2.2      4.5      1.5 versicolor        70      5.6      2.5      3.9      1.1 versicolor        71      5.9      3.2      4.8      1.8 versicolor        72      6.1      2.8      4.0      1.3 versicolor        73      6.3      2.5      4.9      1.5 versicolor        74      6.1      2.8      4.7      1.2 versicolor        75      6.4      2.9      4.3      1.3 versicolor        76      6.6      3.0      4.4      1.4 versicolor        77      6.8      2.8      4.8      1.4 versicolor        78      6.7      3.0      5.0      1.7 versicolor        79      6.0      2.9      4.5      1.5 versicolor        80      5.7      2.6      3.5      1.0 versicolor        81      5.5      2.4      3.8      1.1 versicolor        82      5.5      2.4      3.7      1.0 versicolor        83      5.8      2.7      3.9      1.2 versicolor        84      6.0      2.7      5.1      1.6 versicolor        85      5.4      3.0      4.5      1.5 versicolor        86      6.0      3.4      4.5      1.6 versicolor        87      6.7      3.1      4.7      1.5 versicolor        88      6.3      2.3      4.4      1.3 versicolor        89      5.6      3.0      4.1      1.3 versicolor        90      5.5      2.5      4.0      1.3 versicolor        91      5.5      2.6      4.4      1.2 versicolor        92      6.1      3.0      4.6      1.4 versicolor        93      5.8      2.6      4.0      1.2 versicolor        94      5.0      2.3      3.3      1.0 versicolor        95      5.6      2.7      4.2      1.3 versicolor        96      5.7      3.0      4.2      1.2 versicolor        97      5.7      2.9      4.2      1.3 versicolor        98      6.2      2.9      4.3      1.3 versicolor        99      5.1      2.5      3.0      1.1 versicolor       100      5.7      2.8      4.1      1.3 versicolor       101      6.3      3.3      6.0      2.5 virginica       102      5.8      2.7      5.1      1.9 virginica       103      7.1      3.0      5.9      2.1 virginica       104      6.3      2.9      5.6      1.8 virginica       105      6.5      3.0      5.8      2.2 virginica       106      7.6      3.0      6.6      2.1 virginica       107      4.9      2.5      4.5      1.7 virginica       108      7.3      2.9      6.3      1.8 virginica       109      6.7      2.5      5.8      1.8 virginica       110      7.2      3.6      6.1      2.5 virginica       111      6.5      3.2      5.1      2.0 virginica       112      6.4      2.7      5.3      1.9 virginica       113      6.8      3.0      5.5      2.1 virginica       114      5.7      2.5      5.0      2.0 virginica       115      5.8      2.8      5.1      2.4 virginica       116      6.4      3.2      5.3      2.3 virginica       117      6.5      3.0      5.5      1.8 virginica       118      7.7      3.8      6.7      2.2 virginica       119      7.7      2.6      6.9      2.3 virginica       120      6.0      2.2      5.0      1.5 virginica       121      6.9      3.2      5.7      2.3 virginica       122      5.6      2.8      4.9      2.0 virginica       123      7.7      2.8      6.7      2.0 virginica       124      6.3      2.7      4.9      1.8 virginica       125      6.7      3.3      5.7      2.1 virginica       126      7.2      3.2      6.0      1.8 virginica       127      6.2      2.8      4.8      1.8 virginica       128      6.1      3.0      4.9      1.8 virginica       129      6.4      2.8      5.6      2.1 virginica       130      7.2      3.0      5.8      1.6 virginica       131      7.4      2.8      6.1      1.9 virginica       132      7.9      3.8      6.4      2.0 virginica       133      6.4      2.8      5.6      2.2 virginica       134      6.3      2.8      5.1      1.5 virginica       135      6.1      2.6      5.6      1.4 virginica       136      7.7      3.0      6.1      2.3 virginica       137      6.3      3.4      5.6      2.4 virginica       138      6.4      3.1      5.5      1.8 virginica       139      6.0      3.0      4.8      1.8 virginica       140      6.9      3.1      5.4      2.1 virginica       141      6.7      3.1      5.6      2.4 virginica       142      6.9      3.1      5.1      2.3 virginica       143      5.8      2.7      5.1      1.9 virginica       144      6.8      3.2      5.9      2.3 virginica       145      6.7      3.3      5.7      2.5 virginica       146      6.7      3.0      5.2      2.3 virginica       147      6.3      2.5      5.0      1.9 virginica       148      6.5      3.0      5.2      2.0 virginica       149      6.2      3.4      5.4      2.3 virginica       150      5.9      3.0      5.1      1.8 virginica      Example PCA biplot using the iris data set in R (code below): ![enter image description here](http://i.stack.imgur.com/urswP.jpg) This figure indicates that Petal length and Petal width are important in determining PC1 score and in discriminating between Species groups. setosa has smaller petals and wider sepals. Apparently, similar conclusions can be drawn from plotting linear discriminant analysis results, though I am not certain what the LDA plot presents, hence the question. The axis are the two first linear discriminants (LD1 99% and LD2 1% of trace). The coordinates of the red vectors are "Coefficients of linear discriminants" also described as "scaling" (lda.fit$scaling: a matrix which transforms observations to discriminant functions, normalized so that within groups covariance matrix is spherical). "scaling" is calculated as `diag(1/f1, , p)` and `f1 is sqrt(diag(var(x - group.means[g, ])))`. Data can be projected onto the linear discriminants (using predict.lda) (code below, as demonstrated http://stackoverflow.com/a/17240647/742447). The data and the predictor variables are plotted together so that which species are defined by an increase in which predictor variables can be seen (as is done for usual PCA biplots and the above PCA biplot).: ![Example LDA biplot using the iris data set in R](http://i.stack.imgur.com/5xWcW.jpg) From this plot, Sepal width, Petal Width and Petal Length all contribute to a similar level to LD1. As expected, setosa appears to smaller petals and wider sepals. There is no built-in way to plot such biplots from LDA in R and few discussions of this online, which makes me wary of this approach. **Does this LDA plot (see code below) provide a statistically valid interpretation of predictor variable scaling scores ?** Code for PCA:               require(grid)            iris.pca <- prcomp(iris[,-5])       PC <- iris.pca       x="PC1"       y="PC2"       PCdata <- data.frame(obsnames=iris[,5], PC$x)            datapc <- data.frame(varnames=rownames(PC$rotation), PC$rotation)       mult <- min(         (max(PCdata[,y]) - min(PCdata[,y])/(max(datapc[,y])-min(datapc[,y]))),         (max(PCdata[,x]) - min(PCdata[,x])/(max(datapc[,x])-min(datapc[,x])))       )       datapc <- transform(datapc,                           v1 = 1.6 * mult * (get(x)),                           v2 = 1.6 * mult * (get(y))       )            datapc$length <- with(datapc, sqrt(v1^2+v2^2))           datapc <- datapc[order(-datapc$length),]            p <- qplot(data=data.frame(iris.pca$x),                      main="PCA",                      x=PC1,                      y=PC2,                      shape=iris$Species)       #p <- p + stat_ellipse(aes(group=iris$Species))       p <- p + geom_hline(aes(0), size=.2) + geom_vline(aes(0), size=.2)       p <- p + geom_text(data=datapc,                           aes(x=v1, y=v2,                              label=varnames,                              shape=NULL,                              linetype=NULL,                              alpha=length),                           size = 3, vjust=0.5,                          hjust=0, color="red")       p <- p + geom_segment(data=datapc,                              aes(x=0, y=0, xend=v1,                                 yend=v2, shape=NULL,                                  linetype=NULL,                                 alpha=length),                             arrow=arrow(length=unit(0.2,"cm")),                             alpha=0.5, color="red")       p <- p + coord_flip()                 print(p)      Code for LDA               #Perform LDA analysis     iris.lda <- lda(as.factor(Species)~.,                      data=iris)          #Project data on linear discriminants     iris.lda.values <- predict(iris.lda, iris[,-5])          #Extract scaling for each predictor and     data.lda <- data.frame(varnames=rownames(coef(iris.lda)), coef(iris.lda))          #coef(iris.lda) is equivalent to iris.lda$scaling          data.lda$length <- with(data.lda, sqrt(LD1^2+LD2^2))     scale.para <- 0.75          #Plot the results     p <- qplot(data=data.frame(iris.lda.values$x),                    main="LDA",                    x=LD1,                    y=LD2,                    shape=iris$Species)#+stat_ellipse()     p <- p + geom_hline(aes(0), size=.2) + geom_vline(aes(0), size=.2)     p <- p + theme(legend.position="none")     p <- p + geom_text(data=data.lda,                        aes(x=LD1*scale.para, y=LD2*scale.para,                            label=varnames,                             shape=NULL, linetype=NULL,                            alpha=length),                        size = 3, vjust=0.5,                        hjust=0, color="red")     p <- p + geom_segment(data=data.lda,                           aes(x=0, y=0,                               xend=LD1*scale.para, yend=LD2*scale.para,                               shape=NULL, linetype=NULL,                               alpha=length),                           arrow=arrow(length=unit(0.2,"cm")),                           color="red")     p <- p + coord_flip()          print(p)      The results of the LDA are as follows               lda(as.factor(Species) ~ ., data = iris)          Prior probabilities of groups:         setosa versicolor  virginica       0.3333333  0.3333333  0.3333333           Group means:                Sepal.Length Sepal.Width Petal.Length Petal.Width     setosa            5.006       3.428        1.462       0.246     versicolor        5.936       2.770        4.260       1.326     virginica         6.588       2.974        5.552       2.026          Coefficients of linear discriminants:                         LD1         LD2     Sepal.Length  0.8293776  0.02410215     Sepal.Width   1.5344731  2.16452123     Petal.Length -2.2012117 -0.93192121     Petal.Width  -2.8104603  2.83918785          Proportion of trace:        LD1    LD2      0.9912 0.0088


title: Linear Regression with Outlier accounting in Bugs

text: I'm trying to redo an exercise in BUGS from this webpage: a linear regression over a data set with some outliers, using a model that accounts for them. This model uses a mixture of signal and noise (the outliers). The likelihood is: $$p(x_i,y_i,e_i | \theta, g_i,\sigma_B) = \frac{g_i}{\sqrt{2\pi e_i^2}} \exp \Big\\{ -\frac{1}{2e_i^2} (y - \hat{y}(x|\theta))^2 \Big\\} \\\ + \frac{1-g_i}{\sqrt{2\pi \sigma_B^2}} \exp \Big\\{ -\frac{1}{2 \sigma_B^2} (y - \hat{y}(x|\theta))^2 \Big\\}$$ where ($x_i,y_i$) is the available data, and $e_i$ the respective error, $\hat{y}(x|\theta)$ is the estimate $\theta_0 + \theta_1 x$. The nuisance parameters $g_i$ range from $0$ to $1$ and state the degree to which the i-th data point fits the model (so, $g_i \approx 0$ indicates an outlier) and $\sigma_B$ is the variance of the Gaussian that models the noise (herein, just some big number, say $50$). I've tried to code this in BUGS:               model {       for (i in 1:n) {         tau[i] <- 1/pow(e[i],2)         part1[i] ~ dnorm(mu[i], tau[i])              part2[i] ~ dnorm(mu[i], tauOutliers)              y[i] <- g[i] * part1[i] + (1-g[i]) * part2[i]              mu[i] <- theta0 + theta1 * x[i]         g[i] ~ dunif(0,1)       }            theta0 ~ dflat()       theta1 ~ dflat()       tauOutliers <- 1/sigmaB     }      However I get a `multiple definitions of node y[1]` error. I think BUGS does not accept deterministic assignments to observable data. It also does not accept the more direct `y[i] ~ g[i] * dnorm(mu[i], tau[i]) + (1-g[i]) * dnorm(mu[i], tauOutliers)`. Does anyone know how to solve this problem in BUGS? Thanks,


title: How can I feed the randomized tables into the linear model as its 'expected tables'

text: I want to estimate all the possible interactions of three variables with linear model in R. My dataset is like this:                   a   b   c   Freq     1   0   0   No  457     2   1   0   No  14     3   0   1   No  17     4   1   1   No  2     5   0   0   Yes 450     6   1   0   Yes 16     7   0   1   Yes 25     8   1   1   Yes 1      And I make a linear model in R:               glm.model = glm (Freq ~ A * B * C, data = mytable, family = poisson)      But I have my own expected tables from randomization, so I want to feed my randomization tables into the linear model as its 'expected tables', instead of letting it generate its own expected tables. I'm still pretty new to linear model, so do you know how can I solve this problem in R? Thank you so much :)


title: Find a projection of a $k$-simplex with minimal "radius"

text: I asked the following question on Math.SE and I've been told that I may get more helpful answers here on Cross Validated. This is my first question here, so I hope the tags I've chosen are accurate enough. > Let $S$ be a $k$-simplex in $\mathbb{R}^k$. I'd like to find a hyperplane > $P$ that passes through the origin such that projections of vertices of $S$ > onto $P$ are as close as possible to the origin (in minimax sense). The > distances are measured in $\ell_2$-norm. > > I could only get trivial results for $k=2$, but even for $k=3$ I couldn't > find a way to find a solution. > > Any hint on whether there is a polynomial-time algorithm that can solve this > would be appreciated. If you know papers that are relevant I'd be thankful > if you let me know.


title: Metropolis Ergodicity

text: I have encountered one last problem with regarding to the Metropolis-Hastings algorithm. I know that ergodicity is needed in the algorithm to imply convergence to a unique stationary distribution. But how is ergodicity proved in the M-H algorithm? Thank you for your help


title: Complex-valued estimator of Real Parameter?

text: I was in statistics class today and the professor wrote "consistent estimator". The first thing I thought of was that the range of the estimator should be a subset of the set of plausible population parameters (this was not the definition). Is there a name for estimators with ranges larger than the set of plausible parameters? Does anyone ever use such strange estimators in applications? For example, we could consider (1-a)X+aV where X, V are unbiased estimators of some real parameter. A complex a-value would yield an unbiased estimator with complex outputs! (Note: a usual definition of an estimator would exclude this; this entails a loosening of the definition of an estimator to simply a function from the data sets.)


title: VAR with dummies, how to analyze dummy effect?

text: I'm currently trying to perform a VAR model analysis. I have the following variables, all on a one-minute scale, with 900 entries. Put in eViews as endogenous: Y1: Amount of tweets per minute Y2: TV ratings Put in eViews as exogenous, because they are dummies: X1: Times when social media encouragement was on-screen (0,1) X2: Emotional moment (0,1) I would like to test the effect of X1 on Y1 and in my VAR I get a coefficient of around 150. Which can be correct, according to this graph. graph Now, most papers do not use these coefficients but test the effect of a shock. Though, this is only possible for the endogenous variables Y1 and Y2. I cannot find any research where the effect of a dummy in a VAR is explained. Can anyone help me get into the right direction? Do I use the VAR coefficient or can I do other tests with these dummies? Thanks so much, I'm pretty much stuck here!


title: A repeated-measures LMM... Am I understanding this correctly?

text: So I've been reading West et al's book on Linear Mixed Models as I'm trying to figure out how I should analyze my data, and I was wondering if I've correctly understood what I've ready so far. So, I have repeated measures on a response variable: for each subject, a measure under the control condition and under a "treatment." The subjects also differ in sex, which we predict to be an important factor in the response to treatment. Importantly, some subjects (not all) are measured multiple times, because they went through "treatment" multiple times. So, in those cases, I have a new control measure and treatment measure for each individual. What I think I is a LMM with 3 fixed effects: my response variable, sex, and treatment condition. In addition, I think I will need 1 random effect: subject identification (to correct for the non-independence of some of my data points because some subjects were included twice). First of all, am I missing anything? The book lost me a bit when it started including all of these interactions in their models... I'm also working with very small sample sizes (because I work with endangered and long-lived animals), so I worry about trying to cram too many variables into the model. Second, do my data need to follow any particular distribution? Or do the residuals need to follow some particular distribution? Thanks for taking the time to read...


title: calculating ANOVA 'by hand'

text: I've calculated an ANOVA in R using the `aov()` function, then calculated it semi-manually using my own code. I'm interested in why my mean sum of squares between groups (0.634443) is not the same as that outputted by `aov()` (0.4145), and yet my mean sum of squares within groups is the same as that outputted by `aov()`?               # create fake data     set.seed(50)     x = data.frame(length=rnorm(9),                     site = c(rep('a', 3),                              rep('b', 3),                              rep('c', 3)), stringsAsFactors=F)     # ANOVA summary     summary(aov(length ~ factor(site), x))     # Df Sum Sq Mean Sq F value Pr(>F)     # factor(site)  2  0.829  0.4145   0.514  0.622     # Residuals     6  4.840  0.8067           # calculate sum of squares within, between and total     sumlength <- sum(x$length)            n = 3         overallMean <- sumlength/n         aMean <- mean(x$length[x$site=="a"])         bMean <- mean(x$length[x$site=="b"])         cMean <- mean(x$length[x$site=="c"])           ssw <- sum((x$length[x$site == 'a'] - aMean) * (x$length[x$site == 'a'] - aMean)) +                  sum((x$length[x$site == 'b'] - bMean) * (x$length[x$site == 'b'] -                  bMean)) + sum((x$length[x$site == 'c'] - cMean) * (x$length[x$site ==                  'c'] - cMean))       sst <- sum((x$length - overallMean) * (x$length - overallMean))       ssb <- ((aMean - overallMean) * (aMean - overallMean) + (bMean - overallMean) *              (bMean - overallMean) + (cMean - overallMean) * (cMean - overallMean)) * 3          meanSsb <- ssb/2     # 0.634443     meanSsw <- ssw/6     # 0.8066979


title: Is there a more precise definition of generative models?

text: Intuitively, a generative model is one that we can generate good data from. What I'm looking for is a formal definition of "good data." For example, for any classification model I could generate data by simply sampling really close to the decision boundary. But this is not "good data" because the data will probably be from a very different distribution as the original data. Is there a more precise definition of what constitutes a generative model?


title: AUC for more than two groups?

text: Standard ROC curves look at how setting various thresholds on a continuous measure can be used to predict a two-level ordinal outcome (example: antibody level -> (not sick, sick) ). This can then be integrated over the whole range of thresholds yielding the classical Area Under the Curve (AUC). **The AUC has an interesting alternative definition** : it is also the probability that, for a given pair of subjects, the continuous measure will rank them correctly with respect to the ordinal outcome. In the case where this ordinal outcome has more than two levels (example: not sick, a little sick, very sick) , this alternative definition can be used to get a **natural generalization of the concept of AUC.** My question is: **does this generalization have a name**? Is it a standard or at least known measure? **Are there any more widely used alternatives?**


title: R caret package question

text: library(car)     library(caret)     trainIndex <- createDataPartition(Prestige$income, p=.7, list=F)         prestige.train <- Prestige[trainIndex, ]         prestige.test <- Prestige[-trainIndex, ]         my.grid <- expand.grid(.decay = c(0.5, 0.1), .size = c(5, 6, 7))         prestige.fit <- train(income ~ prestige + education, data = prestige.train,             method = "nnet", maxit = 1000, tuneGrid = my.grid, trace = F, linout = 1)         prestige.predict <- predict(prestige.fit, newdata = prestige.test)         prestige.rmse <- sqrt(mean((prestige.predict - prestige.test$income)^2))      The above was discussed here How to train and validate a neural network model in R?   1. Does caret package run many times for combinations of decay and size? If so, what is the default # of iterations?   2. What is the final choice of decay and size? When I do `summary(prestige.fit)`, the decay is 0.5 and size is 5. Is that the final combination that caret chose as best option?


title: finding $E[X\;|\; Y]$ that $Y=|X|$ in $N(0,1)$

text: suppose $X\sim N(0,1)$. how can calculate $E[X\;|\; Y]$ that $Y=|X|$ please guide me


title: Closed form solution for posterior with scale mixture normal prior

text: Is there a closed form for the posterior mean and variance of an observation drawn with known measurement error if the prior is a scale-mixture normal? To be more clear, here's the implicit model: Y~normal(Z, sigma_Y) (sigma_Y is known, Y is observed) Z~mixture normal where both components have a known shared mean u and each component has known variance sigma_1^2 and sigma_2^2 (WLOG, sigma_2>sigma_1). Overall probability of being in component 2 is a known number p. I suspect that there's a closed form that would allow me to avoid using EM or MCMC technqiues. Has anyone seen one anywhere?


title: Beyond multiple linear regression for longitudinal data?

text: I have a relatively large dataset from a longitudinal study (~100k subjects and ~10 random effects) where the outcome is a real-valued parameter. A simple run of R's `lm` or one of the robust multiple linear regression models like `rlm` says that X1-X7 are all strongly significant, but the overall fit of the line is very poor at R^2=0.09. From summary plots of the data, it seems like there may be a non-linear relationship between the random effects. My goal here is to understand how important each random effect really is in determining the outcome, thus we want to increase the fit of the model while maintaining interpretability. What should be my next steps? Polynomial regression? Gaussian Processes? More manual data analysis to engineer better features? Something else?


title: Elementary approach to higher order asymptotics

text: I am trying to understand “higher order asymptotics”. I find several texts on Likelihood asymptotics, nothing’s easy to read... if you have any nice pointers on this direction, I’ll be interested; however my main question follows. The following ”roadmap“ to higher order asymptotics seems natural to me, and seems something to do and understand _before_ looking at likelihood theory, but I don’t find anything on these lines: Consider $T = f( X )$ where $X$ is a random vector of known distribution; let's say, for simplicity, that $X$ is normal, or multinomial (or the concatenation of several independent multinomials). The Delta method tells me that asymptotically, $T$ is normal, and how to compute its mean and variance (using linear approximation of $f$). I think one can find a better approximation of the distribution of $T$ by   * computing the first moments of $T$ (using linear, or quadratic, or higher order approximations of $f$)   * ”finding“ a distribution with the same first moments **plus** some other “parsimony” criteria to ensure uniqueness (”finding” is not well-defined: at least, being able to evaluate it numerically) Is that possible? Do you know any textbooks/lecture/article going in that direction? **Edit** fg nu gave me some pointers for the second step, that lead me to Edgeworth series. A few references:   * Rothenberg 1984 (given by fg nu)   * Blinnikov 1998 (from the wikipedia article)   * Chapter 16 of Feller’s Introduction to Probability  The first step is rather elementary, however any good pointer is still appreciated.


title: Couple of linear regressions with heteroskedasticity and correlated errors

text: I have two linear models: $$ y_1 = \mu_1 + \beta_1x + \epsilon_1 $$ $$ y_2 = \mu_2 + \beta_2x + \epsilon_2 $$ where $\epsilon_1, \epsilon_2 \sim \text{MVN}(0, \Sigma(x))$, where $\Sigma(x)$ in NOT diagonal. At the moment I'm regressing $y_1$ and $y_2$ on $x$, then I'm calculating the squared residuals $\epsilon_1^2$, $\epsilon_2^2$ and $\epsilon_1\epsilon_2$ and I'm regressing them linearly on $x$: $$ \epsilon_1^2 = \mu_3 + \beta_3x + z_1 $$ $$ \epsilon_2^2 = \mu_4 + \beta_4x + z_2 $$ $$ \epsilon_1 \epsilon_2 = \mu_5 + \beta_5x + z_3. $$ What I need are estimates of all the $\mu$s and $\beta$s. I think that I'm losing efficiency in both regressions (because of the heteroskedasticity in the first and the chi-squared distributed error in the second). Is there any computationally cheap way to improve things?


title: R square in mixed model with random effects

text: For mixed models with random intercepts only, the statistic for R square is $$R^2 = (V_{int\ only} - V_{full\ model}) / V_{int\ only}$$ My question is: How to estimate R square in mixed models with random effects? And estimate R square for each level? Let's say it's a 3-level model with random effects of time at level 2.


title: Is it preferable to normalize the data between [0-1] before running PCA

text: I'm interested to know if it is recommended to scale a matrix from some integer range (e.g 0-255) to [0.0-1.0] before running a PCA. Does it affect the floating point error?


title: Simultaneous confidence intervals for multinomial proportions: pooling replicates or not?

text: I observed 21 times a response that could be one out of 8 categories (multinomial response). I computed the proportions for each category. I did it independently 12 times. Hence I have 12 replicates for each proportion. Is bootstrapping the mean proportion over the 12 replicates for getting the confidence interval a good idea? I tried (see codes below) and got some very large confidence intervals. Actually, according to the bootstrapped CI, there are no significant differences in proportions, while a chi-square test on pooled proportions (over replicates) confirms that there are differences. So, do I have to pool my proportions or not? If I have to, do you know an R function for computing Simultaneous confidence intervals for multinomial proportions? If I am able to keep replicates, is bootstrapping as done below OK?               # generation of a data set: multinomial réponse variable repeatedly measured     sub<-c("s1","s2","s3","s4")#subject     per<-c("p1","P2","p3")#period, RM within subject     tim<-paste("t",1:21,sep="")#tim, RM within period     mydata<-expand.grid(tim=tim,per=per,sub=sub)     #for simplicity, let's consider that period is not a repeated measure     #within subject. Hence, we have 12 trials ("independent") of 21 observations          #random generation of the response variable     require(plyr)     posl<-paste("z",1:8,sep="")     mydata$pos<-factor(NA,levels=posl)         n=nrow(mydata)         mydata$pos<-replicate(n,sample(posl,1,prob=c(0.05,0.05,0.05,0.4,0.3,0.05,0.05,0.05)))          #strong preference for z3 and z4          ##pre-treatment of the above row data, number of time (out of 21) each position was observed     mydata2=ddply(mydata,.(sub,per,pos),summarize,nobs=length(pos),.drop=F)     mydata2[,1:3]=catcolwise(function(x)as.factor(x))(mydata2)     summary(mydata2)          # boxplot of frequencies of occpupancy     require(ggplot2)     mydata2$fobs=mydata2$nobs/21     ggplot(mydata2)+geom_boxplot(aes(pos,fobs))          ###chi2 test     nobsT=ddply(mydata,.(pos),summarize,T=length(pos))     nobsT=nobsT[,2]     chisq.test(nobsT)          #bootstrapping of the mean proportions over the 12 replicates     require(boot)     maf=function(data,index){o=ddply(data[index,],.(pos),summarize,mmea=mean(fobs),.drop=F);as.vector(o[,2])}     bootobj=boot(mydata2,maf,R=9999,stype="i")          confint=ddply(mydata2,.(pos),summarize,mean.prop=mean(fobs))     confint$boot.LOW=-999         confint$boot.UP=-999     for (posi in 1:8){       try({confint$boot.LOW[posi]<-boot.ci(bootobj,0.95,index=posi,type="bca")$bca[4]},silent=T)       try({confint$boot.UP[posi]<-boot.ci(bootobj,0.95,index=posi,type="bca")$bca[5]},silent=T)     }          ####################### plotting boostraped CI     ggplot(mydata2)+geom_boxplot(aes(pos,fobs))+       geom_errorbar(data=confint,aes(pos,ymin=mean.prop-boot.LOW,ymax=mean.prop+boot.UP),lty=2,width=0.5)+       geom_point(data=confint,aes(pos,mean.prop),shape=24,size=2)      References of the CI estimations I would like to run: Glaz, J. and C.P. Sison. 1999. Simultaneous confidence intervals for multinomial proportions. Journal of Statistical Planning and Inference 82:251-262. Sison, C.P and J. Glaz. 1995. Simultaneous confidence intervals and sample size determination for multinomial proportions. Journal of the American Statistical Association, 90:366-369. Paper available at http://tx.liberal.ntu.edu.tw/~purplewoo/Literature/!Methodology/!Distribution_SampleSize/SimultConfidIntervJASA.pdf


title: Correlated predictors in linear model: where does shared variance go?

text: I analyze a dataset consisting of responses to two conditions with a multiple regression analysis. I dummy code the two conditions to retrieve the contribution of each of them to the data. Up to here, no problem, I know what I am doing. Second, I want to check that the difference between the two conditions (Co1 and Co2) is not confounded by another variable of no interest (NoI). I know that in a multiple regression model, the estimates of each condition correspond to the expected change in the data when a change of 1 is observed in its predictor, while all others are held constant. I thus add a covariate to the analysis with NoI and look at the estimates. _Provided NoI does share some variance with Co1 and Co2_ , should I consider that NoI has taken all of the variability of the data that it can account for? If yes, then, Imagine I'd be looking at it the other way around: I would want to know the impact of NoI on my data, but would have also two confounding conditions Co1 and Co2. I'd run the same model and observe that the estimate for NoI is smaller when Co1 and Co2 are accounted for. I would say with the same model that Co1 & Co2 have taken all of the variability of the data that they can account for... _If NoI shares some variance with Co1 & Co2_, this is contradictory with the previous paragraph. So I guess that each NoI and Co takes its share (half?) of the variability it can account for in the data? but how is that share determined? I've read this question and the answer and thought I understood everything: that it depends on the way the sum of squares is split, until I read this comment to find out that we're talking of sum of squares only when _testing_ for effects, not when _estimating_. This confuses me a lot. I don't understand why the problem of sharing variance or sum of squares is different in ANOVA and when doing multiple regression. I get all the more confused when I read elsewhere that ANOVA and multiple regression are one and the same thing. I also came across the idea that we could orthogonalize regressors sequentially, as described here, and done with this code so that they are completely independent and only account for their own share of variance. This is appealing but does not seem like a wide spread practice. In sum, I need help for   1. Understanding how variability in the data is shared/split/not between two correlated regressors.   2. What aspects of ANOVA and multiple regression are one and the same, and what is different?   3. Knowing if it is best practice to sequentially orthogonalize regressors as mentioned above. I hope this question makes sense. I'm really struggling to understand what I'm doing here...


title: Excel CHISQ.TEST

text: In Excel, if I get a value that is less than alpha (i.e. I get a value near zero) for the CHISQ.TEST function, does this mean that the my observed data follows the expected distribution or is it the other way around?


What makes the results differ for fixed-effects models vis-à-vis random effects models? The Cochrane Collaboration's website indicated that two models can produce different results for a meta analysis. I would like to understand this point better.


title: Variance of $X_i / \sum\limits_{j=1}^n X_j$

text: Suppose that $X_1 ,X_2,...X_N$ are independent identically distributed positive random variables. How to find the variance of $\frac{X_i}{\sum\limits_{j=1}^n X_j}$, $i \in 1, 2, \ldots, n$?


title: Difference of unequal sample size

text: I am a newbie in statistics. I have two columns of numerical data of unequal length, motion data. I need difference of the two columns, by shortening the longer column.               766.5 543.8     764.9 544.3     764.0 544.0     765.5 543.7     764.0 543.2     764.4     765.3     765.8      Now I want a third column of data which is the subtraction of column two from column one. So I want to re-size column one data to the length of column two, either by deleting some values or something like that. i.e there are 8 values in first column and 5 values in second column. So I want to resize column one to length 5, with minimal loss in the SD and other parameters. I want to use the results of the difference of two columns as a input to neural network and try to build a predictive model with this data. Are there any methods / algorithms which can help me achieve this, while keeping the accuracy of the shortened column good enough?


title: Unimodal distribution squared input variable

text: When conducting a GLM, why do I have to square my independent value in order to model an unimodal distribution? for example (Matlab code)               %GLM     [logit_rain,dev,stats] = glmfit([(rain) (rain).^2     ],[dependend_variable],'binomial','logit');


title: Dealing with 0 values when calculating the mle for a Dirichlet distribution

text: I have $N$ pmfs, and for each each $L$ samples. Each sample has a variable amount of $x$ values, but the $x$ values that they have can be matched. So for example: $$sample_1 \rightarrow\ x_1 = 0, x_2 = 0, x_3 = 0.2, x_4 = 0.4, x_5 = 0.4$$ $$sample_2 \rightarrow\ x_1 = 0.3,x_2=0, x_3 = 0.4, x_4 = 0.3,x_5=0$$ I'm using a python program to calculate the mle from the samples which is a port of Thomas P. Minka's Matlab Fastfit code (Estimating a Dirichlet Distribution). The problem is that, for fitting, it sums over $logp+psi(\sum^k a_k) - logp$. Since some of the $x$ values are 0, some logp values are -inf. Therefore, summing over this makes everything -inf. How can I deal with 0 values when calculating the mle for a Dirichlet distribution?


title: Principal component analysis with known uncertainty varying across both samples and variables

text: I want to do dimensionality reduction on a data set $X_{ij}$. In this case, $i$ indexes samples and $j$ indexes a large number of variables (densities at different locations in space). The units of all elements are the same, but the uncertainties (which I can estimate) vary quite a lot, mostly across samples. Each $X_{ij}$ is Gaussian, with standard deviation $\sigma_{ij}$ varying across both samples and variables. My problem is that doing standard PCA without taking into account the uncertainties will tend to produce results which are dominated by variance due to measurement error. What I would like to do is estimate some similar quantities to the principal components, but which will be more robust to outliers with large $\sigma_{ij}$. The variables $X_{ij}$ are drawn from $\mathcal{N}(\mu_{ij},\sigma_{ij})$. I have one such sample for each $i$ and $j$, and I would like to calculate a transformation involving $X_{ij}$ and $\sigma_{ij}$ which will estimate principle components of the matrix $\mu_{ij}$, for example. It is difficult for me to make progress analytically because I cannot calculate how uncertainty in $X$ propagates to its principle components. I can generate surrogate data to confirm that uncertainty dominates the principle component analysis, and I can also do this to validate any alternative techniques. My main goal is to look for independent modes across the variables, so my concerns are with dimensionality reduction for presentation, description and simplification, rather than physical interpretation.


title: Drawing smooth 2D curves with 2D spline in R

text: The problem I am running into here is the following: there is a set of 4 points in 2D, and I need to plot a smooth curve in 2D space that runs through these points. Is there a package in R to do this?


title: Trouble Interpreting Fixed Effects and What Features Are Fixed

text: I have a large panel dataset. It follows a number of schools over 10 years. It is essentially set up as follows: Year School VariableOfInterest ExplanatoryVariable InstrumentForExplanatoryVariable I want to see how well VariableOfInterest is described as a linear function of ExplanatoryVariable. I want to compare the results of a simple regression and instrument variable regression with a few fixed effects models. I'd like to see how well VariableOfInterest is described when fixing school effects and then fixing both school and year effects. I'm unsure of how to look at only school fixed effects. I believe I've already correctly done a correct regression fixing both school and year. In Stata the command is just: xtreg VariableOfInterest ExplanatoryVariable, fe. Does anyone have any insight? I apologize if this is a basic or misguided question.


What are some examples of applied machine learning problems that requires using mixed models? I'm just introduced to the notion of mixed models. As I understand it, it is a combination of parametric model such as logistic function and non-parametric model attempting to model random effects that can be intractable with parametric models. What are some real-world examples that might be useful to develop some intuition or big-picture understanding of this method?


title: Power analysis in meta-analysis: Compute the maximum significance needed by an additional sample for meta-analysis to reach genome-wide significance

text: I have an association study analysis (GWAS) that does not reach genome-wide significance (max _p_ -value < 10E-5), and there is a potential additional sample (so with known N and known effect size/SE for each association) that we want to use together with the GWAS to conduct a fixed effect meta-analysis. However, before including this sample, I am looking to run a power analysis (hopefully with the help of software) that can help me understand the maximum significance level that this additional sample would have to reach for the meta-analysis to reach genome-wide significance. Can anyone please help understand if this is possible, and if you know any R packages or software that can help with this power analysis for fixed-effect meta-analysis? Thanks.


title: Transforming scores between different rating scales – Importance of response options

text: I'm working on a meta analysis of the mean of a single variable (using arithmetic mean as the effect size, rather than the relationship between two variables). The attitude I'm assessing has been measured by dozens of different measures, which have many different rating scales. I'm using the following equation to transform scores between two different measures to put them on the same metric: $$X_2 = \frac{(X_1 - \min_1)(\max_2 - \min_2)}{\max_1 - \min_1}+\min_2$$ Where $X_2$ = the score on the second measure; $X_1$ = score on the first measure I want to transform; $\min_1$ and $\max_1$ = lowest and highest possible scores on measure 1; and $\min_2$ and $\max_2$ = lowest and highest possible scores on measure 2. (Card, N., _Applied Meta-Analysis for Social Science Research_ , 2011, pg. 148). Assume that any two measures I'm trying to compare are assessing the same construct. I know I can use this if the two measures have the same rating scale (1–5, strongly agree to strongly disagree). My questions are:   1. What about two measures that have the same number of response options, but slightly different wording? For example, both are 5 point scales, but measure 1 is:  > 1 = strongly disagree, 2 = disagree, 3 = neutral, 4 = agree, 5 = strongly > agree while measure 2 is: > 1 = completely disagree, 2 = disagree up to a point, 3 = neutral, 4 = agree > up to a point,   >  5 = agree completely How critical is it that they have exactly the same wording? (I'm inclined to think these are roughly the same and could be combined).   2. What about measures with different numbers of responses? Ex: measure 1 is the same as the example above, but measure 2 has a 1–7 rating scale with only endpoint anchors of strongly disagree/strongly agree, measure 3 has a 1–6 rating scale of:  > 1 = strongly disagree, 2 = moderately disagree, 3 = slightly disagree, 4 = > slightly agree, 5 = moderately agree, 6 = strongly agree and measure 4 has a 1–7 rating scale of: > 1 = strongly agree, 2 = agree, 3 = slightly agree, 4 = neutral, 5 = slightly > disagree, 6 = disagree, 7 = strongly disagree Converting between 5 and 7 point scales seems less problematic than converting between 5 and 6, since you lose the neutral midpoint. However, I welcome any advice on this as it will inform my methods moving forward. Thank you!


title: Moments and exponential distribution

text: Stuck on how to solve this, can't seem to find the answers online at all a. The random variable $X$ has an exponential distribution with parameter and probability density function $f_X(x) = \theta e^{-\theta x}$ , where $x>= 0$. Obtain the moment generating function $m_X(t) = E[e^{tX}]$ of $X$. For what values of t is this defined? Suppose that $X_1, . . . ,X_n$ are mutually independent and identically distributed random variables each having an exponential distribution with parameter $\theta$. Let $S_n = X_1 + · · · + X_n$. Write down, together with a brief statement of what results you have used, the moment generating function of $S_{n}$ . b. Suppose that $y=\frac{\Theta S_{n}}{\sqrt{n}}-\sqrt{n}$ has moment generating function mY (t). Show that $m_{y}(t)=e^{-t\sqrt{n}}\left ( 1-\frac{t}{\sqrt{n}} \right )^{n}$ Obtain $logmY (t)$ and deduce what happens as n -> $\infty $. What do you conclude about the distribution of $Y$ for large $n$? By expanding $mX(t)$ as a power series in $t$, or otherwise, deduce that $E[Xr] = r!/$$\theta$$^{r}$ for r = 1, 2, 3, . . ..


title: What particular measure to use? Multiple regression or MANOVA?

text: I've been doing a research paper - Effects of Workplace Bullying on Employees' Productivity, Self-Confidence and Self-Esteem. So I have 1 IV (Workplace bullying) and 3 DVs. what can I use to interpret it? Another thing: if it's MANOVA, can I execute it without having to use factors?


title: Calculating Prediction Interval

text: I have the following data located here. I am attempting to calculate the 95% confidence interval on the mean purity when the hydrocarbon percentage is 1.0. In R, I enter the following.               > predict(purity.lm, newdata=list(hydro=1.0), interval="confidence", level=.95)        fit      lwr      upr     1 89.66431 87.51017 91.81845      However, how can I derive this result myself? I attempted to use the following equation. $$s_{new}=\sqrt{s^2\left(1+\frac{1}{N}+\frac{(x_{new}-\bar x)^2}{\sum(x_i-\bar x)^2}\right)}$$ And I enter the following in R.               > SSE_line = sum((purity - (77.863 + 11.801*hydro))^2)     > MSE = SSE_line/18     > t.quantiles <- qt(c(.025, .975), 18)     > prediction = B0 + B1*1     > SE_predict = sqrt(MSE)*sqrt(1+1/20+(mean(hydro)-1)^2/sum((hydro - mean(hydro))^2))     > prediction + SE_predict*t.quantiles     [1] 81.80716 97.52146      My results are different from R's predict function. What am I misunderstanding about prediction intervals?


title: How to compute a limit of $L_p$ norms?

text: If $\mu$ is a probability measure then $$\left\| f \right\|_{\infty} = \lim_{p \rightarrow \infty} \left\| f \right\|_{p}$$ where $\left\| f \right\|_{p} = \left( \int f^p d\mu\right) ^{1/p}$ and $\left\| f \right\|_{\infty}$ is the essential supremum of $f$ with respect to $\mu$. Would someone be so kind as to suggest ways to prove this?


title: Visualizing data spread along both y- and x-axis in 2D scatterplot

text: I'm currently in the process of doing a literature review of a brain atrophy measurement, called BPF, commonly used in research regarding neurodegenerative disorders. I have searched for publications having investigated BPF and for each study population I have come across I have noted population size, mean BPF, mean age and the standard deviations of these two means. I do not have access to source data. I would like to visualize all these populations in a 2D bubble-plot (scatterplot with variable dot-sizes), with BPF on one axis and age on the other, each data point representing one population. I would like the size of the points (bubbles) to be representing population size. I would like to include error bars for each point representing +-SD of age (bar running along the age axis) and +-SD of BPF (bar running along the BPF-axis). I am currently using SPSS21 and have myself come so far as to be able to create the bubble-plot (using size(PopulationSize) for the dots) with error bar for one of the variables (using SPSS high-low-close graph), but I can not seem to get a second error bar along the other axis in there. Any suggestions on how to accomplish making this graph, using SPSS or other software, would be greatly appreciated. Thank you for your help!


title: Using orthoginal comparison for non parametric data

text: Can I use an orthoginal comparison with non parametric data? What about using a post hoc test with my non parametric data? I am counting the number of objects within a nucleus and based on my situation the data should follow a Poisson distribution. This assumption is validated with a normality test.


title: Comparing effectiveness of various treatments over control via matched pairs

text: Imagine this experiment: assume that I have two sets of twins. For the first pair of twins, I keep one as a control, and the second twin I give drug A. For the second pair of twins, I keep one as a control and I give the second twin drug B. Now, I want to show that the effectiveness of drug A over control is bigger than the effect of drug B over control. In each case, I measure a response in each twin like say body temperature. My goal would be to prove that drug A lowered body temperature more than drug B compared to the control, a single p-value summarizing the experiment for instance. I replicate the whole experiment three times. So, for drug A, I have three pairs of numbers, each making a comparison between a pair of twins, one given the drug and one left as a control; and the same for drug B. How do I go about analyzing an experiment like this? (For a lot of technical reasons related to my application, I have to work with matched pairs; and I'm not really experimenting on humans.)


title: How to stabilize features over time?

text: I have in my model some features, which seem to have some super power to explain the training data. To be more specific, the model is to predict clicks and the features are the click-through-rate(CTR) of some granularity. The CTRs of certain granularity are the same in both training and testing set. Althought we gained a lot of precision in training set by adding such features, the performance gain in testing set was limited. Because I seperated training/testing dataset over time, so I suspected that it was caused by feature being unstable over time. (and I also verified this hypothesis by calculating granularity CTR over different period). So my question is how to stablize such features, and if it is not possible, how shall we draw a conclusion whether such feature is only noise?


title: Quadratic weighted kappa strength of agreement

text: In the case of the kappa-value there are some attempts to qualify how good or bad the agreements are. For example Landis & Koch in the article _The Measurement of Observer Agreement for Categorical Data_ talks about "strength of agreement" based on kappa values:                Kappa       Strength of agreement     =====       =======================     0.0-0.20     Slight     0.21-0.40    Fair     0.41-0.60    Moderate     0.61-0.80    Substantial     0.81-0.90    Almost perfect      My question is if there are some attempts to qualify the strength of agreement" based on "quadratic weighted kappa values" in the same way. Any references about attempts to define "generic strength of agreement" using quadr. weighted kappa? My assumption is that it would not be meaningful to define "strenght of agreement" based on kappa and to use the same qualification using quadratic weighted kappa. They are different values.


title: How can I measure the effects certain events have on the frequency of other events over time?

text: **EDIT: added more details following @kjetil comments** I have the following problem:   * I monitor one stream of events of type A - those events can be considered instantaneous.   * I also monitor additional streams of events of types B1, B2,... Bn - those have a time interval associated with them.   * I want to determine if the frequency of type A events _increases_ significantly during any specific B types (in other words, I want to determine exactly which B types have a significant _positive_ influence during the intervals of their events)   * Many instances of the same type of event can occur simultaneously   * The effect is likely to be linear in the number of concurrent events of the specific B type   * There are probably many B types that don't effect A-events   * B types of events are not necessarily the only thing that influence type A events   * All events probably have some complex distribution that changes depending on the day/hour/week/month/previous events - not a simple Poisson   * I _constantly_ monitor all streams - I want to discover when a B-type starts to influence the frequency of A-event and I would also want to know if at some point in the future it stopped influence it (and this should continue indefinitely)) For example:   * A-type events are power outages    * B1-type events correspond to opening and closing of air conditioners.   * B2-type events correspond to opening and closing of televisions. Do you have any idea how I should approach this problem? (general references are fine) Thanks!


title: Estimated critical value for hypothesis testing

text: For the classical simple linear regression model I have derived an hypothesis test for $H_0\colon \left\\{\frac{y^*(x^*)-x^*}{\sigma}>1 \right\\}$ where $x^*$ is a given value of the covariate $x$ and $y^*(x^*)=a + b x^*$ is the theoretical mean of the response $y$. I use the test statistic $t=\frac{\hat{y}(x^*)-x^*}{\hat\sigma}$. Something very nice happens: at the boundary $\left\\{\frac{y^*(x^*)-x^*}{\sigma}=1 \right\\}$ of $H_0$, there is only one possible distribution of the random variable $t$. Hence for a given $\alpha \in ]0,1[$ I can find the critical value $C$ such that the type 1 error of the hypothesis test is exactly $\alpha$ when the rejection rule is $t<C$. Very nice, but now I want to do the similar hypothesis test for Deming regression, and things are not so nice : the asymptotic distribution of $t$ has several possible distribution at the boundary of $H_0$. Hence I have derived an estimated critical value $\hat C$ and I use the rejection rule $t<\hat{C}$. Simulations show that the type 1 error is approximately well controlled: for a desired type 1 error $\alpha$ the effective type 1 error is close to $\alpha$. But I wonder whether there are some pitfalls with my procedure ? Do you know other examples where one similarly uses an estimated critical value ? And are there some known pitfalls with these examples ? I think this is not rigorously correct to use a test statistic whose distribution is not uniquely determined at the boundary of $H_0$, but I don't know how to do otherwise.


title: Using priors in classification for regression

text: I am ultimately trying to perform a regression task - for this example, let's say I'm trying to determine the height (in pixels) of a person in an image. However, rather than doing regression, I am solving this problem by hypothesising a few different possible heights (e.g. by detecting parts of the image that look like 'head' and 'feet') and then finding which hypothesis is the most likely. For want of a better term, I am calling this 'classification for regression'. I am using training data to learn a classifier, and my question is how the _likelihood_ and _prior_ can be used in this context. In my mind, there are several pieces of information that can be levered:   1. Goodness of fit of the hypotheses. In this case, this might be how much the detected feet look like feet and how much the detected heads look like heads   2. Overall likelihood of matching. e.g. how much does the bit between the detected feet and head look like a body.   3. Relation to prior knowledge about body size. Assuming absolute sizes are known, this could be a normal distribution over the range of peoples' heights. Which of these pieces of information is a prior, and which can be used to infer the likelihood? (1) and (2) seem to be likelihood data, while (3) seems to be a prior, but why not use (3) as part of the likelihood? Are there any more general points and tips for this sort of task?


title: How does Logistic regression classifier modelize the dataset?

text: I'm working on a system that be able to detect the hand contour. So I have 270 instance in my dataset: 7 class of hand contour, 8 feature vectors of each instance.   Firstly, I used Weka to determine which is the best classifier for my dataset. And I found out Logistic regression classifier.   But, the problem is I don't know 'How does the Logistic regression classifier work?' or 'WHAT IS THE MAIN PONIT (IDEA) OF LOGISTIC REGRESSION CLASSIFIER?'. I've been searching on google for a long time but I still can't find the answer.   Here's what I've found that should be useful:   1 1992-JSTOR-Cessie-Logistic_Regression   [2] [Chapter 4 Logistic Regression as a Classiﬁer]1 Thanks in advance!


title: Negative value in shape parameter of gev.fit

text: I got negative value for shape parameter by using `gev.fit` in `ismev` package of R. What does it indicate? Does it mean that it belongs to another distribution instead of generalized extreme value (GEV) distribution? here is the output of `gev.fit`:               $conv [1] 0          $nllh [1] 12.96933          $mle [1] 13.796261 1.847000 -1.084086          $se [1] 0.001307101 NaN NaN      The value of shape parameter=-1.084086, it is too negative which doesn't approximate zero.


title: How to obtain detrended data and demeaned data in Stata?

text: I would like to run KSS test (KAPETANIOS, G. – SHIN, Y. – SNELL, A. 2003: Testing for a UnitRoot in the Nonlinear STAR Framework) , but I don't know how to obtain demeaned data and detrended data. Could you describe me the process in STATA?


title: Inclusion of lagged dependent variable in regression

text: I'm very confused about if it's legitimate to include a lagged dependent variable into a regression model. Basically I think if this model focuses on the relationship between the change in Y and other independent variables, then adding a lagged dependent variable in the right hand side can guarantee that the coefficient before other IVs are independent of the previous value of Y. Some say that the inclusion of LDV will biase downward the coefficient of other IVs. Some others say that one can include LDV which can reduce the serial correlation. I know this question is pretty general in terms of which kind of regression. But my statistical knowledge is limited and I really have a hard time figuring out if I should include a lagged dependent variable into a regression model when the focus is the change of Y over time. Are there other approaches to deal with the influence of Xs on the change of Y over time? I tried different change scores as DV as well, but the R squared in that situation is very low. Thanks a lot!


title: Interpretation of binomial parameters in zeroinfl in R

text: In the help page of `zeroinfl`, it says " **a binary model is used that captures the probability of zero inflation**." But it seems it is modeling the opposite of that (see the code below). Am I interpreting the verbal description wrong? Because arguably, if p is modelled, 1-p is automatically modelled, and vice versa.               n <- 10000     y <- rpois(n,3)     p <- 0.3        # probability of zero inflation     zi <- rbinom(n,size=1,prob=p)==0     y <- y*zi       # zero inflated poisson          library(pscl)     (model <- zeroinfl(y~1|1,dist="poisson"))     exp(coef(model)[2])/(1-exp(coef(model)[2])) # this seems 1-p instead of p


title: Is it wrong to jitter before performing Wilcoxon test?

text: I wrote a script tests the data using the `wilcox.test`, but when I got the results, all the p-values where equal to 1. I read in some websites that you could use jitter before testing the data (to avoid ties as they said), I did this and now I have an acceptable result. Is it wrong to do this?               test<- function(column,datacol){         library(ggplot2)     t=read.table("data.txt", stringsAsFactors=FALSE)     uni=unique(c(t$V9))         for (xp in uni) {                for(yp in uni) {                 testx <- subset(t, V9==xp)                 testy <- subset(t, V9==yp)                 zz <- wilcox.test(testx[[datacol]],jitter(testy[[datacol]]))                 p.value <- zz$p.value               }             }     }      * * * This is the output of `dput(head(t))`               structure(list(V1 = c(0.268912,     0.314681, 0.347078, 0.286945,      0.39562, 0.282182), V2 = c(0.158921, 0.210526, 0.262024, 0.322006,      0.133417, 0.283025), V3 = c(0.214082, 0.166895, 0.132547, 0.147361,      0.09174, 0.169093), V4 = c(0.358085, 0.307898, 0.258352, 0.243688,      0.379224, 0.2657), V5= c(-0.142223, 0.010895, 0.14655,      0.08152, 0.02116, 0.030083), V6 = c(0.096408, -0.091896,          -0.331229, -0.446603, -0.088493, -0.262037), V7` = c(1.680946,      1.649559, 1.534401, 1.130529, 3.441356, 1.211815), V8 = c("NC_000834",  "NC_000844",     "NC_000845", "NC_000846", "NC_000857",     "NC_000860" ), V9 = c("Chordata",     "Arthropoda", "Chordata", "Chordata",      "Arthropoda", "Chordata"), V10 =     c("???:???", "Diplostraca",      "???:???", "Rheiformes", "Diptera",     "Salmoniformes"), V11 = c("???:???",     "Branchiopoda", "Mammalia", "Aves",     "Insecta", "Actinopterygii" )), .Names     = c("V1", "V2", "V3", "V4", "V5", "V6",  "V7",     "V8", "V9", "V10",     "V11"), row.names = c(NA,  6L),     class = "data.frame")      The data is very large, and that's the thread I started and they told me it might be wrong to do this **Note** This question comes from tex.SE: generating PDFcontain R output inside latex table


title: Variance of a distribution's product with itself

text: I have a distribution $X$. By playing around with random samples from $X$, I've determined that $Var(X^i) > Var(X)$ where $i > 1$. However, I can't seem to find a formula for the expected variance of $X^i$, or why it should be greater than $Var(X)$. Moving away from normal distributions, should it generalize that the scale parameter of any $X^i$ will be greater than $X$?


title: OpenBugs vs. JAGS

text: I am about to try out a BUGS style environment for estimating Bayesian models. Are there any important advantages to consider in choosing between OpenBugs or JAGS? Is one likely to replace the other in the foreseeable future? I will be using the chosen Gibbs Sampler with R. I don't have a specific application yet, but rather I am deciding which to intall and learn.


title: over estimating ARIMA model

text: When an ARIMA model over estimates, i.e if the forecast from the model is always higher than the actual value, what's the cause? The model i used was SARIMA (1,1,0)x(1,0,0). Below is the time series data, the ACF & PACF and the time series plot. Please help me out with identifying the model.The 95% confidence limit is 0.183. I have the actual values for the year 2011 and 2012. The excel file of my data ![enter image description here](http://i.stack.imgur.com/RbrNO.png) ![enter image description here](http://i.stack.imgur.com/JUbXo.png) ![enter image description here](http://i.stack.imgur.com/1ebP1.png)


title: How to decide on the MLE when pmf is 0?

text: Suppose you have $\theta=\\{1,2\\}$ and the sample of (0,1,2) with the task of finding MLE: \begin{array} {|c|c|c|} \hline x & p(x|\theta=1) & p(x|\theta=2) \\\ \hline 0 & 1/2 & 1/4 \\\ \hline 1 & 0 & 1/4 \\\ \hline 2 & 1/3 & 0 \\\ \hline \end{array} I can't compare $\prod_i p(x_i|\theta)$, because of the zeroes. Even though I can compare probability for a single observation, I'm still having trouble quantifying the difference. $p(x=1|\theta=2)>p(x=1|\theta=1)$ for x=1, but can I do a similar comparison for the whole vector? ### Edit The original question is > One observation is taken on a discrete random variable X with pmf f(x|θ), > where θ∈{1,2,3}. Find the MLE of θ. This is followed by a 4x5 table which includes zeroes in every column. I think the idea behind this question is that we have a distribution that doesn't describe population well. Do we just say that the method of MLE is not applicable here?


title: What is the relation between statistics theory and decision theory?

text: I was wondering how statistics and decision theory are related? It looks to me all the statistics problems/tasks can be formulated in decision theory. Also problems in decision theory can be formulated in statistics/probability problems, or in deterministic way. So is statistics just a part of the problems studied in decision theory? Or are they just two theories with overlapping and neither falls completely inside the other? But I have to admit that I don't have a systematic big picture of what topics statistics theory and decision theory are covered respectively, and would like to here some of your point of view. Thanks and regards!


title: Which method(s) for forecasting time series of event durations

text: I have the $N$ individuals each observed for $T$ days. For each individual I have some basic demographic data. Each $n$ individual, during the observed time $T_n$ may experience event $E$ which is recorded as the duration of E on each of the observed day. For example for a given individual $n$ I have "$D_1n,D_2n,D_3n$ which are valid for the whole observed time, plus a vector of duration of $E$ events on each observed day as follows:               Day Event duration(min)     1   0     2   10     3   50     ...      No distribution assumption can be made for $E$ durations, either in patient or in population level. I would like to to investigate if a predictive model of Event duration can be constructed for each individual and for the population overall. For example, **given** an already observed time frame (the event duration history that we have already observed) and the observed demographics (of the patient which are valid for the complete timem past and future), with what confidence can we forecast what each individual will experience the next X days? I am not so familiar with that kind of analyses, I have a biostatistic background and from what I have researched thus far this is more of econometrics problem). I would appreciate any help or any references to start with. Many thanks.


title: Value between bigger and smaller set

text: I must measure a value. To do so, I only have samples with known value and a measuring tool that can tell if the sample is bigger or smaller than the value to measure. Once the measurement is complete, I've got two set: one which contains values that are smaller and one that contains values that are bigger. For example, with those set [-5,-3,-5,-1,-2,-3] and [1,4,3,2,5,3] I can assume that the value has a high probability to be between -1 and 1. But with the two sets [-5,2,-3,1,-2,-6,-4] and [7,3,-1,2] it is less clear. Which tools can I use to find the range or the value and it's degree of confidence? Thanks you. Edit: as my question was not well formulated, here is how the data is gathered: There is a reference sound with an unknow frequency. Then people are given sounds of other frequencies and they must tell if the tone is higher or lower. Thus the two sets. With the data gathered, I need to find the unknown frequency and how much I can trust the result.


title: Is the generalized Pearson’s chi-square statistic a pseudo $R^2$?

text: Is the generalized Pearson's chi-square statistic viewed as a pseudo $R^2$? I think yes, because a pseudo $R^2$ is a generalization of the form of $R^2$, and the generalized Pearson's chi-square statistic has a form with the residuals in $R^2$ in its numerators. But I am not sure. Thanks.


title: Is random forest applied only to continuous response variable?

text: I am trying to apply random forest on a binary response variable, but it's saying the response variable has 5 or fewer unique values. Was it happening because random forest works only with the continuous response variable?


title: Equivalence of a slope coefficient in a regression using individual and aggregate data

text: I want to compare the coefficients of the two linear regressions: $$ \begin{array}{c} y_{ij} &=& \beta_0 + \beta_1 X_{j} + u_{ij} \\\ y_{j} &=& \alpha_0 + \alpha_1 X_{j} + u_{j} \\\ y_{j} &=& \frac{1}{N_j} \sum_{i=1} ^ {N_j} y_{ij} \end{array} $$ so that the second equation differs only by using the mean of group $j$ instead of individual observations $y_{ij}$ as the dependent variable. Is there an easy way to show that $\hat{\beta}_1 = \hat{\alpha}_1$, where hat denotes OLS estimate?


title: Significance of a source phase in a complex ICA

text: I am using complex-valued ICA to extract sources for complex-valued sensor data. One of the three ambiguities for complex ICA is phase ambiguity, i.e., phase rotation $\exp(i\theta_k)$ of the sources $s_k$ (in addition to permutation and scaling ambiguities, as for real ICA). Yet the complex linear ICA model $\bf{z}=A\cdot \bf{s}$ starts from the general form that all three quantities $\bf{z},A,\bf{s}$ are complex, so the merit of this assumption for $\bf{s}$ eludes me. Is there any meaning that can be given to the phases of the $s_k$, or is $|s_k|$ all that can be reported?


title: Bayesian analysis of data

text: I have a big dataset in the form: $X_1, X_2, X_3, X_4, Y$. All the $X_i, i \in {1,...,4}$ come from different unknown distributions and $Y$ follows a bernoulli distribution, so it can take only values in $\\{0, 1\\}$. The problem is how I can predict the value of $Y$ given the observations for $X_i$. From what I know, this can be solved using binary logistic regression. Also, we can set a threshold of 0.95, lets say, for the probability of Y=1 in order to estimate the value of Y (i.e. if P(Y=1)>=0.95, then we estimate that Y is 1) with some given confidence.   The question that arises here is whether we can use (hierarchical?) Bayesian models to approach the same problem of estimating the value of Y. Any suggestions? Also, is there any paper, book chapter, or tutorial that addresses the same problem? From what I have seen, the statistics textbooks cover the trivial case of well-known distributions where we can calculate everything analytically. But what about my case where I have multidimensional experimental data that can follow an arbitrary distribution? Thanks.


title: Getting marginal effects after a panel oprobit regression in Stata (using gllamm package)

text: I am trying to estimate the number of companies entering certain markets using panel data. To do so, I ran an ordered probit regression in `Stata` using the `gllamm` package. The number of companies goes from 0 to 5. The command is as follows (where `id` = a fixed market): `gllamm NCompanies independent_variables, i(id) link(oprobit)` The outcome, if I understood it correctly, gives the cut points of a latent variable and the impact each independent variable has on the latent variable. So far so good to me, I ran into troubles when I tried to estimate the marginal effects using the commands `mfx` or `margins`. When I type `mfx, dydx` I get exactly the same results as displayed by the `gllamm` (?!) minus the cut points. So, I don´t know if the `gllamm` was already showing the marginal effects or if it was showing the impact on the latent variable... if I type `mfx, eydx` the outcome differ from the previous `mfx` \- so, I am tending to believe that it was not a bug on the `mfx` command, but that `gllamm` was already yielding the marginal effects. To check this, I tried to the `margins` command: `margins, atmeans continuous` this gives a completely different answer then the previous ones - and the command `margins, atmeans dydx(independent_variables) continuous` does not converge. Does any one knows the answer to the following questions?   1. Is the `gllamm` already giving the marginal effects in the above case?   2. Am I using the `margins` command in the wrong way if I want to compute the marginal effects of each independent variable (they are all continuous)? If so, what should I do instead?


title: Q&A Sites for machine learning and NLP

text: So my question is asking for recommending good Q&A sites about machine learning, NLP or data mining. I think guys here may often go to these relevant sites.


title: How to correctly use the GPML Matlab code for an actual (non-demo) problem?

text: I have downloaded the most recent GPML Matlab code GPML Matlab code and I have read the documentation and ran the regression demo without any problems. However, I am having difficulty understanding how to apply it to a regression problem that I am faced with. The regression problem is defined as follows: * * * Let $\mathbf{x}_i \in \mathbb{R}^{20}$ be an input vector and $\mathbf{y}_i \in \mathbb{R}^{25}$ be its corresponding target. The set of $M$ inputs are arranged into a matrix $\mathbf{X} = [\mathbf{x}_1, \dots, \mathbf{x}_M]^\top$ and their corresponding targets are stored in a matrix $\mathbf{Y} = [\mathbf{y}_1 - \mathbf{\bar{y}}, \dots, \mathbf{y}_M-\mathbf{\bar{y}}]^\top$, with $\mathbf{\bar{y}}$ being the mean target value in $\mathbf{Y}$. I wish to train a GPR model $\mathcal{G} = \lbrace \mathbf{X}, \mathbf{Y}, \theta \rbrace$ using the squared exponential function: $k(\mathbf{x}_i, \mathbf{x}_j) = \alpha^2 \text{exp} \left( - \frac{1}{2\beta^2}(\mathbf{x}_i - \mathbf{x}_j)^2\right) + \gamma^2\delta_{ij}$, where $\delta_{ij}$ equals $1$ if $i = j$ and $0$ otherwise. The hyperparameters are $\theta = (\alpha, \beta, \gamma)$ with $\gamma$ being the assumed noise level in the training data and $\beta$ is the length-scale. To train the model, I need to minimise the negative log marginal likelihood with respect to the hyperparameters: $-\text{log}\, p(\mathbf{Y} \mid \mathbf{X}, \theta) = \frac{1}{2} \text{tr}(\mathbf{Y}^\top\mathbf{K}^{-1}\mathbf{Y}) + \frac{1}{2}\text{log}\mid\mathbf{K}\mid + \,c,$ where c is a constant and the matrix $\mathbf{K}$ is a function of the hyperparameters (see equation k(xi,xj) = ...). * * * Based on the demo given on the GPML website, my attempt at implementing this using the GPML Matlab code is below.               covfunc = @covSEiso;     likfunc = @likGauss;     sn = 0.1;     hyp.lik = log(sn);     hyp2.cov = [0;0];     hyp2.lik = log(0.1);     hyp2 = minimize(hyp2, @gp, -100, @infExact, [], covfunc, likfunc, X1, Y1(:, n));     exp(hyp2.lik)     nlml2 = gp(hyp2, @infExact, [], covfunc, likfunc, X1, Y1(:, n));     [m s2] = gp(hyp2, @infExact, [], covfunc, likfunc, X1, Y1(:, n), X2);     Y2r(:, n) = m;      X1 contains the training inputs X2 contains the test inputs Y1 contains the training targets Y2r are the estimates from applying the model n is the index used to regress each element in the output vector Given the problem, is this the correct way to train and apply the GPR model? If not, what do I need to change?


title: Standardized regression coefficient ($\beta$) in multi-linear regression by groups?

text: As you see, I have more than 500 groups of data. how can I get the std.coef in mult-linear regression for each group? Until now, I only know how to get one group by using `lm.beta`. For example:               library(QuantPsyc)     p <- read.csv('Data1.csv')     model1 <- lm(YCOORD110 ~ TEMMIN + TEMVAR)     lm.beta(model1)       For groups, several codes like `by`, `apply` can be used, but seem much more complicated beyond my control. very appreciate for your answers! ![enter image description here](http://i.stack.imgur.com/VdwVd.png)


title: Dummy Variable Trap in Stata

text: So basically I am using the `wine.dta` file and have a range of dummy variables for 6 regions. I have excluded 1 of the regions (regressing on 5 of the dummy variables) but for some reason _Stata_ says there is still a collinearity problem and only gives 4 coefficients... Why is there still a collinearity problem when I have already taken out one of the dummy variables to avoid the dummy variable trap.


title: Path analysis or a full SEM?

text: I have been researching this topic for some time now and was after some opinions or comments. Are there any situations where a path analysis (a form of SEM) might be more appropriate than a full SEM? Specifically, if measures used are standardised and often very long making a measurement model extremely complex and subsequently the increase in the parameters to be estimated places a greater demand on sample size, I wondered whether it would be justifiable to use path analysis as a compromise? The focus is on the relationships between the 'latent variables' Which I understand are treated as observed in a path analysis. Measurement models account for the fact that measurement will never be perfect but many tests performed in SPSS, correlational and tests of difference, I assume are on this basis I.e. that there is no error in measurement? I have mixed feelings on this and the literature is equally mixed. Many articles say they are using SEM but do not discuss the specification of measurement models and makes me wonder whether they are using path analysis instead. A recent course I have been on suggests that this is acceptable but again I appreciate there are limitations and wanted to reach out to others that have used these techniques. I hope this makes sense. I suppose I want to make sure that I have considered all avenues and feel confident in my decisions.


title: Average in complexity

text: Lets say that complexity if individual is ranked in ordinal scale as : 1 = simple 2= moderate 3= complex 4 = very complex, how would i indicate the average of complexity by taking from 100 respondents? Thanks


title: use spss to analyse a cluster randomized data

text: I have data from a cluster randomized trial, that is, each subject in the same cluster has received the same treatment. Do chi square test and t-test need some adjustments? More precisely, which analysis have I to carry out in Spss to compare frequency in contingency tables and to compare means for quantitative variables?


title: Investigating robustness of logistic regression against violation of linearity of logit

text: I am conducting a logistic regression with a binary outcome (start and not start). My mix of predictors are all either continuous or dichotomous variables. Using the Box-Tidwell approach, one of my continuous predictors potentially violates the assumption of linearity of the logit. There is no indication from goodness-of-fit statistics that fit is problematic. I have subsequently run the regression model again, substituting the original continuous variable with: firstly, a square root transformation and secondly, a dichotomous version of the variable. On inspection of the output, it seems that goodness-of-fit improves marginally but residuals become problematic. Parameter estimates, standard errors, and $\exp(\beta)$ remain relatively similar. The interpretation of the data does not change in terms of my hypothesis, across the 3 models. Therefore, in terms of usefulness of my results and sense of interpretation of data, it seems appropriate to report the regression model using the original continuous variable. I am wondering this:   1. When is logistic regression robust against the potential violation of the linearity of logit assumption?   2. Given my above example, does it seem acceptable to include the original continuous variable in the model?   3. Are there any references or guides out there for recommending when it is satisfactory to accept that the model is robust against the potential violation of linearity of the logit?


title: What is the equivalent of a repeated ANOVA using a mixed effect model with lme?

text: I have a data frame with this structure:               'data.frame':   39 obs. of  3 variables:      $ topic   : Factor w/ 13 levels "Acido Folico",..: 1 2 3 4 5 6 7 8 9 10 ...          $ variable: Factor w/ 3 levels "Both","Preconception",..: 1 1 1 1 1 1 1 1 1 1 ...      $ value   : int  14 1 36 17 5 9 19 9 19 25 ...      and i want to test the effect value ~ variable, considering that observation are grouped in topics. So I thought to use a repeated measure ANOVA, where "variable" is considered as a repeatead measure on every topic. the call is aov(value ~ variable + Error(topic/variable)). My question is: if i want to do the same with lme how should I model the random effect?               ~ 1 | topic/variable      or               ~ variable | topic      the anova returns a significant F. for the lme, the first version gives a non significant effect, while the second is significant; so I suppose that the anova is more similar to the second. Did I understood right? Thanks


title: Why is the expectation maximization algorithm used?

text: From what little I know the EM algorithm can be used to find the maximum likelihood when setting to zero the partial derivatives with respect to the parameters of the likelihood gives a set of equations that cannot be solved analytically. But is the EM algorithm needed instead of using some numerical technique to try to find a maximum of the likelihood with respect to the constraint of the set of equations mentioned.


title: ANOVA and Levene test problem (2)

text: I want to reword my previous question. I'm running ANOVA for 4 clusters to understand is there any difference between these clusters in terms of their reasons for travel. After I checked Levene test, one of my variables had a p-value lower than 0.05, but that variable has F(3,324) = 6.35, p=0.000 for the ANOVA. Now, I have 2 questions :   1. When the Levene test is rejected, but on the other hand, the ANOVA test is also rejected, i.e. there is a statistical difference, can we accept this ANOVA analysis?   2. When we used ANOVA for all of the variables and the Levene test is not rejected for all of them, with regard to this fact that the Levene test of one variable is rejected, can we use Kruskal-Wallis test just for that problematic variable or not? in other words, is it true that I just use a Kruskal-Wallis test for one variable?


title: Suitable graph to visualize the spread of data

text: I have some data like               A  2     A  4     A  76     B  8     B  13     .. ....      Basically, what is the spread of A, B, etc. What would be a suitable graph to visualize such information in R?


title: Difference of means -- Testing two algorithms

text: I have to perform 2 tests of the difference of means for the results obtained by 2 learning algorithms. **TEST 1**   1. I have 2 algorithms   2. I take a data set (with N examples)   3. I execute both algorithms 10 times   4. Each time, each algorithm is trained on a set of 9/10*N examples and tested against a set of N/10 examples (i.e. a 10-fold cross validation is done), and the result is a XX% value of correct classified examples in such a test set. Here I use the t-test. I know that the t-test makes assumptions of normality. I use it because the data set is quite big ( N>100 ), so I can assume a normal distribution of data (and if data is normally distributed, the means of a single algorithm are normally distributed, and so the differences of the means of the two algorithsm are normally distributed). Or should I consider only N/10 as the sample size from which the means are drawn? **TEST 2**   1. I have two algorithms   2. I take 10 instances of a combinatorial optimization problem (a Travelling Salesman Problem)   3. I run each algorithm on the 10 instances   4. For each run, a write down the result value (that is at some distance from the optimal value for that instance) How can I perform a test on the differences of means? Which means can I consider? Which data should I check for normality?


title: SNEP distribution fitting and parameter estimates

text: I'm looking to try and fit a SNEP distrib to some data but am not sure how to go about this - can anyone give me some derivations for parameter estimates etc?


title: Random Forest: IncNodePurity and Feature Selection for Binary Logistic Regression

text: After creating a Random Forest object using `randomForest` with around 500 candidate variables, I used `importance(object)` to display _IncNodePurity_ for each of the candidate variables in relation to the binary outcome of interest (Payment/No Payment). I am aware that _IncNodePurity_ is the total decrease in node impurities, measured by the Gini Index from splitting on the variable, averaged over all trees. What I don't know is what should be the cutoff for candidate variables to be retained after making use of `randomForest` for feature selection in regards to binary logistic regression models. For example, the smallest _IncNodePurity_ among my 498 variables is 0.03, whereas the largest _IncNodePurity_ is 96.68. In summary, I have one main question: Is there a cutoff for _IncNodePurity_? If yes, what is it? If no, how do you determine the cutoff? Do you simply take 10 candidate variables with the largest _IncNodePurity_ if you want a model with only 10 predictor variables? Any thoughts or references are greatly appreciated. Thanks!


title: What's the difference between Normalization and Standardization?

text: At work we were discussing this as my boss has never heard of normalization. In Linear Algebra, Normalization seems to refer to the dividing of a vector by its length. And in statistics, Standardization seems to refer to the subtraction of a mean then dividing by its SD. But they seem interchangeable with other possibilities as well. When creating some kind of universal score, that makes up 2 different metrics, which have different means and different SD's, would you Normalize, Standardize, or something else? One person told me it's just a matter of taking each metric and dividing them by their SD, individually. Then summing the two. And that will result in a universal score that can be used to judge both metrics. For instance, say you had the number of people who take the subway to work (in NYC) and the number of people who drove to work (in NYC). Train ---> x Car ---> y If you wanted to create a universal score to quickly report traffic fluctuations, you can't just add mean(x) and mean(y) because there will be a LOT more people who ride the train. There's 8 million people living in NYC, plus tourists. That's millions of people taking the train everyday verse hundreds of thousands of people in cars. So they need to be transformed to a similar scale in order to be compared. If mean(x) = 8,000,000 and mean (y) = 800,000 Would you normalize x & y then sum? Would you standardize x & y then sum? Or would you divide each by their respective SD then sum? In order to get to a number that when fluctuates, represents total traffic fluctuations. Any article or chapters of books for reference would be much appreciated. THANKS! Also here's another example of what I'm trying to do. Imagine you're a college dean, and you're discussing admission requirements. You may want students with at least a certain GPA and a certain test score. It'd be nice if they were both on the same scale because then you could just add the two together and say, "anyone with at least a 7.0 can get admitted." That way, if a prospective student has a 4.0 GPA, they could get as low as a 3.0 test score and still get admitted. Inversely, if someone had a 3.0 GPA, they could still get admitted with a 4.0 test score. But it's not like that. The ACT is on a 36 point scale and most GPA's are on 4.0 (some are 4.3, yes annoying). Since I can't just add an ACT and GPA to get some kind of universal score, how can I transform them so they can be added, thus creating a universal admission score. And then as a Dean, I could just automatically accept anyone with a score above a certain threshold. Or even automatically accept everyone whose score is within the top 95%.... those sorts of things. Would that be normalization? standardization? or just dividing each by their SD then summing?


title: Sample size required for mixed design ANOVA to achieve adequate statistical power

text: ### Context I ran an experiment with `3 x 2` design with three levels of within subjects factor (repeated measures) and two levels to the between subjects factors. I am interested in examining the changes from baseline and the interaction effect. ### Question   * How do I compute the required sample size for a `3 x 2` design in order to achieve adequate statistical power?


title: Machine learning algorithm for ranking

text: I have got a set of elements $X$ which I can describe according to $n$ characteristics. Thus: $$x_i: \\{c_{i1}, c_{i2}, \ldots, c_{in}\\} \mid x_i \in X $$ where $c_{ij}$ is the (numerical) evaluation for element $i$ according to the characteristics $j$. So my elements can be viewed as points in a $n$ dimension space. According to my readings, there exists algorithms like the "Bayes Classifier" which can provide me with a "yes" or "no" type of answer on any element of my set, provided that I did use a "training set" consisting of some elements of my set, and the expected outcome of the algorithm. Based on that data, the algorithm should be able to take any other element, not part of the training set, and provide a "yes" or "no" answer based on what it learnt thanks to the training set. This is great if you have some kind of idea of what your expect (training set) but you are unsure of the specific rules that yields to that result. What I would like to do with my data is not get a "yes" or "no" type of answer, but I would like to introduce a ranking within the elements. Some of them are "better" than others. Just like for the Bayes filter, I have got a generic idea of what I expect. I could thus generate a "training ranking" taken from a subset of my elements, that I would feed into the MLA. Based on that training it would be able to rank my whole set. In order to do that I see two approach :   1. Each element would be given a score by the MLA, then rank the elements according to the score.   2. The MLA would be able to take two elements $x_i$ and $x_j$ and determine wich one is better (pairwise comparisons). Use quicksort using that comparison operation. Note : based on a score the pairwise function is trivial to implement, and based on a pairwise function, it is trivial to generate a score, so these are just two approches to yield to the same results. Are there examples out there of MLA's that can provide a scoring function, or a pairwise comparison function ? EDIT : In order to add more context : currently my items are ranked according to an algorithm that generates a score (real number) te each items by making calculations on the $c_{ij}$. While the generated ranking is quite correct, I often have to modify the algorithm to tweak it in some way because I can clearly see some items that are not ranked according to what I would have expected. So currently my design process is :   1. Get an idea of what a perfect ranking would be   2. Try to (manually) derive an algorithm that would rank the items like that   3. Observe the results   4. Adapt the algorithm So I thought about MLA's since the starting point of my process is what could be used as the training data. I would probably start by taking the current rank, swap items according to my needs and feed that.


title: Estimate error of prediction from R-square

text: What I have:   1. a linear model $y=a_0+a_1x$ with given parameter estimates,   2. the number of values used for fitting the model,   3. the Pearson R² value. I need to estimate errors of prediction. I don't see a way to calculate it, but is there a way to at least get a rough estimate?


title: How to plot decision boundary of a k-nearest neighbor classifier from Elements of Statistical Learning?

text: I want to generate the plot described in the book ElemStatLearn "The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Second Edition" by Trevor Hastie & Robert Tibshirani& Jerome Friedman. The plot is: ![enter image description here](http://i.stack.imgur.com/oY7hr.png) I am wondering how I can produce this exact graph in `R`, particularly note the grid graphics and calculation to show the boundary.


title: Forecast pricing data with lead times?

text: I have some "big data" containing event pricing information with 45 days of lead time quotes up to the day before the event. Thus, it's structured like so: -45 100.00 -44 120.00 ... -1 110.00 We have events every day of the year, and have year-over-year data. I'm a relative newbie to forecasting models, so this is a basic question: is it possible to use this data to perform forecasting? What would be a good model to use to do this?


title: What is the difference between the sum of two covariance matrices and the covariance matrix of the sum of two variables?

text: I'm wondering if someone could help to explain the difference between two covariance matrices. Suppose that ${\bf K}_X$ and ${\bf K}_Y$ are two covariance matrices of real random vectors. What is the difference between ${\bf K}_X+{\bf K}_Y$ and ${\bf K}_{X+Y}$?


title: Problem with comparing GLM models having a different link function

text: Given the same set of covariates and distribution family, how can I compare models having different link functions? I think the correct answer here is "AIC/BIC", but I am not 100% sure. Is it possible to have nested models if they have a different link?


title: What to do when sample covariance matrix is not invertible?

text: I am working on some clustering techniques, where for a given cluster of d-dimension vectors I assume a multivariate normal distribution and calculate the sample d-dimensional mean vector and the sample covariance matrix. Then when trying to decide if a new, unseen, d-dimensional vector belongs to this cluster I am checking its distance via this measure: $$\left(X_i-\hat{\mu}_X\right)'\hat{\sigma}_X^{-1}\left(X_i-\hat{\mu}_X\right)>B_{0.95}\left(\frac{p}{2},\frac{-p}{2}\right)$$ Which requires me to calculate the inverse of the covariance matrix $\hat{\sigma}_X$. But given some samples I cannot guarantee that the covariance matrix will be invertible, what should I do in the case that it is not? Thanks