title: Combining p-values from different statistical tests applied on the same data

text: Although the title of the question appears trivial, I would like to explain that it is not that trivial in the sense that it's different from the question of applying the same statistical test in similar datasets to test against a total null hypothesis (meta-analysis, e.g. using Fisher's method for combining p-values). What I am looking for, is a method (if it exists and if the question is valid in statistical terms) that would combine p-values from two different statistical tests (e.g. a t-test and a u-test, even if one is parametric and the other is not), applied to compare the centers of two samplings from two populations. So far I have searched a lot on the web without a clear answer. The best answer I could find was based on game theory concepts by David Bickel (http://arxiv.org/pdf/1111.6174.pdf). A very simplistic solution would be a voting scheme. Suppose that I have two vectors of observations $A=[a_1, a_2, ..., a_n]$ and $B=[b_1, b_2, ..., b_n]$ and I want to apply several t-like statistics (t-test, u-test, even 1-way ANOVA) to test the hypothesis that the centers (means, medians etc.) of the two undlerlying distributions are equal against the hypothesis that they are not, at a significance level of 0.05. Suppose that I run 5 tests. Would it be legitimate to say that there is sufficient evidence to reject the null distribution if I have a p-value < 0.05 in 3 out of 5 tests? Would another solution be to use the law of total probability or this is completely wrong? For example, suppose that $A$ is the event that the null distribution is rejected. Then, using 3 tests, $T_1$, $T_2$, $T_3$ (meaning that $P(T_1)=P(T_2)=P(T_3)=1/3$), would a possible value for $P(A)$ be $P(A) = P(A|T_1)P(T_1) + P(A|T_2)P(T_2) + P(A|T_3)P(T_3)$, where $P(A|T_i)$ is the probability that the null distribution is rejected under the test $T_i$. I apologize if the answer is obvious or the question too stupid


title: Estimated effects may be unbalanced. What does it mean?

text: I am running a summary(aov(…)) in R and I got this message: `Estimated effects may be unbalanced` What does it mean? How may I solve this problem?


title: Confusion: different definitions of MAP estimation in Graphical Models (MRFs)

text: 1. The "classical" MAP estimation: $$\hat\theta = \arg\max_{\theta}P(\theta|\mathbf{x})$$ where $\mathbf{x}$ are the observations and $\theta$ are the parameters.    2. In this book chapter (page 6, second item), MAP estimation for a MRF is to maximize $P(\mathbf{x}|\mathbf{z},\theta)$ w.r.t $\mathbf{x}$, where $\mathbf{x}$ is the sequence of states, $\mathbf{z}$ is the set of observed data.   3. In this paper (page 2), that is to maximize $P(\mathbf{x}|\theta)$.  I would appreciate it if somebody could make the connection among these three clearer. Thanks in advance.


title: What is the optimal sample size for fitting a GARCH model?

text: I tried fitting an ARMA(1,1)/GARCH(1,1) model to my data consisting of around 5000 data points but I got significant results in Ljung Box test on standardized residuals and squared residuals. However when I used only the last 3000 data points the model showed much better results with non- significant standardized residuals and squared residuals. My question is why is this the case?Isn't more data supposed to give better models?If not what is the optimal sample size? Also please see my unanswered question: Procedure for fitting an ARMA/GARCH Model


title: How can I classify images as "normal" or "lesions" using SVM?

text: I want to classify images as "normal" or "lesions", and then rate the images in the lesions category from mild and severe. This classification would be based on edge, color and texture features. For this concept I want to use a binary SVM classifier using threefold cross validation, in Matlab. Where do I start?


title: How to encode factors as dummy variables when using stepPlr?

text: When using the `step.plr()` function in the stepPlr package, if my predictors are factors, do I need to encode my predictors as dummy variables manually before passing it to the function? I do know that I can specify "level", but how the "level" parameter works is confusing to me. My understanding is that I need to tell `step.plr()` explicitly which factors should be encoded as dummy variables and thus leaving one factor out intentionally. Let's consider a simple example. Suppose I have 1 categorical predicator with 4 levels and binary response. Normally, if I use `glm()` to fit a logistic regression model, `glm()` would automatically convert the categorical predicator into 3 dummy variables. Now in `stepPlr()`, do I specify the "level" parameter for that predictor with 4 levels or 3 levels? The "Help" section is vague, and says: > If the j-th column of x is discrete, level[[ j ]] is the set of levels for > the categorical factor. Does it mean I should tell `step.plr()` about all 4 levels, or I should make an intelligent decision myself and tell `step.plr()` to use only 3 levels? ==============UPDATE (16 Oct 2012)============= The following example will demonstrate what is the problem with `step.plr()`'s automatic dummy variable encoding. It is a slight modification of the code in the function's help section. set.seed(100)               n <- 100     p <- 3     z <- matrix(sample(seq(3),n*p,replace=TRUE),nrow=n)     x <- data.frame(x1=factor(z[ ,1]),x2=factor(z[ ,2]),                     x3=factor(sample(seq(3), n, replace=TRUE, prob=c(0.2, 0.5, 0.3))),                     x4=factor(sample(seq(3), n, replace=TRUE, prob=c(0.1, 0.3, 0.6))))     y <- sample(c(0,1),n,replace=TRUE)     fit <- step.plr(x,y, cp="aic")     summary(fit)      And here's an excerpt of the result:               Call:     plr(x = ix0, y = y, weights = weights, offset.subset = offset.subset,          offset.coefficients = offset.coefficients, lambda = lambda,          cp = cp)          Coefficients:           Estimate Std.Error z value Pr(>|z|)     Intercept  0.91386   5.04780   0.181    0.856     x4.1       1.33787   4.61089   0.290    0.772     x4.2      -1.70462   4.91240  -0.347    0.729     x4.3       0.36675   3.18857   0.115    0.908     x3.1:x4.1  7.04901  14.35112   0.491    0.623     x3.1:x4.2 -5.50973  15.53674  -0.355    0.723     x3.1:x4.3 -0.50012   7.95651  -0.063    0.950      You can see that all levels, that is, (1,2,3), are used to fit the model. But normally you only need two dummy variables to encode a predictor with 3 levels. On the other hand, if you use `glm()`:               glm(y~.^2, data=x, family=binomial)      you will get the correct dummy variable encoding.


title: What is a second order exponential decay function?

text: Can you explain what is a second order exponential decay function: $$ y(x) = y_0+A_{1}e^{-\frac{x}{t_1}}+A_{2}e^{-\frac{x}{t_2}} $$ (the $t_i$, $A_i$, and $y_0$ are constants and, presumably, the "decay constants" $t_i$ are positive)?   * Qualitatively, what is the difference between "first order" and "second order"? (A first order exponential function has the form $y(t)=y_0 + A_1 e^{-\frac{x}{t}}$.)   * How can we estimate $t_1$ and $t_2$ from data?


title: Do you have a global vision on those analysis techniques?

text: I'm currently on a project where I basically need, like we all do, to understand how output $y$ is related to input $x$. The particularity here is that data $(y,x)$ is given to me one piece at a time, so I want to update my analysis every time I receive a new $(y,x)$. I believe this is called "on- line" processing, as opposed to "batch" processing, where you have all the data needed and do your calculations using all data at the same time. So I've looked around for ideas, and I've finally came with the conclusion that the world is divided in three:   * The first part is the land of statistics and econometrics. People there do OLS, GLS, instrument variables, ARIMA, tests, difference of differences, PCA and whatnot. This land is mostly dominated by linearity and does only "batch" processing.   * The second part is the island of machine learning and other words like artificial intelligence, supervised and unsupervised learning, neural networks and SVMs. Both "batch" and "on-line" processing are done here.   * The third part is a whole continent I've just discovered, mostly populated by electrical engineers, so it seems. There, people often add the word "filter" to their tools, and they invented great stuffs like the Widrow-Hoff algorithm, Recursive least squares, the Wiener filter, the Kalman filter, and probably other things I haven't discovered yet. Apparently they do mostly "on-line" processing as it better fits their needs. So my question is, do you have a global vision on all this? I'm under the impression that these three parts of the world don't talk too much to each other. Am I wrong? Is there a Grand Unified Theory of Understanding How $Y$ Relates to $X$? Do you know any resources where the bases of that theory might be laid down? I'm not sure if this question really makes sense, but I'm a little lost between all those theories. I imagine the answer to the question "should I use this or that?" would be "it depends on what you want to do (and on your data)". However I feel like those three worlds try to answer to the same question ($y=f(x)$?) and so it should be possible to have a higher view on all this, and deeply understand what makes each technique particular.


title: Levels of variable and $\chi^2$

text: Does it matter how many levels the independent variable has for significance? For example, suppose socioeconomic status has 10 levels and the output variable has 2 levels. Suppose this is significant in a $\chi^2$ test. If we increase the number of levels to 20 for the independent variable would this change the significance? The reason I ask is because I can't seem to put my independent variable into 10 levels. Instead there are 30 levels (basically a listing of all the values). Would this matter for significance?


title: Visualizing multiple tables

text: In four different experiments, I ran two programs `A` and `B` with different parameter setteings `a`, `b`, and `c`. Target variable is kind of time -- lower Target value means higher performance. I know the reasons why program `A` is better than `B` or vice-versa. Currently I am representing this data as tables in my latex document. However, I am not happy with this. I would like to provide a single figure which could show the difference between the parameter setteings effectively, and would be easy for the reades to get an idea of why `A` is better than `B` or vice-versa in a particular experiment. Moreover, the units of the parameters `a`, `b`, and `c` are not the same. Thats why I don't like to use barplots for this problem. It would be great if some one help for drawing a nice graph in R. Thank you.               Experiment-1     program & a & b & c & Target     A & 0.86  & 151 & 7 & 27      B & 0.82 & 168 & 114 & 48           Experiment-2     program & a & b & c & Target     A & 0.88 & 23530 & 67 & 154      B & 0.79 & 23090 & 34 & 82           Experiment-3     program & a & b & c & Target     A & 0.19  & 7855 & 172 & 53     B & 0.16 & 7398 & 189 & 104              Experiment-4     program & a & b & c & Target     A & 0.30  & 5062 &301&  250      B & 0.21 & 8797 & 443 & 198


title: If correlation doesn't imply causality then what's the value of knowing the correlation between two variables?

text: Let's say as a business owner (or marketing or anyone who understands a scatter plot) is shown a scatter plot of two variables: number of advertisements vs number of product sales per month for the past 5 years (or another time-scale so that you have more samples. I just made this one up). Now he/she see's the scatter plot and is told that the correlation coefficient (corr) is:   1. 1 or   2. 0.5 or   3. 0.11 or   4. 0 or   5. -0.75 or   6. -1 Basically any valid value for `corr` Question: What does this even mean to a decision maker or any consumer of the scatter plot? What decisions can one take just based on this? I.e.: What is the use of seeing correlation between any two variables and what can one do with that information in isolation? Is it only to see what to and not to consider for inclusion in regression analysis or is there a more practical use? Just curious, I've always worked with this technique, but I've been told that correlation by itself is not of much use - so what "IS" the use?


title: Comparing two pearson coefficents w/ r^2

text: I have three variables, x,y, z. My sample n is about 17,000. My pearson coeficients & r^2: for x -> z it is 0.187, 0.034 for y -> z it's 0.311, 0.096 testing for significance, I receive a z-score of 12.3 & p < 0.00 So, I know the difference is significant, but I want to say in (somewhat layman's terms) how strong it is. So, I have looked at coefficient of determination, r^2, to get the ratio: 0.096:0.034, or 2.8. Is it then correct/accurate to say the correlation between y & z is 2.8 times as strong as the correlation between x & z?


title: Looking for formalization of the idea of low-variance predictors

text: In baseball, Bill James suggested using, to predict next season's winning percentage `runs_scored^2/(runs_scored^2 + runs_allowed^2)`, rather than this season's winning percentage, because there is a great deal of luck involved in winning games. Similarly, I think the number of points a football (American) team may be expected to allow is better predicted with at what rate they have given up yards, rather than at what rate they have given up points, because the former is subject to less luck.   My question is how would one go about building a predictive model that would select the low variance predictors, rather than the high-variance ones, which tend to do better on the training set? Is there anything better than just training a number of models and testing their performance on independent data?


title: Classification vs clustering

text: I am beginner to data mining. This is my understanding, regression is used to predict continuous values. It is a type of classification. Classification is Supervised and Clustering is unsupervised. In classification we have predefined classes, whereas in clustering there are no classes defined. What does this exactly mean? **But in both of them we do have training data, which tells which output the set of features should belong to. Then how can we say that clustering does not have per-defined classes?** Please explain when to use classification and clustering. Can clustering also predict continuous values? Thanks. Please help. This post has been edited.


title: How is AUC of decision tree calculated?

text: I have a dataset which only has one continuous variable, and I try to use decision tree algorithm to build a model which classify the +ve and -ve label from the dataset. I run 10-fold cross-validation. How is the AUC calculated for the decision tree classifier? Will the algorithm check different threshold value of the classifier, and determine the AUC? What about I have more than 2 continuous variable?


title: RMSE for negative binomial hurdle models?

text: I am working in Program R. I am conducting analysis on zero inflated over dispersed data. I am leaning towards hurdle regression with the count aspect modeled in a negative binomial framework. I am interested in assessing the model fit and stumbled upon this paper which uses root mean square error (RMSE). I calculated the RMSE from my model using the following equation and was less than pleased with the results.               Hurdle.RMSE <- sqrt(mean(Hurdle.mod$residuals^2))      As I said this data is zero inflated (82% zeros in fact). When I calculate the number of predicted zeros in the data I get 465 (which is 82% as hurdle partitions out all zeros). I looked at the predict values from the model and non of them are estimated as zero. So after looking into that further I find that the predicted values are the mean value for the particular covariate environment associated with that observation as specified by the coefficients of the model. I interpret this as the model accounting for the over abundance of zeros. So when comparing the raw residuals (predicted - observed) from a zero-inflated model one would expect the residuals to be large as the majority of the observed values are zero and there is a positive estimate for every predicted value. I thought the model was performing well as the predicted count from the model is fairly close to the mean of the original observations (matching the covariates of course). So I am now under the impression that I am either 1. mis-calculating the RMSE for this type of model, or 2. RMSE should not be used with this type of model. I am open to suggestions on other types of model fit diagnostics but other correlative assessments (Pearson's r, Spearman's p) should also indicate a lack of fit if I use the predicted values as they are spit out from the hurdle framework. Thanks.


title: K-medians, formula to compute the median

text: If you are running K-medians, and your distance metric is the L1 norm, how do you derive that the center of each centroid is the median of the data points assigned to it? Second, how do you compute the geometric median? Third, are there any implementations of k-medians algorithm?


title: Getting the rule from cross validation

text: I've got a question. Let's say I have a medical data representing 2 classes of patients (healthy and unhealthy) and some number of predictors which characterize these patients. Choosing different combinations of predictors I may construct logistic regression model (for distinguishing healthy and unhealthy patients) and then compare prognosis with these combination based on some particular measure, Area under roc curve, for example. In order the results to be valid I need to somehow split my data into training and test sets. Let's say that my data is of a small size and therefore I want to use some leave-k-out cross validation. This allow me to compare different combinations of predictors according to their prognostic possibilities. The question is: how can I get a classification rule from cross validation? If I have 1 particular split then everything is obvious: based on training set I have a particular rule which is then applied to a test set. But how should one deal in the situation of leave-k-out cross validation where there are many "trainings" each with different coefficients in the logistic regression. Should this coefficients be averaged or something else? Is there some substantiated procedure for getting the rule from cross validation? Thanks!


title: Calculating if population mean differs from a given value

text: I have an exam in two days and can't figure this out at all: > In the local cola canning factory, the mean fill of cans is set at 300 > milliliters but there is concern that the population mean fill may not in > fact be 300ml. Assume that the standard deviation of the amount of cola in a > random can is 1. 2 ml. A random sample of 100 cans showed a mean fill of x̄ >  299. 64 ml. > > Is there evidence at the 1% significance level that the population mean fill > differs from 300 ml? To answer this question, carry out in detail a z-test > but do not use the p-value method. > > Calculate the p-value of the test in (i) above and explain why your decision > would be the same as that which you reached in (i) if you had been asked to > conduct the test in (i) by the p-value method.


title: Standard normal distribution on a subspace

text: Let $U \subset \mathbb{R}^n$ be a vector space with $\dim(U)=d$. A standard normal distribution on $U$ is the law of a random vector $X=(X_1, \ldots, X_n)$ taking values in $U$ and such that the coordinates of $X$ in one ($\iff$ in any) orthonormal basis of $U$ is a random vector made of $d$ independent standard normal distributions ${\cal N}(0, 1)$. When reading this question I asked myself the following question. Let $Y=(Y_1, \ldots, Y_n)$ be a standard normal distribution on $\mathbb{R}^n$. Is is true that the conditional distribution of $Y$ given $Y \in U$ is the standard normal distribution on $U$ ? The squared norm ${\Vert X \Vert}^2$ of $X$ has a chi-square distribution $\chi^2_d$. Thus, if this is true, that would explain @Argha's claim. Sorry if the LaTeX is mistyped, I don't see the LaTeX rendering :( **EDIT 01/10/2012:** Ok I see. Write $y=u+v$ the orthogonal decompostion of $y$ in $U\oplus U^\perp$. Then $$\Pr(Y\in \mathrm{d}y \cap Y \in U)=\Pr(P_U Y \in \mathrm{d}u)$$. That shows that $(Y \mid Y \in U) \sim P_U Y$. This is little bit heuristic but morally correct. Finally it is clear from the definition that $P_U Y$ is standard normal on $U$.


title: Kalman Filter in PROC UCM

text: Under what circumstances will the Kalman filter implementation in proc UCM fail? High variance of the constituents is one case. Are there other cases you've encountered? Anything other than :http://support.sas.com/documentation/cdl/en/etsug/63939/HTML/default/viewer.htm#etsug_ucm_sect031.htm and http://support.sas.com/documentation/cdl/en/etsug/63939/HTML/default/viewer.htm#etsug_ucm_sect029.htm


title: How to communicate regression model performance to non-stats people?

text: In R I created a gradient boosted random forest from 100,000 records with 10 cross validation folds using the `gbm` library. I want to communicate the strength and accuracy of the model, because it worked very well (far surpassing the random forest, kNN or SVM). There are plenty of strength and accuracy tests (RMS, MAE, etc.) but the one I favor in the work place is correlation, `cor()`, because it has proven to be the most intuitive for non-stats people. The correlation between the actual and predicted values was 0.80. I understand that correlation is a measure of linear association, but in this example, it is also a measure of predictive ability. If cor = 1, we are 100% accurate. But is it fair to say we are 80% accurate if cor = 0.80? Or if cor = -1, could we say it is -100% accurate; that doesn't seem right. **So my question is:**   * For predictive regression models, how are you communicating the strength and accuracy of the model to non-stats people?


title: How do I interpret the covariance matrix from a curve fit?

text: I'm not too great at statistics, so apologies if this is a simplistic question. I am fitting a curve to some data, and sometimes my data best fits a negative exponential in the form $a * e^{(-b * x)} + c$, and sometimes the fit is closer to $a * e^{(-b * x^2)} + c$. However, sometimes both of those fail, and I would like to fall back to a linear fit. My question is, how can I determine which model fits a particular data set the best from the resulting variance-covariance matrix that is returned from the `scipy.optimize.curve_fit()` function? I believe the variance is on one of the diagonals of this matrix, but I'm not sure how to interpret that. UPDATE: Based on a similar question, I'm hoping that the variance-covariance matrix can tell me which of the three models I am attempting best fits the data (I am trying to fit many datasets to one of these three models). The resulting matrices look like this for the given example:               pcov_lin      [[  2.02186921e-05  -2.02186920e-04]      [ -2.02186920e-04   2.76322124e-03]]     pcov_exp     [[  9.05390292e+00  -7.76201283e-02  -9.20475334e+00]      [ -7.76201283e-02   6.69727245e-04   7.90218415e-02]      [ -9.20475334e+00   7.90218415e-02   9.36160310e+00]]     pcov_exp_2      [[  1.38338049e-03  -7.39204594e-07  -7.81208814e-04]      [ -7.39204594e-07   8.99295434e-09   1.92970700e-06]      [ -7.81208814e-04   1.92970700e-06   9.14746758e-04]]      Here is an example of what I am doing:               import numpy as np     import matplotlib.pyplot as plt     import scipy as sp     import scipy.optimize          def exp_func(x, a, b, c):         return a * np.exp(-b * x) + c          def exp_squared_func(x, a, b, c):         return a * np.exp(-b * x*x*x) + c          def linear_func(x, a, b):         return a*x + b          def main():         x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], np.float)         y = np.array([1, 1, 1, 1, 0.805621, 0.798992, 0.84231, 0.728796, 0.819471, 0.570414, 0.355124, 0.276447, 0.159058, 0.0762189, 0.0167807, 0.0118647, 0.000319948, 0.00118267, 0, 0, 0], np.float)              p0 = [0.7746042467213462, 0.10347274384077858, -0.016253458007293588]         popt_lin, pcov_lin      = scipy.optimize.curve_fit(linear_func, x, y)         popt_exp, pcov_exp      = scipy.optimize.curve_fit(exp_func, x, y)         popt_exp_2, pcov_exp_2  = scipy.optimize.curve_fit(exp_squared_func, x, y)              plt.figure()         plt.plot(x, y, 'ko', label="Original data")         plt.plot(x, linear_func(x, *popt_lin), 'r-', label='linear')         plt.plot(x, exp_func(x, *popt_exp), 'b-', label='exponential')         plt.plot(x, exp_squared_func(x, *popt_exp_2), 'g-', label='exponential squared')         plt.legend()         plt.show()          if __name__ == '__main__':         main()


title: Getting different results when plotting 95% CI ellipses with ggplot or the ellipse package

text: I want to visualize the results of a clustering (produced with `protoclust{protoclust}`) by creating scater plots for each pair of variables used for classifying my data, colouring by classes and overlapping the ellipses for the 95% confidence interval for each of the classes (to inspect which elipses-classes overlap under each pair of variables). I have implemented the drawing of the ellipses in two different ways and the resulting ellipses are different! (bigger ellipses for first implementation!) A priori they differ only in size (some diferent scaling?), as the centers and angle of axes, seem to be the similar in both. I guess I must be doing something wrong by using one of them (hope not with both!), or with the arguments. Can anyone tell me what I am doing wrong? Here the code for the two implementations; both are based on the answers to How can a data ellipse be superimposed on a ggplot2 scatterplot?               ### 1st implementation      ### using ellipse{ellipse}     library(ellipse)     library(ggplot2)      library(RColorBrewer)     colorpal <- brewer.pal(10, "Paired")          x <- data$x         y <- data$y     group <- data$group     df <- data.frame(x=x, y=y, group=factor(group))          df_ell <- data.frame()      for(g in levels(df$group)){df_ell <- rbind(df_ell, cbind(as.data.frame(with(df[df$group==g,], ellipse(cor(x, y),scale=c(sd(x),sd(y)),centre=c(mean(x),mean(y))))),group=g))}           p1 <- ggplot(data=df, aes(x=x, y=y,colour=group)) + geom_point() +        geom_path(data=df_ell, aes(x=x, y=y,colour=group))+scale_colour_manual(values=colorpal)          ### 2nd implementation      ###using function ellipse_stat()      ###code by Josef Fruehwald available in: https://github.com/JoFrhwld/FAAV/blob/master/r/stat-ellipse.R          p2 <-qplot(data=df, x=x,y=y,colour=group)+stat_ellipse(level=0.95)+scale_colour_manual(values=colorpal)      Here is the two plots together (left graph is `p1` implementation (`ellipse()`): ![enter image description here](http://i.stack.imgur.com/Qf49A.jpg) The data are available here: https://www.dropbox.com/sh/xa8xrisa4sfxyj0/l5zaGQmXJt


title: Relative variances of higher-order vs. lower-order random terms in mixed models

text: _**TL, DR summary:_** Is there any theoretical or empirical basis to support the following statement being true as a general rule of thumb? "When estimating a mixed model, typically the estimated variances/standard deviations of random effects associated with 'higher-order' terms (e.g., random effects of two-way, three-way, and beyond interaction terms) turn out to be _smaller_ than the estimated variances/standard deviations of random effects associated with 'lower-order' terms (e.g., the residual variance, variances associated with simple effects of grouping factors)." The source of this claim is me. ;) * * * Okay, now for the longer version... Typically when I sit down to start analyzing a new dataset which I know will call for a mixed model, one of the first models that I fit (after the statistical foreplay of looking through the observations in the dataset, plotting various things, cross-tabulating different factors, etc.) is one that is pretty close to the "maximal" random effects specification, where every random effect that is in-principle possible to estimate from the data, is estimated. Naturally, it is not uncommon that this nearly-maximal model will have some computational problems (convergence errors, or wacky variance/covariance estimates, or etc.) and that I have to trim back this model to find one that my data can more easily support. Fine. In these situations, the method I have come to prefer for trimming random terms is not to rely on significance tests or likelihood ratios, but rather to just identify the random effects that seem to have the smallest standard deviations (which can admittedly be a little tricky when predictors are on very different scales, but I try to take account of this in my appraisal) and remove these terms first, sequentially in an iterative process. The idea being that I want to alter the predictions of the model as little as possible while still reducing the complexity of the model. One pattern that I seem to have noticed after a pretty good amount of time spent doing this is that following this method very often leads me to trim random effects associated with higher-order terms (as defined above) of the model first. This is not always true, and occasionally some of the higher- order terms explain a lot of variance, but this doesn't seem to be the general pattern. In sharp contrast, I usually find that lower-order random terms -- particularly those associated with simple effects of the grouping factors -- explain a pretty good amount of variance and are fairly essential to the model. At the extreme, the residual term commonly accounts for close to the most variance, although of course removing this term wouldn't be sensible. This entirely informal observation leads me to form the hypothesis that I stated at the beginning of this question. If it is true, then it constitutes a useful piece of advice that might be passed down to people who are less experience with this kind of model selection process. But before I begin doing so, I want to check with other, more experienced users of mixed models about their reactions to this observation. Does it seem more or less true to you? Is it roughly consistent with your experience fitting many different mixed models to many different datasets? Do you know of any sensible, theoretical reasons why we might actually _expect_ this to be true in a lot of cases? Or does it just seem like bullshit? One possible answer here is that it is not true even in my own case, and I have simply deceived myself. Certainly a possibility that I am open to. Another possibility is that it might be true in my own case, but that this could simply be a kind of coincidence having to do with the kinds of datasets that I tend to work with routinely (which, FYI, are datasets in psychological / social sciences, a slight majority being experimental in origin, but also a fair proportion of non-experimental stuff). If this is the case then there is probably no good reason for expecting my observations to hold in general in other fields that handle very different kinds of data. Still, if there is a coherent non-coincidental reason for why this might be expected to be true, even if only for these particular kinds of datasets, I would love to hear it. And of course another possibility is that others _have_ noticed similar patterns in their own data, and that it represents some kind of general rule of thumb that people find useful to keep in mind as they fit mixed models to various different data. If this is the case then it seems like there must be some compelling statistical-theoretical reason for why this pattern arises. But I really don't know what that reason would look like. I welcome anyone's thoughts and opinions about this. Note that as far as I'm concerned, totally legitimate responses to this question might be as simple as comments like "Yeah I have noticed something similar in the data I've worked with, but I have no idea why it should be true" or conversely "I have noticed nothing remotely like this in the data I've worked with." Of course I also welcome longer and more involved discussions...


title: How do I interpret the credibility interval in a Bayesian Regularized Regression?

text: A penalized regression provides biased estimates of the regression coefficients (bias-variance trade-off principle). Therefore, standard errors and confidence intervals are regarded as not very meaningful for those biased estimates arising from (frequentist) penalized regression method, see e.g. the discussion Estimating R-squared and statistical significance from penalized regression model . I would assume that the same problems exists in an Bayesian approach but Kyung, Gill, Ghaosh and Casella (2010) say that the Bayesian formulation produces valid standard errors. Does it mean that a 95% credibility intervals includes with 95% probability the true biased estimate and if yes, is this a useful information?


title: What ways are there to show two analytical methods are equivalent?

text: I have two different analytical methods that can measure the concentration of a particular molecule in a matrix (for instance measure the amount of salt in water) The two methods are different, and each has it's own error. What ways exist to show the two methods are equivalent (or not). I'm thinking that plotting the results from a number of samples measured by both methods on a scatter graph is a good first step, but are there any good statistical methods ?


title: Modeling a Poisson process with 10000 events per minute

text: I am trying to model a system that generates events modeled by a Poisson process. I am using the following ruby code:               INTERVAL = 0.005     LAMBDA = 167.0     events = Hash.new(0)          def f(x, lambda)       1 - Math.exp(-lambda * x)     end          random_gen = Random.new     start = Time.now.to_f          while Time.now.to_i - start < 60       if random_gen.rand < f(INTERVAL, LAMBDA)         bucket = (Time.now.to_f - start).round.to_i         events[bucket] += 1       end            sleep(INTERVAL)     end      I am trying to generate about 10000 events per minute, or 167 events/second, so I am using $\lambda = 167$. I am using function `f` as $1 - e^{(-\lambda x)}$ where $x$ is the interval I am sleeping, I am using this function so that the inter-arrival time of events follows an exponential distribution. However, I am not getting the expected results, this code generates about 7194 events per minute, with a mean of about 117 events per second. I would expect this code to generate 10000 events per minute with an mean of about 167 events per second. What am I doing wrong? Thank you for your help. **UPDATE** Fixed typo with sampling time and added `random_gen` and `start` definitions


title: Is is correct to compare t-statistics of different pairs of cointegrated timeseries?

text: I am testing for cointegration all the pairs from a set of 100 stocks. I run an Ordinary Least Square Regression on each pair and then I test for the existence of unit roots in the residuals. I am using the ADF test with no trend and no intercept. I am surprised that I found a big amount of pairs which the test on the residual rejected the null hypothesis of having a unit root, with p-values smaller than 0.01. Can I use the Dick Fuller Statistic to rank those pairs? All the series have the same length. Thank you!


title: How to establish relationships using box plots?

text: I am examining the relationship between a categorical variable and a quantitative one. For the test of linearity, I was told that I must use scatter plots. However, this is no longer applicable if one of the variables is categorical. So I resorted to box plots. Now, I have difficulty identifying/checking if there is a linear relationship among males and females (gender is the categorical variable) in terms of their side by side box plots. Can anyone help me with this?


title: Forecasting hierarchical time series R package

text: I have to forecast a large set of (hierarchical) time series and since the R package hts allows for confidence intervals for their ensemble, I'd like to use it. I haven't found an example of how to use it, yet. How should I forecast the following series with it: The time series have a simple form. The sum of $n_3$ level-$3$-series make up a level-$2$-series and the sum of $n_2$ level-$2$-series make up a level-$1$-series and the sum of $n_1$ level-$1$-series make up the top-most level-$0$-series (where sums are taken over series in one point in time, $n_j$ natural numbers for $j=0,1,2,3$).


title: Look and you shall find (a correlation)

text: I have a several hundreds of measurements. Now I am considering utilizing some kind of software to correlate every measure with every measure. This means that there are hundreds of correlations. Among these there should (statistically) be a high correlation, even if the data is completely random (each measure has only about 100 datapoints). When I find a correlation, how do I include the information about how hard I looked for a correlation into it? I am not at a high level in statistics, so please bear with me.


title: When is a deviation statistically significant?

text: In a double-blind study, when are deviations from the control group considered statistically significant? And is this related to the number of samples? I realise that every experiment is different and statistical significance should depend on the deviations in measurements and the size of the sample group, but I'm hoping there is an intuitive rule-of-thumb formula that can "flag" an event as interesting. I have a light stats background in engineering, but am by no means a stats guru. A worked example would be much appreciated to help me understand and apply it to every-day things. * * * **Update with example:** OK, here's a (not so simple) thought experiment of what I mean. Suppose I want to measure the toxicity of additives in a village's water supply by comparing mortality rates over time. Given the village's population, natality and mortality rates over several years and the date upon which an additive was introduced into the water supply (disregard quantity), when would a rise in the mortality rate become interesting? Intuitively, if the mortality rate remains between 0.95% and 1.25% for 10 years, and suddenly spikes to 2.00%, then surely this would be an interesting event _if_ an additive was added that year (assume short-term toxic effects). Obviously there could be other explanations for the rise, but let's focus on statistical significance. Now, how about if it rises to 1.40%? Is that statistically significant? Where do you draw the line? I'm starting to get the feeling that I need to "choose" a critical region, which feel less authoritative. Can a Gaussian distribution guide me on this? What other information do I need to determine statistical significance?


title: How do you model a scalar response using a scalar IV, a categorical IV and categorical covariate w/ interactions?

text: So I originally thought I would be able to solve this problem fairly easily with ANCOVA, but now that I've thought more deeply about it I realized it is not so simple. Essentially we have a scalar response variable, $y$ which we want to predict using three variables $x$, $c$ and $g$. $x$ is a scalar variable (which takes three values -3, 0, 3), $c$ is a categorical variable (which takes three values C1, C2, and C3) and $g$ is a categorical variable which takes two values (Y or N). The issue is that $g$ was not an independent variable in the design of the experiment, but a covariate that was measured during the experiment. Further, we have reason to believe that this covariate may interact with the other two variables in how it predicts $y$. How do we design a model that is able to both consider the fact that $g$ is a covariate and allow us to analyze interactions between our three predictors. Ideally, if whoever is capable of answering this can provide matlab code to solving this problem, that would be wonderful. However, even a way to create the design matrix for a general linear model would be immensely helpful. Also, if I am approaching this entire problem from an incorrect perspective, I'm very open to other approaches to analyzing this data. Thank you!


title: Dirichlet distribution plot in R

text: I want to plot $15$ random draws of Dirichlet distribution with $\alpha = 1$ and dimension $n=10$ in R.               require(MCMCpack)     alpha <- 1     draws <- 15     dimen <- 10     x <- rdirichlet(draws, rep(alpha, dimen))      I want an output similar to the following image. This image is from Prof. David Blei's Topic Modeling tutorial at KDD 2011. I think it has been done using ggplot2 but I don't know how to generate it. Any help would be greatly appreciated. ![enter image description here](http://i.stack.imgur.com/2mScL.png)


title: Topic prediction using latent Dirichlet allocation

text: I have used LDA on a corpus of documents and found some topics. The output of my code is two matrices containing probabilities; one doc-topic probabilities and the other word-topic probabilities. But I actually don't know how to use these results to predict the topic of a new document. I am using Gibbs sampling. Does anyone know how? thanks


title: Error propagation and under-represented variable

text: I am running a relatively simple propagation of error analysis on an equation: $$S = \frac{ABC}{XY}$$ I have independent random error terms calculated for each variable (for the purposes of this analysis I am assuming systematic errors to be negligible), but I am having a specific issue with the $C$ variable in that, after some testing, we are certain that the quantity it represents is being under- represented. Restated, during the data collection process, some data loss is occurring and $C$ represents a numeric quantity and is capturing roughly only 75% of the true value according to our best estimates (which is to say, if we have, say, $C=5$, our data indicates that this is only 75% of the true value of C). Is there a specific way to deal with this situation? My first inclination is to run some simple algebra and calculate the true value for the variable (in the above example of $C=5$, quick math gives me a "true" value of roughly $C=6.6667$) but I am unsure if this is the best or most rigorous way to go about handling this situation. Any advice or citations to assist would be greatly appreciated.


title: Orthogonal contrasts in R

text: I am having somewhat of a problem setting up contrasts from an anova in R. I have done this in the past, but at the moment I don't seem to be able to get R to use my contrast matrix instead of the treatment contrasts. Any help on this would be greatly appreciated, I guess its possible I am making some basic mistake that I can not see. (I am using R2.15.1) My data consists of growth rate estimates for a 37 of bacterial clones. For these clones we measured growth rate multiple times (4-20) and we were wanting to simply know whether some of the the clones had a different growth rate (mu) from an ancestral clone (Anc). So I set up an anova as:               m1<-aov(mu~clone)      which gave me a significant effect of clone:                            Df Sum Sq Mean Sq F value   Pr(>F)         clone        36  32.52  0.9033   2.986 3.98e-06 ***     Residuals   123  37.21  0.3025                          ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       I then set up my contrast matrix               levels(clone)      [1] "A21"  "A22"  "A23"  "A24"  "A25"  "A26"  "A27"  "A28"  "A29"  "A30"      [11] "A31"  "A32"  "A33"  "A34"  "A35"  "A36"  "A37"  "A38"  "A39"  "A40"      [21] "Anc"  "T126" "T127" "T128" "T129" "T130" "T131" "T132" "T133" "T134"     [31] "T135" "T35"  "T36"  "T37"  "T38"  "T39"  "T40"       As I wanted to compare all other clones to the Ancestral (Anc), I created a matrix with 36 columns in which I assigned 1 to Anc, -1 to the clone I wanted to compare and 0 to all others. Thus the matrix was as:               contrasts(clone)<-cbind(     c(-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0)      ,c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1))      Now, the problem is that when I rerun the anova and do summary.lm, I still get the treatment contrasts being used and obtain:               Call:     aov(formula = mu ~ clone)          Residuals:           Min        1Q    Median        3Q       Max      -0.059077 -0.010799  0.000524  0.013324  0.057926           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept) 0.111175   0.011905   9.339 5.35e-16 ***     cloneA22    0.057165   0.016836   3.395 0.000923 ***     cloneA23    0.034798   0.016836   2.067 0.040843 *       cloneA24    0.009232   0.016836   0.548 0.584436         cloneA25    0.047776   0.016836   2.838 0.005316 **      cloneA26    0.013474   0.016836   0.800 0.425091         cloneA27    0.012300   0.016836   0.731 0.466445         cloneA28    0.047080   0.018185   2.589 0.010787 *       cloneA29    0.036722   0.018185   2.019 0.045626 *       cloneA30    0.035574   0.016836   2.113 0.036625 *       cloneA31    0.060275   0.016836   3.580 0.000493 ***     cloneA32    0.045941   0.016836   2.729 0.007290 **      cloneA33    0.017751   0.016836   1.054 0.293792         cloneA34    0.013581   0.016836   0.807 0.421418         cloneA35    0.018089   0.016836   1.074 0.284750         cloneA36    0.035200   0.016836   2.091 0.038607 *       cloneA37    0.048030   0.016836   2.853 0.005086 **      cloneA38    0.046683   0.016836   2.773 0.006423 **      cloneA39    0.040877   0.016836   2.428 0.016632 *       cloneA40    0.051376   0.016836   3.052 0.002790 **      cloneAnc    0.020183   0.013161   1.533 0.127727         cloneT126   0.020191   0.016836   1.199 0.232738         cloneT127   0.051375   0.016836   3.051 0.002790 **      cloneT128   0.043587   0.016836   2.589 0.010789 *       cloneT129   0.038895   0.016836   2.310 0.022543 *       cloneT130   0.017442   0.016836   1.036 0.302238         cloneT131   0.020714   0.016836   1.230 0.220928         cloneT132   0.038433   0.016836   2.283 0.024159 *       cloneT133   0.028807   0.016836   1.711 0.089601 .       cloneT134   0.037626   0.016836   2.235 0.027233 *       cloneT135   0.017300   0.016836   1.028 0.306186         cloneT35    0.050487   0.016836   2.999 0.003282 **      cloneT36    0.040415   0.016836   2.401 0.017871 *       cloneT37    0.054457   0.016836   3.235 0.001565 **      cloneT38    0.041835   0.016836   2.485 0.014304 *       cloneT39    0.038966   0.016836   2.314 0.022301 *       cloneT40    0.025847   0.016836   1.535 0.127292      So here clone A21 is still being used for comparison. Thanks in advance for any possible help, Ricardo


title: Is it advantageous to use dummy variables when learning a regression model?

text: Say I have 3 random variables ${X, Y, Z}$, and I have collected an iid sample of size $N$ from them: ${\cal D} = \{ (x_i, y_i, z_i), i = 1,\dots,N \}$. The conditional expectation $E[ X | y ]$ can also be written as \begin{eqnarray*} E[ X | y ] = \int E[ X | y, z ] \,\, dF_{ Z | Y }(z) \end{eqnarray*} My objective is to learn $E[ X | y ]$ from ${\cal D}$, and my question is: For finite sample size, should I expect any numerical advantage from learning $ E[ X | y, z ], dF_{ Z | Y }(z)$ first, and then integrating out the dependence on z? Or, can I discard the data collected about $z$, and directly learn $E[ X | y ]$. To make things simpler, I'm ready to assume that my learning strategy is consistent.


title: Choosing the right correlation coefficient

text: Let $X \in \{0,1\}^{n \times m}$ be a matrix with $n,m \geq 600$. I want to calculate the correlations between the columns of the matrix $X$. To be more precise, I need the correlation matrix $C \in [-1,1]^{m \times m}$. My question is, is a rank correlation coefficient like Spearman's $\rho$ or Kendall's $\tau_B$ suitable and if yes, which one more?


title: "Intervals" for a scale-based answer

text: I have to analyse data that look at the relevance of some posts. Basically each post is rated by 100 persons, and the scale looks something like this (like a Likert scale):   1. Irrelevant   2. Slightly relevant   3. Relevant   4. Useful   5. Extremely Useful But each person can choose a certain position with a 0.25 increment on this scale. So, more specifically, the scale looks like this [1, 1.25, 1.50, 1.75, 2, 2.25,...3,...4,...,4.75,5] I was asked to compute some intervals for each post. My first reaction was to compute the mean and the 95% CI, which is evidently wrong since the data is not normally distributed. Which is the best approach for this problem, where should I start looking? I read something about Bayesian Approach on this site, will that do? Thanks


title: How to interpret p-values of 0 or 1?

text: I ran an ANOVA finding for example an interaction between gender and grade than I want to know in what grades boys and girls differ, but in many cases I find (adjusted) p-values of 0 and 1. How / why is this possible? Doesn't seem right...               as.factor(gender)                     1     16    16.2    2.6377  0.104396         as.factor(grade)                      7  50077  7153.9 1165.4184 < 2.2e-16 ***     as.factor(gender):as.factor(grade)    7    132    18.9    3.0795  0.003056 **      Residuals                          7747  47555     6.1                             ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1             Tukey multiple comparisons of means         95% family-wise confidence level          Fit: aov(formula = rating ~ as.factor(gender) * as.factor(grade), data = users_c[users_c$grade %in% 1:8, ])          $`as.factor(gender)`                diff        lwr        upr     p adj     m-f -0.09135851 -0.2016276 0.01891058 0.1043964          $`as.factor(grade)`              diff        lwr       upr     p adj     2-1 0.3823566 -0.5454435  1.310157 0.9169296     3-1 1.9796023  1.1649854  2.794219 0.0000000     4-1 3.9558543  3.1534606  4.758248 0.0000000     5-1 5.7843111  4.9829529  6.585669 0.0000000     6-1 7.0752044  6.2708610  7.879548 0.0000000     7-1 8.4868609  7.6776332  9.296089 0.0000000     8-1 9.3867231  8.5626511 10.210795 0.0000000     3-2 1.5972457  1.0395026  2.154989 0.0000000     4-2 3.5734976  3.0337642  4.113231 0.0000000     5-2 5.4019544  4.8637616  5.940147 0.0000000     6-2 6.6928478  6.1502200  7.235476 0.0000000     7-2 8.1045042  7.5546625  8.654346 0.0000000     8-2 9.0043665  8.4329024  9.575831 0.0000000     4-3 1.9762520  1.6694948  2.283009 0.0000000     5-3 3.8047088  3.5006705  4.108747 0.0000000     6-3 5.0956021  4.7837806  5.407424 0.0000000     7-3 6.5072586  6.1830461  6.831471 0.0000000     8-3 7.4071208  7.0474558  7.766786 0.0000000     5-4 1.8284568  1.5588754  2.098038 0.0000000     6-4 3.1193501  2.8410202  3.397680 0.0000000     7-4 4.5310066  4.2388618  4.823151 0.0000000     8-4 5.4308688  5.0998193  5.761918 0.0000000     6-5 1.2908933  1.0155630  1.566224 0.0000000     7-5 2.7025498  2.4132612  2.991838 0.0000000     8-5 3.6024120  3.2738803  3.930944 0.0000000     7-6 1.4116565  1.1141985  1.709114 0.0000000     8-6 2.3115187  1.9757711  2.647266 0.0000000     8-7 0.8998622  0.5525763  1.247148 0.0000000          $`as.factor(gender):as.factor(grade)`                     diff         lwr        upr     p adj     m:1-f:1  0.005917865 -1.77842639  1.7902621 1.0000000     f:2-f:1  0.318074165 -1.28953805  1.9256864 0.9999988     m:2-f:1  0.442924925 -1.11597060  2.0018205 0.9998619     f:3-f:1  1.769000750  0.35262166  3.1853798 0.0020136     m:3-f:1  2.174229216  0.76569156  3.5827669 0.0000147     f:4-f:1  3.738998543  2.34268666  5.1353104 0.0000000     m:4-f:1  4.163719997  2.77146170  5.5559783 0.0000000     f:5-f:1  5.769586591  4.37599400  7.1631792 0.0000000     m:5-f:1  5.816721075  4.42497532  7.2084668 0.0000000     f:6-f:1  7.169439003  5.77317769  8.5657003 0.0000000     m:6-f:1  7.000924045  5.60308216  8.3987659 0.0000000     f:7-f:1  8.330142924  6.92683436  9.7334515 0.0000000     m:7-f:1  8.674488370  7.26930678 10.0796700 0.0000000     f:8-f:1  9.535307293  8.11198164 10.9586329 0.0000000     m:8-f:1  9.251081088  7.82191240 10.6802498 0.0000000     f:2-m:1  0.312156300 -1.12690148  1.7512141 0.9999959     m:2-m:1  0.437007060 -0.94741539  1.8214295 0.9995001     f:3-m:1  1.763082885  0.54136279  2.9848030 0.0000892     m:3-m:1  2.168311350  0.95569081  3.3809319 0.0000001     f:4-m:1  3.733080678  2.53468294  4.9314784 0.0000000     m:4-m:1  4.157802132  2.96412989  5.3514744 0.0000000     f:5-m:1  5.763668726  4.56844048  6.9588970 0.0000000     m:5-m:1  5.810803210  4.61772882  7.0038776 0.0000000     f:6-m:1  7.163521138  5.96518233  8.3618599 0.0000000     m:6-m:1  6.995006180  5.79482611  8.1951862 0.0000000     f:7-m:1  8.324225059  7.11768240  9.5307677 0.0000000     m:7-m:1  8.668570505  7.45984987  9.8772911 0.0000000     f:8-m:1  9.529389428  8.29962271 10.7591561 0.0000000     m:8-m:1  9.245163223  8.00863850 10.4816879 0.0000000     m:2-f:2  0.124850760 -1.02282435  1.2725259 1.0000000     f:3-f:2  1.450926585  0.50586965  2.3959835 0.0000172     m:3-f:2  1.856155050  0.92289131  2.7894188 0.0000000     f:4-f:2  3.420924378  2.50621691  4.3356318 0.0000000     m:4-f:2  3.845645832  2.93713824  4.7541534 0.0000000     f:5-f:2  5.451512425  4.54096139  6.3620635 0.0000000     m:5-f:2  5.498646910  4.59092496  6.4063689 0.0000000     f:6-f:2  6.851364838  5.93673457  7.7659951 0.0000000     m:6-f:2  6.682849880  5.76580854  7.5998912 0.0000000     f:7-f:2  8.012068759  7.08671595  8.9374216 0.0000000     m:7-f:2  8.356414205  7.42822339  9.2846050 0.0000000     f:8-f:2  9.217233128  8.26179669 10.1726696 0.0000000     m:8-f:2  8.933006923  7.96888762  9.8971262 0.0000000     f:3-m:2  1.326075825  0.46649985  2.1856518 0.0000150     m:3-m:2  1.731304290  0.88471145  2.5778971 0.0000000     f:4-m:2  3.296073618  2.46998162  4.1221656 0.0000000     m:4-m:2  3.720795071  2.90157332  4.5400168 0.0000000     f:5-m:2  5.326661665  4.50517434  6.1481490 0.0000000     m:5-m:2  5.373796150  4.55544575  6.1921465 0.0000000     f:6-m:2  6.726514078  5.90050756  7.5525206 0.0000000     m:6-m:2  6.557999120  5.72932364  7.3866746 0.0000000     f:7-m:2  7.887217999  7.04935402  8.7250820 0.0000000     m:7-m:2  8.231563445  7.39056617  9.0725607 0.0000000     f:8-m:2  9.092382368  8.22140761  9.9633571 0.0000000     m:8-m:2  8.808156163  7.92766524  9.6886471 0.0000000     m:3-f:3  0.405228465 -0.13578346  0.9462404 0.4221367     f:4-f:3  1.969997793  1.46166478  2.4783308 0.0000000     m:4-f:3  2.394719246  1.89762897  2.8918095 0.0000000     f:5-f:3  4.000585840  3.49977062  4.5014011 0.0000000     m:5-f:3  4.047720325  3.55206739  4.5433733 0.0000000     f:6-f:3  5.400438253  4.89224417  5.9086323 0.0000000     m:6-f:3  5.231923295  4.71940255  5.7444440 0.0000000     f:7-f:3  6.561142174  6.03389412  7.0883902 0.0000000     m:7-f:3  6.905487620  6.37327442  7.4377008 0.0000000     f:8-f:3  7.766306543  7.18788499  8.3447281 0.0000000     m:8-f:3  7.482080337  6.88942637  8.0747343 0.0000000     f:4-m:3  1.564769328  1.07871270  2.0508260 0.0000000     m:4-m:3  1.989490781  1.51520464  2.4637769 0.0000000     f:5-m:3  3.595357375  3.11716862  4.0735461 0.0000000     m:5-m:3  3.642491860  3.16971239  4.1152713 0.0000000     f:6-m:3  4.995209787  4.50929846  5.4811211 0.0000000     m:6-m:3  4.826694830  4.33626022  5.3171294 0.0000000     f:7-m:3  6.155913709  5.65010831  6.6617191 0.0000000     m:7-m:3  6.500259155  5.98928021  7.0112381 0.0000000     f:8-m:3  7.361078078  6.80213257  7.9200236 0.0000000     m:8-m:3  7.076851872  6.50319055  7.6505132 0.0000000     m:4-f:4  0.424721453 -0.01192015  0.8613631 0.0668946     f:5-f:4  2.030588047  1.58971048  2.4714656 0.0000000     m:5-f:4  2.077722532  1.64271796  2.5127271 0.0000000     f:6-f:4  3.430440460  2.98119847  3.8796825 0.0000000     m:6-f:4  3.261925502  2.80779484  3.7160562 0.0000000     f:7-f:4  4.591144381  4.12045589  5.0618329 0.0000000     m:7-f:4  4.935489827  4.45924616  5.4117335 0.0000000     f:8-f:4  5.796308750  5.26892973  6.3236878 0.0000000     m:8-f:4  5.512082545  4.96913148  6.0550336 0.0000000     f:5-m:4  1.605866594  1.17800058  2.0337326 0.0000000     m:5-m:4  1.653001078  1.23118920  2.0748130 0.0000000     f:6-m:4  3.005719006  2.56923916  3.4421989 0.0000000     m:6-m:4  2.837204048  2.39569420  3.2787139 0.0000000     f:7-m:4  4.166422928  3.70789927  4.6249466 0.0000000     m:7-m:4  4.510768373  4.04654394  4.9749928 0.0000000     f:8-m:4  5.371587296  4.85503631  5.8881383 0.0000000     m:8-m:4  5.087361091  4.55492128  5.6198009 0.0000000     m:5-f:5  0.047134485 -0.37906079  0.4733298 1.0000000     f:6-f:5  1.399852412  0.95913504  1.8405698 0.0000000     m:6-f:5  1.231337454  0.78563790  1.6770370 0.0000000     f:7-f:5  2.560556334  2.09799705  3.0231156 0.0000000     m:7-f:5  2.904901779  2.43669086  3.3731127 0.0000000     f:8-f:5  3.765720703  3.24558412  4.2858573 0.0000000     m:8-f:5  3.481494497  2.94557538  4.0174136 0.0000000     f:6-m:5  1.352717928  0.91787572  1.7875601 0.0000000     m:6-m:5  1.184202970  0.74431204  1.6240939 0.0000000     f:7-m:5  2.513421849  2.05645683  2.9703869 0.0000000     m:7-m:5  2.857767295  2.39508230  3.3204523 0.0000000     f:8-m:5  3.718586218  3.20341827  4.2337542 0.0000000     m:8-m:5  3.434360013  2.90326187  3.9654582 0.0000000     m:6-f:6 -0.168514958 -0.62249009  0.2854602 0.9968060     f:7-f:6  1.160703921  0.69016548  1.6312424 0.0000000     m:7-f:6  1.505049367  1.02895400  1.9811447 0.0000000     f:8-f:6  2.365868290  1.83862318  2.8931134 0.0000000     m:8-f:6  2.081642085  1.53882109  2.6244631 0.0000000     f:7-m:6  1.329218879  0.85401081  1.8044269 0.0000000     m:7-m:6  1.673564325  1.19285330  2.1542753 0.0000000     f:8-m:6  2.534383248  2.00296656  3.0657999 0.0000000     m:8-m:6  2.250157043  1.70328327  2.7970308 0.0000000     m:7-f:7  0.344345446 -0.15203755  0.8407284 0.5648416     f:8-f:7  1.205164369  0.65953016  1.7507986 0.0000000     m:8-f:7  0.920938164  0.36023867  1.4816377 0.0000022     f:8-m:7  0.860818923  0.31038540  1.4112524 0.0000101     m:8-m:7  0.576592718  0.01122178  1.1419637 0.0401330     m:8-f:8 -0.284226205 -0.89329509  0.3248427 0.9688007


title: Data sources related to Federal regulation

text: I'm looking for sources of data related to Federal regulation and enforcement. Things like the following:   * Size of Federal regulatory agencies (number of workers)   * Number of pages of Federal regulations [total pages that are in force; not just new ones]   * Number of enforcement actions brought by various agencies   * Cost of complying with regulations Can anyone tell me where I might find such information?


title: In natural language parsing, what is the feature function?

text: I'm working in the area of natural language processing, to be specific I'm reviewing the parsers that take advantage of data-mining techniques. I've read an introduction to natural language processing, specifically the maximum entropy approach by Della Pietra and others. In that paper, the authors explain that **features** are indicator functions that take into account certain aspects of the text. The problem is that this aspects affect the maximum entropy function. In the case of parsing, I have problems understanding the concept of features. For parsing one must have:   1. Two sets, $X$ for the sentences and $Y$ for the parses   2. For each sentence $x \in X$, one must have the N-best list of candidate parses; denoted $GEN(x) \subset Y$.   3. Training samples $(x_j,y_j)$. The $x_j$ are sentences and the $y_j$ are parses.   4. And a **feature** mapping $\Phi(X,Y)$ that maps each pair (x,y) to a vector of feature values. The problem of parsing will consist of finding the weight vector such that the following expression is minimized $$ -\sum log P(y_j|x_j,w) + \lambda ||w||_1 $$ where $$ P(y_j|x_j,w) = \frac{exp \{w^T \Phi(x_j,y_j)\} }{ \sum_{u \in GEN(x)} exp\{ w^T\Phi(x_j,u) \} } $$ My question is, what are the features in this case? Some codified version of the parsing tree? Are the features a statistical concept? I will also appreciate any bibliography that can be useful in tackling this issue. **UPDATE:** I've talked to a professor and he had just told me that the features are obtained after doing a classification algorithm based on the parse-tree. Do you know what kind of classification is done with the parse tree?


title: How to interpret coefficients in a Tobit regression?

text: I've run a Tobit regression but I'm not sure of how I should interpret the coefficients. Do you have any suggestions? Thank you very much!


title: C# ARMA library

text: Do you know any C# library or source code that can be used for ARMA / ARIMA forecasting?


title: Post hoc chi-square test

text: I am looking at a $\chi^2$ (crosstabs) test looking at interactions between types of dyads in an animal group. My categories are Male-Male, Male-Female, and Female-Female dyads in the rows and "Were Seen to Interact" and "Were Never Seen to Interact".   The matrix looks something like:               21   7       85   19       77   1      There is a very significant p value. My question is, when running the post hoc tests, is it appropriate to compare each category to the sum of the other two (i.e. Male-Male, vs Male-Female+Female-Female) when looking for which category is driving the difference, or should I compare the categories to each other (Male-Male vs Male-Female, then Male-Male vs Female-Female)? I'm looking for an answer to the question of whether one of the types of dyads are more likely to fail to interact. I think the statement "There are more Males-Female dyads that were never seen to interact than would be expected by chance" would be easier to comprehend than if I had to write out the comparisons of the three pairs. I just don't know if it's allowable. Thanks for your time.


title: Correct way to combine bigrams and unigrams in naive bayes classifier?

text: I'm writing naive bayes classifier for spam detection. I'd like it to recognize difference between "Paris Hilton" and "Hilton Paris", so I'm using bigrams in addition to unigrams. Currently I throw both bigrams and unigrams in one bag as independent words, but obviously they're not independent. What is the correct way to combine observed data from both bigrams and unigrams?


title: Controlling p-value in SPSS

text: I have corrected my alpha level using Šidák correction to take into account the experiment-wise error rate, as I am doing multiple correlation tests using the same data set. My new alpha is 0.00056. I tried using SPSS but the p threshold (alpha) seems to be either 0.01 or 0.05 and the p value of the correlation is given only in three decimal places. For example, SPSS returns the p value of one of my correlations as 0.000 and indicated it is significant at 0.01. However, my p threhold is not 0.01 but 0.00056. How do I determine if this correlation is significant? How do I work out the extact p value at 4 decimal points, which will then enable me to compare it with 0.00056. Is there an easy to use but reliable online calculator, or is this adjustment possible in SPSS?


title: Test whether distribution of random values converges to fixed value

text: I conduct a simulation experiment which produces values with some random effects involved. In theory, the results should all be some fixed value $c$. What is an appropriate test to check whether the distribution converges to the fixed value? I saw others using a binomial test for that purpose, but it seem inapprobiate to me here. Maybe a one-sided $t$-test can be used? Thank you for your clarification.


title: Assessing group differences on multiple outcomes

text: I have one independent variable with two levels: Gender (Male and Female), and dependent variables will be composed of three test scores: Impulsivity, Assertiveness, and Social Desirability. My research question is to know the difference between males and females on all three test scores. I'm confused about the kind of analysis I should use, and whether or not I've actually labelled the dependent and independent variables correctly? Should I be using a MANOVA?


title: Is concept of similarity objective?

text: Imagine following example: We have two pairs of points (i.e. 4 objects in some space) and two similarity measures. According to first similarity measure, objects from first pair are more similar then objects from second pair. On the contrary, according to second similarity measure, objects from the second pair are more similar than objects from the first pair. Can this happen (i.e. do similarity measures preserve ordering)? If yes, then does clustering make sense (since we can manipulate with similarity concept by using suitable measure)?


title: Can we predict whether or not persons will favor legalization of marijuana based on sex and three categories of SES?

text: I am trying to figure out the relationship between a dichotomous independent variable sex, an ordinal three level SES independent variable predicting a dichotomous variable for legalization of pot or against. I am confused about what test to run.


title: How to model repeated measures with time-varying covariates in SPSS?

text: I wanted to test the effectiveness of a particular type of "talking" therapy on depression. I envisaged selecting ONE group of people and measuring their heart rate and scores on depression scale for a period of 8 weeks. So I would have 8 measurements for the heart rate (session heart beat average) and 8 measurements for depression scale. How do I go about choosing the right statistical test, and how do I get my head around this in SPSS?!


title: Which is better, SD or SE in line charts?

text: Among the literature, some line charts use mean/SE but others use mean/SD. Which is correct or better?


title: Derivation of M-step for pLSA

text: I was looking at section 6 of these notes and trying to understand the derivation of the M-step at the top of page 10. I understood the derivation for the model without background, but I do not understand where the $(1-P(Z_{d,w}=\theta_B \mid d,w))$ terms come from. Intuitively, I understand that they are down-weighting words that are mostly background and not topic- related. But, how do you get the expressions in the M-step from the derivatives?


title: Why is it difficult to incorporate uncertainty in random effects when making predictions from mixed models?

text: There are several threads on R-sig-ME about obtaining confidence intervals for predictions using `lme4` and `nlme` in R. For example here and here in 2010, including some commentary by Dougals Bates, one of the authors of both packages. I hesitate to quote him verbatim, for fear of them being taken out of context, but anyway, one comment he makes is > "You are combining parameters and random variables in your predictions and > I'm not sure what it would mean to assess the variability of those > predictions. A Bayesian may be able to make sense of it but I can't get my > head around it." https://stat.ethz.ch/pipermail/r-sig-mixed- > models/2010q1/003447.html I know that the Bayesian glmm package `MCMCglmm` can produce credible intervals for predictions. Lately, the development version of `lme4` on github has been given a `predict` method, but it is accompanied by the following comment: > " @note There is no option for computing standard errors of predictions > because it is difficult to define an efficient method that incorporates > uncertainty in the variance parameters; we recommend \code{\link{bootMer}} > for this task." https://github.com/lme4/lme4/blob/master/R/predict.R So, why is it difficult to incorporate uncertainty in random effects when making predictions from mixed models in a frequentist setting ?


title: Method/workflow for analyzing data with changing structure?

text: I am analyzing data relating to networks that evolve over times (more precisely, a snap shot of the network at every discrete time step). Each node of the network denotes a person who perform some action at each time step. As the network evolves, communities may emerge and we run some detection algorithm to identify these communities at each time step. As part of our analysis, we want to focus on these communities and the nodes within each. A difficulty here is that these communities are changing (e.g. the number of communities may be different at each step, so are the memberships of these communities), so there is no way to fix the number of columns/variables needed to account for all possible number of communities (and their various structural properties) emerged at each time step. In step 1, there might be just one community, while in step 6 there might be three. I would like to be able to analyze these data on two levels:   1. For the detected communities, analyze their network structural properties and see whether they stabilize over time (which would imply that the number of detected communities stabilize too).   2. For each node of the network, analyze how their actions (and other characteristics) relate to the communities to which they belong to. If possible, I would like to utilize Stata for such analysis but other tools are OK too. More generally, I am interested in how I should tackle such scenario of analysis properly.


title: How to find the "optimal" cutoff-points in a non-linear relationship?

text: I have the following challenge: The dataset has one dependent and one independent variable which are connected in a non-linear fashion. I am trying to give a more qualitative picture here because I am unsure how to handle the dataset. Basically there are five regions from "one end to the other":   1. "Unstable" relationship with "too few" data points   2. Negative relationship with the dependent variable in negative territory   3. Negative relationship with the dependent variable around "0" (seems also "unstable")   4. Negative relationship with the dependent variable in positive territory   5. "Unstable" relationship with "too few" data points I differentiated regions 2-4 because it makes a difference concerning the actions to be taken for the different regions. **My question**   I need an algorithm (best an implementation/tool) with which I can make the above mentioned regions exact (it is ok to stay with 5 regions). The width of the regions is not equal so the cutoff-points have to be optimized on basis of "distinctiveness". I thought of piecewise regression or optimal width histogram or even some more elusive clustering algorithms but I haven't gotten very far with these approaches (might be my fault). Perhaps you can share some ideas where to start.


title: How to test whether a covariance matrix has changed over two time points?

text: My task is to test if there's change in covariance matrix of 6 variables. Values of 6 variables are measured twice from same subjects (3 years between measurements). How can I do that? I've been doing most of my work using SAS.


title: Help with finding proportion

text: serve speed of a particular player was 102 miles per hour (mph)and the standard deviation of the serve speeds was 13 mph. Assume that the statistician also gave us the information that the distribution of the serve speeds was mound-shaped and symmetric. What proportion of the player's serves was between 89 mph and 128 mph? This is what I have so far: z=(89-102)/13 = -1 z = (128-102)/13 = 2 So, P(-1 < z < 2)? Am I on the right track? If so how do I solve P(-1 < z < 2)? Any help would be appreciated. Thank you


title: Sampling distributions help with questions

text: > **Possible Duplicate:**   >  Probability of mean of random sample being in a certain range When a pizza restaurant’s delivery process is operating effectively, pizzas are delivered in an average of 45 minutes with a standard deviation of 6 minutes. Suppose each night, 9 deliveries are selected at random and their delivery times are recorded. For the sake of argument, assume that the pizza deliveries are normally distributed. (3 pts each)   1. What is the probability that the average time for the sample mean each night is within 2 minutes of the population mean?   2. What is the probability that the average time for the sample mean each night is within 2.5 minutes of the mean?   3. What is the probability that the average time for the sample mean is within 4 minutes of the mean? I'm doing some research and teaching myself how to solve these, any help on the answers?


title: Estimating random effects and applying user defined correlation/covariance structure with R lme4 or nlme package

text: I have following type of data. I have evaluated 10 individuals each repeated 10 times. I have 10x10 relation matrix (relationship between all combination of the individuals).               set.seed(1234)     mydata <- data.frame (gen = factor(rep(1:10, each = 10)),                           repl = factor(rep(1:10, 10)),                           yld = rnorm(10, 5, 0.5))      This gen is different varieties of plant, so each can be repeatedly grown and yield is measured. The covariance matrix is relatedness measure by genetic similarity calculated by ibd probabilities in seperate experiments.               library(lme4)     covmat <- round(nearPD(matrix(runif(100, 0, 0.2), nrow = 10))$mat, 2)     diag(covmat) <- diag(covmat)/10+1     rownames(covmat) <- colnames(covmat) <- levels(mydata$gen)          > covmat                        10 x 10 Matrix of class "dgeMatrix"                               1    2    3    4    5    6    7    8    9   10     1  1.00 0.08 0.06 0.03 0.09 0.09 0.10 0.08 0.07 0.10     2  0.08 1.00 0.08 0.09 0.04 0.12 0.08 0.08 0.11 0.09     3  0.06 0.08 1.00 0.10 0.05 0.09 0.09 0.07 0.04 0.13     4  0.03 0.09 0.10 1.00 0.02 0.11 0.09 0.06 0.04 0.12     5  0.09 0.04 0.05 0.02 1.00 0.06 0.07 0.05 0.02 0.08     6  0.09 0.12 0.09 0.11 0.06 1.00 0.12 0.08 0.07 0.14     7  0.10 0.08 0.09 0.09 0.07 0.12 1.00 0.08 0.03 0.15     8  0.08 0.08 0.07 0.06 0.05 0.08 0.08 1.00 0.06 0.09     9  0.07 0.11 0.04 0.04 0.02 0.07 0.03 0.06 1.00 0.03     10 0.10 0.09 0.13 0.12 0.08 0.14 0.15 0.09 0.03 1.00      My model is:               yld = gen + repl + error       both gen and repl are considered random and I want to get the random effect estimates associated with each gen, however I need to consider the relationship matrix. If it is too complex to fit nested models, I would just remove repl from the model, but ideally I will keep it.               yld = gen +  error       How can I achieve this using R packages, perhaps with nlme or lme4? I know that ASREML can do it but I do not have hold and I love R for being robust as well as free.


title: What type of regression analysis is best to model the relationship between self-efficacy and activity level over three time points?

text: I am conducting an orthopedic study to predict improvements in activity level (dependent variable) based on a type of self-efficacy scale (the independent variable). There are, however, other predictors that I wish to include in the model. All scales are continuous. Since this study will be conducted longitudinally - in other words, I wish to collect this scale data from patients at three different time points (baseline, 3-months, and 6-months) - I am wondering what type of regression analysis is best.


title: Survival Analysis and Predictor Variables

text: Suppose that we are modeling the log-hazard of death using age, gender, bmi, and weight. BMI and weight are correlated. Would it make sense to drop one of these variables from out analysis?


title: How to sample uniformly from an intersection of simplices?

text: The title says it all. Specifically: given strictly positive real numbers $a_1,\dots,a_T$ and $b_1,\dots,b_T$, I want to sample from $$\mu:=\text{Uniform}(\{p\in[0,\infty)^T : \sum_{t=1}^T a_t p_t = 1\}\cap\{p\in[0,\infty)^T : \sum_{i=1}^T b_t p_t = 1\}),$$ assuming that the intersection is nonempty. My best attempt so far is the following rejection sampling algorithm:   1. Set $u_t = \min\{\frac{1}{a_t},\frac{1}{b_t}\}$ for $t=1,\dots,T-1$.   2. Sample $p_t \sim \text{Uniform}(0,u_t)$ for $t=1,\dots,T-1$.      * (Justification: if $p_t>u_t$ then $\max\{a_t,b_t\} p_t > \max\{a_t,b_t\}\min\{\frac{1}{a_t},\frac{1}{b_t}\} = 1$, so $p_t$ cannot be the $t$th coordinate of a sample from $\mu$.)   3. Set $s_a = \sum_{t=1}^{T-1} a_t p_t$ and $s_b = \sum_{t=1}^{T-1} b_t p_t$.   4. If $s_a>1$ or $s_b>1$, reject.   5. If $\frac{1}{a_T}(1-s_a) \neq \frac{1}{b_T}(1-s_b)$, reject.   6. Else, set $p_T = \frac{1}{a_T}(1-s_a)$ and accept $p=(p_1,\dots,p_T)$. The trouble is that in step 5, we will reject with probability 1. Does anyone know how to do this properly, either using something like the above or via a totally different approach? Thanks in advance for any ideas!


title: Assuming two Gaussian distributions of equal mean and variance, then how different can we expect the top X members of each group to be?

text: Here's the thread I got the idea from: http://www.quora.com/Do-men-have-a- wider-variance-of-intelligence-than-women/answer/Ed-Yong Basically, this is a model that might be able to explain why there aren't more females in prestigious math/science competitions - it might be a statistical artifact arising from the simple fact that there are far more males than females in math/science. If this model applies, then we may not need to assume that male intelligence has higher variance than female intelligence. The question I'd like to see addressed: If we **assume** equal means and equal variances (but different sample sizes), then is the model in the paper still the best model when used for predicting, say, the gender composition of the team of the 5-10 best players? Rather than just the gender composition of the grandmaster? http://rspb.royalsocietypublishing.org/content/276/1659/1161.full#sec-3 has the diagram and use of the model They basically used pairing between the top 100 males and top 100 females. Is that a valid assumption to make though? It works for grandmasters - that's true - but would it work if we're trying to select the top 10 people in any field? It's entirely possible, after all, that the expected distributions would be different if we're trying to select from a random distribution of the top 5 players of each gender, rather than the n-th ranked player of each gender. As you increase the number of players you select for a "winning" team, for example, then maybe the distributions play out in a different way. I would expect the smaller group to have higher variance in mean than the larger group. We know that to be true when averaging over the entire population distribution (as a consequence of the central limit theorem). But what if we just want 10 people from each population instead? The fact is that a lot of "potentially" top people will end up dropping out because they would do something other than spend hours a day to practice for a "winning team" High variability of the extreme value though - that makes sense if we're talking about the very top. In a large population, the extreme value is going to be very consistent. Whereas in a small population, the extreme value is going to have A LOT of variability - but that extreme value spends _far_ more time in the left part of the (mean of extreme values) as compared to the right part of it. So if you had a head-to-head match up most years, the population with the larger sample size will win. The thing is, what about a head-to-head matchup of the top 10 members of each distribution? It would be some sort of average between the model the paper used (1 to 1 matchups) and the model where we simply had matchups of the two entire populations with each other.


title: Addressing model uncertainty

text: I was wondering how the Bayesians in the CrossValidated community view the problem of **model uncertainty** and how they prefer to deal with it? I will try to pose my question in two parts:   1. How important (in your experience / opinion) is dealing with model uncertainty? I haven't found any papers dealing with this issue in the machine learning community, so I'm just wondering why.   2. What are the common approaches for handling model uncertainty (bonus points if you provide references)? I've heard of Bayesian model averaging, though I am not familiar with the specific techniques / limitations of this approach. What are some others and why do you prefer one over another?


title: Understanding variance estimators

text: I am having trouble understanding the following. Let $\mu$ and $\sigma^2$ be the true mean and variance, $\bar{x}$ and $s^2$ the measured mean and variance for a random variable $x$, where $$\displaystyle s^2 = \frac{1}{N+k}\sum_i (x_i-\bar{x})^2.$$   * If $s^2$ is an unbiased estimate of variance then $k=-1$.    * If $s^2$ has the smallest mean square spread from true variance then $k=1$.


title: Can I have an example of pairwise comparisons of data against a within-subject factor?

text: Here is an ANOVA model with one between-subject factor (condition; 4 levels) and one within-subject factor (trial_seq; 20 levels).               amod = aov(decision_quality ~ condition*trial_seq +            Error(user_id/trial_seq), data = d.task1)     summary(amod)      I want to do a pairwise analysis for both the between-subject factor and the within-subject factor. I found that it is not straightforward to use a function like TukeyHSD() for a within-subject factor. I researched this problem and found some suggestions in another thread. However, my lack of statistical background prevented me from fully understanding vignettes coming with the multcomp package. Anyway, I tried some random statements like:               summary(glht(amod, linfct = mcp(trial_seq = "Tukey")))      However, it generates the following error:               Error in model.matrix.aovlist(model) :        ‘glht’ does not support objects of class ‘aovlist’     Error in summary(glht(amod, linfct = mcp(trial_seq = "Tukey"))) :        error in evaluating the argument 'object' in selecting a method for          function 'summary': Error in factor_contrasts(model) :        no ‘model.matrix’ method for ‘model’ found!        * What did I do wrong?


title: Assigning class to the cases after K means cluster analysis (SPSS)

text: I have carried out PCA and then clustered the 6 resultant components using K-means clustering technique using SPSS. Normally SPSS adds a class variable for each case indicating its assigned group. Is there any other method I can "calculate" the class variable (i.e. using component scores and REGR factor scores for the K means analysis given in the "Final Cluster centers" table????


title: How to determine the confidence interval or significance of a covariance estimate

text: I was wondering if there is a way to determine the significance of a covariance? So, if have two vectors of returns, and then calculate the covariance, how do I determine if the sample covariance is significantly different from zero? A professor at the university told me to look for the confidence interval? Can use the Wald test for this?


title: Finding optimal combination of parameters for clustering

text: I have a spreadsheet with one object per line. Each column contains values that are parameters of my objects (let's say length, width, height, weight, color). I can classify objects based on color and define classes. Let's say I am blind, and I therefore don't have access to the color parameter. I want to find a combination of the other parameters that can categorize my objects by color with decent accuracy. Is there a way to find this combination of parameters, or do I have to brute force it and test every combination?


title: How to set parameters for SMOTE in WEKA?

text: How I can choose the most optimal parameters for SMOTE in WEKA. More specifically, how should I pick the nearest-neighbours and the SMOTE to apply? I have been using the amount of SMOTE that gives the highest geometric mean and F1 measure but I'm not sure if this is the correct approach.


title: Figuring out quantiles in quantile regression

text: Suppose I have a dataset $\{y_i,x_i\}$ $i=1,2,...n$. For the response variable, $y_i$ as per quantile regression I have the following likelihood: $$p(y_i|\beta,\alpha_i,\sigma) =\frac{\alpha_i(1-\alpha_i)}{\sigma}\exp\left(-\frac{\epsilon_i(\alpha_i-I(\epsilon_i<0))}{\sigma}\right) $$ where, $\epsilon_i=y_i-\beta^Tx_i$, and $\alpha_i$ is the percentile. This is an asymmetric laplacian. Usually for a synthetic dataset you can generate the data and know what the $\alpha_i$'s are (because you generated them).   1. Question is for a real dataset how would you know what the $\alpha_i$ is. My first thought was to put a uniform prior on $\alpha_i$ and look at the posterior, however this would symmetricise the likelihood (atleast intuitively it does). Is ML an option? This would highly depend on $\beta$ which in turn would depend on $\alpha_i$'s.   2. If above estimation of $\alpha_i$ is not possible what datasets would give you the quantile directly (with corresponding $\alpha_i$).


I want to perform test to a 10-fold cross validated result like the following to prove my proposed model is better than the baseline.               Baseline,Proposed     0.807,0.901 - won     0.851,0.861 - won     0.840,0.799 - lose     0.817,0.951 - won     0.907,0.951 - won     ...     0.827,0.911 - won      So my question is, _"Is it right to perform t-test to compare the two cross validated results?"_ The distribution of accuracy is not follow to normal distribution because the range is in [0,1] (finite). So it not seems to be right to use t-test because it assumes that the distribution follows to normal distribution. What is the right method to compare the two accuracy from cross validation?


title: 1-day VaR scaling for higher horizons

text: I have been reading about VaR scaling in "Market Risk Analysis" by Carol Alexander. It talks about square root of time (horizon) scaling is valid only if returns distribution is normal. For stable distributions there is a method, Power Law of exponents, to empirically find out exponent of horizon instead of Square root. I didn't find explanation if returns are not Stable distributions. **What exponent to use if returns are not Stable distribution?** Thanks


title: DLM Trigonometric seasonality vs trend

text: I am working with a structural time series model in R using dlm that has both a stochastic (slope and intercept) local trend and two trigonometric seasonal components (sin and cos pairs). I have difficulty balancing the two. The main frequencies in my problem arise from several diurnal and semi-diurnal tides. I am using recursive forms of the trig components (e.g. from Harvey) and MLE in the package for estimation. All the innovations are assumed uncorrelated and I assume the cos/sin components for each frequency have the same variance. The MLE estimate is not robust to starting guess, which I think has to do with the ambiguity between trend and seasonality I am seeing in my results. Just from looking at the data, there is clear growth and decay of the main frequencies. The model doesn't catch this though -- instead, the periodic component stays stationary and an oscillating trend picks up the changes in amplitude as well as some fluctuations at higher frequncies (double the main trig components). I would like to force the model to evolve the trig components more. I believe this would happen if the trend were more band limited and stayed out of the way. My biggest priority would be representing frequencies lower than the trigonometric ones because these have more predictive value. Is there an appropriate modification of the trend that would do this? Some kind of differencing? Is there a principled way to constrain the innovation variances so that this happens? Thanks.


title: How to find the optimum ordering of a set for a nonlinear equation in decent time

text: Assuming I have an unordered set of tuples of the form: $(a, b)$ How can I order the set so as to optimize the following equation: $\sum_{i=1}^n\sum_{j=1}^{i-1}(a_i*b_i)(1-b_j)$ For small numbers of $n$ I can simply calculate the result for each ordering. However, as this is $O(n!)$, it's not feasible for anything except very small values of $n$. Is there either a closed form solution or a programmatic solution that runs in linear time?


title: Difference between shifted distribution and zero-truncated distribution

text: How would one derive the density or pmf of a distribution that is shifted to the right ? For example, a Poisson distribution that is truncated at zero , is shifted to the right . And it only takes values from $x = 1$, to $x = \infty$. How can I derive its probability mass function ? Do I derive this from the poisson distribution ? I would like to see how it is done mathematically.


title: Implementation of prototype selection in application to nearest neighbor classifier

text: I am looking for an implementation of prototype selection algorithm and it is highly preferable that implementation would be in Python. I am using sikit-learn, but there is none.


title: Using GLM on a continuous response variable

text: Stats noob here who's trying to learn. Let's say that I am developing a glm on a continuous response variable. I've read a number of tutorials on glm and the estimation that it utilizes. However, I'm a little lost on the specifications that are required for developing a glm in R. According to the glm help page, family is supposed to specify the residual distribution. However, I find that my error distribution never matches any of the provided "options" (gamma, gaussian, etc). For example, how would I specify that family if my error ditribution looked as follows. Feel free to provide me with any other information that I seem to be misunderstanding. Thanks! ![enter image description here](http://i.stack.imgur.com/ueRyC.png)


title: Artificial neural network analysis

text: I've designed and programmed an ANN, So now I need to do an analysis about results, I've seen some examples where they do graphics about surface error, average of true outs, and another interesting graphics so I've decided do it, I think to use Gnuplot to do this analysis, but I need some advices or things that i do with the data before plot these. I don't know if I could plot the data as I got or I've to preprocessing before. I've seen this page where are some scripts of data representation (example). More specifically I need to know what kind of analysis I have to do for get information about the efficiency of my network.


title: Is a MANOVA justified?

text: I have to analyze the data from a study in which two groups of subjects (total sample = 30) were tested on 4 cognitive tasks, all realtive to executive functions. For each task, I have accuracy score and reaction times. So, the analysis will be made on 8 dependent variables. In such situation, **is it a good practice to conduct a MANOVA with all the DVs** (or, in alternative, two separate MANOVAs, one for accuracy scores and one for reaction times), **or is it better to analyze the data by a number of t tests**? Mainly, my doubts arise from the fact that I think I don't feel understand when (from a theoretical point of view) a series of variables can be considered as a whole group. Thanks.


title: Are these populations non-random and different?

text: What would be the simplest, straightest-forwardest way to determine the following:   1. Whether A,B,C is non-randomly distributed in each single group (IE: are there more A's in group 1 than random chance would predict?)   2. Whether the distribution of A,B,C is different between the groups. (IE: Group 1 is has more Cs than Bs, so does Group 2 - but is the difference between the two groups significant?) Data:               Group,  A,   B,   C,        1,   126, 357, 348   N=382        2,   86,  196, 139   N=207        3,   63,  185, 162   N=193        4...        5...      And in handy R-ready format:               A <- c(126,86,63,54,47,40,32,32,29,29,27,26,20,18,14)     B <- c(357,196,185,137,95,74,45,69,64,49,54,80,62,41,56)     C <- c(348,139,162,126,82,69,35,63,40,42,40,55,44,29,35)     N <- c(382,207,193,143,100,80,45,70,70,53,55,84,67,42,57)      A,B,C are counts of presence/absence of traits within each group, so lots of As in group 1 are not . Groups are each separate non-related populations. Example: lets say each group represents a species of lizard, while A, B, C indicate whether the lizard has spots on its head(A), Body(B), or Tail(C). For species 1 (Group 1) 382 lizards were examined, 126 had spots on their head, 357 on their body, 348 on their tail... Of 207 lizards in species 2, 86 had spotty heads, 196 spotty bodies, 139 spotty tails. So, is spottiness non-random for members in a group? And does spottiness vary significantly between the groups? I think this is a basic fundamental question, but while I have been bashing through countless pages of theory explaining different models looking at more complex situations, I have totally lost sight of the basics -- which I have only just recently learned. Thanks for helping me get back on track.


title: Collapsing regression of nominal data to single significance value

text: I'm working with a dataset of combined continuous and categorical data, and I'm constructing a regression model based on both types of data simultaneously, using dummy variables. The model will look something like this: $y = \alpha + \beta_1 x_1 + \beta_2 x_2$ where, for example, $x_1 \in [0,1]$ is continuous, and $x_2 \in \textrm{{a, b, c}}$ is nominal. Because of the dummy variables, I will obtain different significance values for each of the three nominal categories. My question is, **is there a simple, commonly used method for collapsing the three significances into a single significance value, taking into account both the p-value and coefficient value?** I say both p-value and coefficient value because if one term has a very high significance but a very $\beta$, I imagine that should be factored in. From some discussion with others in my group, we determined that this method (if extant) would have to be some non-linear transformation of the $\beta$ and p-value. The simplest method we came up with was simply select the minimum p-values with an optional weighting term to account for $\beta$. Are there any accepted ways of accomplishing this?


title: How to compare two series of gene expression levels?

text: I have a data set with gene expression levels measured for three different treatments, and for each treatment we have three biological replicates. Treatments are: Control, 24h at $4$°C, and 24h at $-3$°C. This analysis was carried out on two genotypes.                                      CTRL1 CTRL2 CTRL3 4°C_1 4°C_2 4°C_3 -3°C_1 -3°C_2 -3°C_3     Gene1  Genotype 1       20    22    25     85    92    120   450    380    372            Genotype 2       50    72    78     130   122   184   250    490    700      I'd like to know how to compare `gene1` expression level in `genotype1` vs. `genotype2`. The aim is to find a similarity between expression trend in `genotype1` and `genotype2`.


title: How to calculate Cohen's d effect size intervals for a within-subject design?

text: I have looked at different web tools for calculating confidence intervals for effect size like Cohen's d for within-subject design. My experience from this survey tells me that there are different ways to estimate within SD because results I got differs on both effect size and confidence intervals. They usually use correlations between variable to estimate pooled variance. Can anyone point me to adequate way of calculating these things? In addition I mainly use effect size confidence intervals to determine when intervals are non-overlapping the zero effect as hypothesis testing. Does this produce the same results as using a standard confidence interval for the difference directly? Thanks!