title: Why z=f(x) does not imply E[z]=f(E[x]) when f is not linear?

text: Why $z=f(x)$ does not imply $E[z]=f(E[x])$ when f is not linear? I can give an example, but I couldn't derive the general form. Let $z = {x^2}$ Let $g(x)$ be the pdf of $x$. Then: $ E\left[ z \right] = E\left[ {{x^2}} \right] = \int\limits_{ - \infty }^\infty {\left( {{x^2}} \right)g\left( x \right)dx} $ $f\left( {E\left[ x \right]} \right) = {\left( {\int\limits_{ - \infty }^\infty {\left( x \right)g\left( x \right)dx} } \right)^2}$ $ \Rightarrow E\left[ z \right] \ne f\left( {E\left[ x \right]} \right)$


title: Negative R squared contradicts ssa/sst?

text: I understand from this question - When is R squared negative? that the R squared value of a linear regression model can be negative if the intercept is constrained. And this makes sense if you define R squared as - $$R^2 = 1-\frac{SSE}{SST}$$ One says $SSE>SST$. But then, $$SST = SSA + SSE $$ Total Sum of squares = Sum of squared errors + Sum of squared residual. And with this we get - $$R^2 = \frac{SSA}{SST}$$ And now it is hard to imagine how $R^2$ can be negative. Aren't SSA and SST >0 always?


title: Multivariate Kolmogorov-Smirnov for two samlpes

text: Let’s say we have data of some fruits and we want to consider different statistics for each fruit (there are various samples for each statistic). For example: Watermelon:   * Sugar: 50, 57, 36…    * Water: 101, 143, 128…   * … Melon:   * Sugar: …   * Water: …   * … (…More fruits…) We want to know if, based on this statistics, fruits are the same. I thought about using Kolmogorov-Smirnov test in order to check it. So we have a CDF for each statistic (sugar, water, etc.) for each fruit. If we want to consider all statistics for each fruit, it would be nice to use the multivariate Kolmogorov-Smirnov test for two samples and compare every pair of fruits. The problem is that I cannot find information about the multivariate test with two samples and I read somewhere that the multivariate test is not that good. Could you tell me a different test I could use? Or the Kolmogorov-Smirnov is a good idea (where can I read about the multivariate for two-samples)?


title: Suitable regression model for limited, discrete dependent variable

text: I have data where the dependent variable is discrete and lies between 20 and 40 (possible values are 20, 20.5, 21, 21.5, ..., 39, 39.5, 40). The variable measures some results from a game which can be between 20 (lowest achievable value) and 40 (highest). After some hours of research on the web, I could not find a regression model that ideally fits those characteristics of the described dependend variable. Given what I found, maybe a multinomial logistic regression fits best although my dependent variable is not nomial. I would be very thankful if you could propose me a regression model that fits my problem best according to your consideration. Thanks a lot!!! :)


title: Logistic regression with an independent variable that only applies to a subset

text: I'm trying to perform a logistic regression with the dependent variable being whether a potential customer who attended an event purchased or did not purchase. A subset of potential customers at the event participated in a meet-and-greet type activity with each other and existing customers. For each potential customer, the percentage of existing customers was recorded for the time a potential customer was in the meet-and-greet area. For example, a potential customer may have stopped at the meet-and greet for 1 hour, and during that time, on average there might have been 20% existing customers, and 80% potential customers. I would like to incorporate this variable into my logistic regression model as a continuous variable without categorizing it, but I do not know how to code that variable for potential customers who were at the event, but not the meet- and-greet. I think that setting the non-meet-and-greet people to zero would not be the right way to go. Do I just need to perform a separate logistic regression on the people who participated in the meet-and-greet? I'm using SPSS if that makes a difference. Thank you!


title: Normalization for pattern classification?

text: I'm working off my first independent project for some pattern classification. I'm utilizing some datasets from UCI machine learning, but am not sure on how to start with data normalization. The data isn't that large (feature vector around 15-20 dimensions), but I'm thinking there still needs to be some type of normalization done in order for my classifiers down the line (SVM + KNN) performs properly I guess my main question is regarding when to normalize. I currently don't have my data split between testing and training data. Thus, should I normalize all data now? Thanks ahead of time


title: How do you explain dimension-reduction with statistics?

text: With statistics: how would you explain the dimensionality reduction and dimensionality addition? Like the conversion of a color picture to gray space so that a color blind person could more easily understand what was going on or a painting to a sculpture so that a blind person might enjoy it? In the first instance i call it a conversion to a subspace containing less information a vector transformation to a lower dimension. In the second instance i call it a transformation and extrapolation into a superspace a great dimension with possibly more information depending on the artist involved. Also how do you explain the transcription of vocal frequency information and recordings into a system of musical representation which is not as accurate like western scales?


title: Convergence in distribution results not deriving from central limit theorem?

text: In the stats class I took, all the results I have encountered about the convergence in distribution of some random variables are in one way or another consequences of the Central Limit Theorem. Out of pure curiosity, I wonder whether there exists useful and commonly used results about the convergence in distribution of some statistic $f(X_1,\dots,X_n)$ (where $(X_1,\dots,X_n)$ could be a random sample)   1. Which do not follow from the Central limit theorem.   2. Which holds for any distribution of the individual $X_i$, or under some "mild" conditions on their distribution.  (sorry for the "mild" I know it means nothing precise, so my appreciation of answers will be partly subjective, but I do not know of a way to formalize this. As an example, the CLT requires that $E(||X_i||^2) < \infty$, that is mild enough for me :).) I found some interesting examples at http://www.math.uah.edu/stat/dist/Convergence.html, but appart from the CLT, they all violate condition 2 above. **Remark :** I realize that "which do not follow from the Central limit theorem" might be a little vague. I do not know how to express this more precisely, but I guess I am curious about results which would not take the form :   * $ X_n $ converges in probability to $X$.   * $Y_n$ converges in distribution to $N(a,V)$ by the CLT.   * So $X_nY_n$ converges in probability to $D \sim X N(a,V)$. (e.g. Wald statistic converges in distribution to $N(0,A)'A^{-1}N(0,A) \sim \chi^2_q$) Hope the question is clear enough, feel free to ask for clarifications.


title: Regularized fit from summarized data

text: I have a multiple linear regression problem $y=X\beta+\epsilon$. The number of observations $m$ is large, so by the time the data gets to me it's been summarized into:   * $m$   * $X^TX$   * $X^Ty$   * $y^Ty$   * $\sum_{i=1}^m{y_i}$ The above list appears sufficient to compute the OLS estimate. But can it be used to compute a regularized estimate of some sort (e.g. ridge regression)? If not, can the list be augmented in the same spirit to enable regularization?


title: Is it ok to use Manhattan distance with Ward's inter-cluster linkage in hierarchical clustering?

text: I am using hierarchical clustering to analyse time series data. My code is implemented using the _Mathematica_ function `DirectAgglomerate[...]`, which generates hierarchical clusters given the following inputs:   * a distance matrix D   * the name of the method used to determine inter-cluster linkage. I have calculated the distance matrix D using Manhattan distance: $d(x,y) = \sum_i|x_i - y_i|$, where $i = 1,\cdots, n$ and $n \approx 150$ is the number of data points in my timeseries. My question is, is it ok to use Ward's inter-cluster linkage with a Manhattan distance matrix? Some sources suggest that Ward's linkage should only be used with Euclidean distance. Note that `DirectAgglomerate[...]` calculates Ward's linkage using the distance matrix only, not the original observations. Unfortunately, I am unsure how _Mathematica_ modifies Ward's original algorithm, which (from my understanding) worked by minimizing the error sum of squares of the observations, calculated with respect to the cluster mean. For example, for a cluster $c$ consisting of a vector of univariate observations, Ward formulated the error sum of squares as: $(\sum_j||c_j - mean(c)||_2)^2$ (Other software tools such as Matlab and R also implement Ward's clustering using just a distance matrix so the question isn't specific to Mathematica). Thankyou very much for your help!


title: Is there an equivalent of ARMA for rank correlation?

text: Hi I am looking at extremely non linear data for which the ARMA/ARIMA models do not work well. Though, I see some autocorrelation, and I suspect to have better results for non linear autocorrelation. 1/ is there an equivalent of the PACF for rank correlation? (in R?) 2/ is there an equivalent of ARMA model for non linear / rank correlation (in R?) fRed


title: How to calculate the number of samples necessary to represent a population distribution?

text: I guess this is similar to this one but not quite the same. Consider I have a distribution that represents a population (e.g., Uniform(0,1)). I separate this interval (0,1) into three equal parts (0-1/3, 1/3-2/3 and 2/3-1). For the population, I will have for each interval 33.333% of data. The graph for the population is presented below.![Uniform distribution divided into three parts](http://i.stack.imgur.com/DzTr6.png) I need samples that represent the populations "well enough". For example, for Uniform(0,1), I could have a sample with 50,000 elements and a graph as presented below. ![Distribution of samples of a Uniform distribution](http://i.stack.imgur.com/hQPsq.png) I could get a huge number of samples (e.g., 100000) if performance was not a significant issue, but it is not the case. I am considering that "well enough" is with 4% margin of error and 95% confidence. So, how to calculate the minimum sample size I need to, with a margin of error of 4% and 95% confidence (for each interval), to represent the population (given that the population size is unknown)? ps.: with "represent the population" I mean having a graph of samples similar to the graph of the population. In other words, for the example and parameters given, with 95% confidence, none of the bars will be below 31.33% or greater than 35.34%


title: Why doesn't it seem to be standard to multiply prior scale by fraction of non-zero predictor observations?

text: When doing multivariate regression, it's often the case that some predictors often have many zero values - dichotomous inputs, dummy coding of polychotomous inputs, interval coding, etc. The fraction of nonzero observations for a predictor is especially decreased when interacting these variables. It seems obvious that only the observations of a predictor for which it is nonzero (after any transformations) will affect the estimation of its coefficient. When doing multivariate regression I like to use weakly informative priors (like Gaussian or Cauchy) for regularization. I've noticed that sometimes, predictors which are usually zero still take on implausible regression coefficients, like ~6 for a logistic regression model. This seems to be because a variable that is nonzero for only a handful of observations is still overcoming the prior because it's getting full credit for all of its zero observations, and thus getting a huge coefficient, which is in many cases a prior obviously overfitting, and hurts generalization performance for prediction. Multiplying the prior's scale by the fraction of observations of a predictor that are nonzero fixes this and greatly improves generalization performance, makes coefficients plausible, and seems to be common sense. So - why don't more people do this? I've tried a cursory search in common textbooks (various machine learning books, applied and theoretical Bayesian statistics), R packages on CRAN, some machine learning packages, and I don't see evidence of anyone else doing this... isn't this just common sense? Or is there something wrong with this approach that I'm missing?


title: Exponential distribution probability

text: Let X and Y be independent random variables such that $X∼Exp(1)$ and $Y∼Exp(2)$. Find the probability that $3X+4Y≤5$. I thought the equation for solving this was the integral from $\int_0^a \int_ 0^{(a-y)} \lambda e^{-\lambda x}*\lambda e^{-\lambda y} dx dy$ where $a = 5$, so I first started with the integral from $\int_0^5 \int_0^{(5-y)} 3*1e^{-1x}*4*2e^{-2y} dx dy$ and got $23.83828927$ which is much larger than a probability can be. Any help?


title: Variance of sample mean in an AR(1) process

text: I learned that for an AR(1) model $x_n = b + \phi x_{n-1} + a_n$ with $|\phi| < 1$ and $a_n \sim WN(0, \sigma_a^2)$, the covariance between $x_n$ and $x_{n+h}$ is $$ \mathrm{Cov} (h) = \frac{\sigma_a^2 \phi^h}{1 - \phi^2}$$ Now I would like compute the variance of the sample mean $$ \mathrm{Var}\frac{\sum_{i=1}^n x_i}{n} = \frac{n \mathrm{Cov}(0) + \sum_{i=1}^{n-1} 2(n-i) \mathrm{Cov}(i)} {n^2} = \frac{\sigma_a^2}{(1 - \phi^2) n^2} [n + \sum_{i=1}^{n-1} 2(n-i) \phi^i] $$ Can it be simplified further? Thanks!


title: Correlation and group difference with unequal sample sizes

text: ![try to do linear regression](http://i.stack.imgur.com/e4SdT.jpg)![enter image description here](http://i.stack.imgur.com/kqEby.jpg)My sample is 60 participants. At first, I did not expect to analyze the correlation, but the result is interesting, so I would like to. The 60 participants are 50 people who have good reading skill AND another 10 people who have poor reading skill. They took a test for which the full score = 7. The good reading skill group got an average score = 6.7. The poor reading skill group got an average score = 5.1. **Questions:** Is there a test to find these relations?   1. Do the good reading skill group and poor reading skill group have significantly different scores?   2. Does reading skill relate to achieving score?   3. Does good reading skill affect score positively? I have SPSS, but do not know how to analyze these relationships. Thank you.


title: Response Surface Data Mining methodology

text: Recently, I used the methodology of response surface data mining proposed by Lee and Shin (2010) for job stress evaluation. in their study they considered about 80 factors affecting job stress and their purpose is to select the significant factors. in my study i considered 30 quality management practices and i studied their effect on knowledge creation in service provider corporation in jordan. first questions: is 30 factors enough for efficient use of data mining? second question: according to lee and shin results, after dimensionality reduction using data mining filtering (weka software), seven factors were considered significant, and after integrating data mining results into response surface methodology, the 7 practices also were considered significant. in my case study among the 30 practices examined 8 were considered significant using data mining procedure, and after using response surface methodology only 5 were considered significant. how can i discuss this result??? is there a problem with my data? or my results contradict their results??


title: Good resources on Aalen additive models (survival)

text: I'm trying to do some modelling of time to event data and for various reasons, the Aalen additive model seems promising. I'm using the `aareg` function in the `survival` package, but I'm struggling to find much literature (that isn't paywalled) to get more of a feel for the procedure. Does anyone have any suggestions of some anything reasonably in depth, but accessible to someone without expert level knowledge of newer survival analysis techniques?


title: Distant supervision: supervised, semi-supervised, or both?

text: "Distant supervision" is a learning scheme in which a classifier is learned given a weakly labeled training set (training data is labeled automatically based on heuristics / rules). I think that both supervised learning, and semi- supervised learning can include such "distant supervision" if their labeled data is heuristically/automatically labeled. However, in this page, "distant supervision" is defined as "semi-supervised learning" (i.e., limited to "semi- supervision"). So my question is, **does "distant supervision" exclusively refer to semi- supervision?** In my opinion it can be applied to both supervised and semi- supervised learning. Please provide any reliable references if any.


title: Advantages and limitations of Gwet’s AC1 statistic and PABAK

text: I am looking for alternatives to the kappa to assess inter-rater agreement. I've come across two hopefulls: Gwet’s AC1 statistic and PABAK. I wondered what the advantages and disadvantages of each were?


title: Biased estimates when intercept is included in a linear regression

text: I am simulating 10000 data-sets, each of length 20, that follow an autoregressive model with lag 1, using the following code:               set.seed(1)     N = 20      n.reps <- 10000     burn.in = 50     total <- N + burn.in     x <- matrix(NA, n.reps, total)     x[ , 1] <- 0.1     noise <- rnorm(n.reps*total, 0, 0.1)          for(j in 1:n.reps)        for(i in 2:total)         x[j, i] <- 0.5*x[j, i-1] + noise[j + i]          x <- x[ , -(1:burn.in) ]      Then I estimated the autoregressive coefficients (0.5 in this case) for each observed path, using two linear regressions:               res1 <- res2 <- numeric(n.reps)           for(i in 1:n.reps)     {     res1[i] <- lm(x[i, 2:N] ~ -1 + I(x[i, 1:(N-1)]))$coef     res2[i] <- lm(x[i, 2:N] ~ I(x[i, 1:(N-1)]))$coef[2]     }          mean(res1) # 0.4687619     mean(res2) # 0.3845817      In the first line I'm fitting the right linear model, while in the second I'm including an intercept that is not present in the data generating process. I expected both methods to give me unbiased estimates of the autoregressive coefficient, but it looks like including the intercept makes the estimates be biased downward. The bias disappears as I increase the sample size N. The second linear model nests the real one, so I expected higher variance in the estimates when the larger model is fitted, but not bias. So my question is: am I doing something wrong? Thanks!


title: How to calculate mutual information?

text: I am a bit confused. Can someone explain to me how to calculate mutual information between two terms based on a term-document matrix with binary term occurrence as weights? $$ \begin{matrix} & 'Why' & 'How' & 'When' & 'Where' \\\ Document1 & 1 & 1 & 1 & 1 \\\ Document2 & 1 & 0 & 1 & 0 \\\ Document3 & 1 & 1 & 1 & 0 \end{matrix} $$ ![Mutual Information](http://i.stack.imgur.com/xMbTr.png) Thank you


title: lift/uplift/netlift modeling: CART analysis maximizing a t-statistic from a multivariate regression

text: For tree based modeling most packages use some data element as the dependent variable. In uplift modeling the goal is to maximize the difference between treatment and control groups. One way to implement this is with a CART which maximizes the t-statistic for an interaction term between the split (1 = left, 0 = right) and the treatment variable (1 = treated, 0 = control) b_3 below, y = b_1 * T + b_2 * P + b_3 * T * P as I understand it, the interaction term (T*P) captures difference in treatment effects between the two sides of a split. being new to CART, I don't know which software packages would enable this analysis. I program in python (currently exploring orange), r and sas, so I would appreciate any pointers to example/documentation on similar projects.


title: Omega H asymptotic

text: In the omega function in the psych package for R, we have several outputs including Hierarchical omega, Asymptotic omega and Total omega. I was reading about the outputs, but I find nothing about what the asymptotic omega means. For an analysis, a friend of mine calculated an omega H = 0.70 and an asymptotic omega of 0.98. I'd rather not ignore this big asymptotic omega. Can somebody explain what the asymptotic omega is and how it is different than the hierarchical omega?


title: Putting stationary variables through Johansen procedure

text: Is it okay to feed $I(0)$ variables into the Johansen procedure? I've read three sources that seem to state that this is not what you're supposed to do. However, whenever I've done this, I notice that $\Pi$ is full rank and so it leads me to a VAR and therefore I don't see any problem with this.


title: Time-series autocorrelations all positive

text: I've got 36 months of timeseries data, and eyeballing it, it has a linear trend upward. I wanted to do a little more than just eyeball it though. So I put together a correlogram of autocorrelation values for lags 1 to 12 on my observations. They're all positive. Would you say that's good evidence that I've got a linear trend? See a picture of my autocorrelations here. And would that spike back up at the end there likely be due to some seasonality (that'd make total sense for my business's data...we do heavy sales at Xmas and Valentines day)?


title: How can I get cut-off point in multivariated ROC analysis

text: If I have 1 independent variable (continues) and 1 dependent variable (binary), I can conduct logistic regression and ROC analysis, and I can get a cut-off point of independent variable using ROC analysis. But, If I have 2 independent variable, how can I get two cut-off point in multivariate ROC analysis. In STATA, I can get just maximal AUC score but I can't get the cut-off points each independent variables.


title: How to handle leverage values?

text: I have a dataset with 1747 observations. Outcome variable is categorical, while independent variables are continuous, so I decided to use logistic regression for my analysis. I built the model using backward elimination algorithm, and resulting model seemed to be quite stable. Checking VIF values provided very good results, less than 5% of residuals were large and DFBeta values are less than 1. However, 76 cases were identified as large leverage values (only one of them as large residual as well). When I remove these cases, results are somewhat different. What would be the best approach to handle these values? Thanks,


title: Update rule for beta distribution with fixed K/confidence/sample size

text: Normally you have a beta distribution with shape parameters $a$ and $b$. The mean of this distribution is $a / (a + b)$ and the sample size, or the confidence (or K) is $a + b$. Now, if you do some trials, with let's say N positive outcomes and M negative outcomes, you end up with a posterior distribution that is $\text{Beta}(a + N,~b + M)$. So, now your mean is $\frac{a + N}{a + N + b + M}$ and the sample size / confidence / K is $a + N + b + M$. Now, my question is: **what if you want to keep the K / confidence level fixed?** So let's say K should always be 10. So $\text{Beta}(5,5)$ is fine, as is $\text{Beta}(9,1)$, or $\text{Beta}(1.23, 8.77)$. In other words: in this case I would like the mean of the posterior to be able change to reflect the evidence found in the new data, but the confidence level should remain the same (rather than increasing all the time). Is there a simple update rule for this scenario as well?


title: How to check if removing a sample makes a difference in mean and stdev values?

text: I'd like to ask if someone could help me with the following problem: we have measured the same sample 5 times and we would like to check if there are significant differences in mean and stdev values if we use:   * All 5 datapoints   * Only the last 4 datapoints   * Only the last 3 datapoints We have performed ANOVA analysis but we are not sure about the results because we might not have homocedasticity. Which tests would you do to investigate this issue? Thanks in advance for your help.


title: Why are the Lagrange multipliers sparse?

text: I've read that for the Maximal Margin Classifier SVM, after solving the dual problem, most of the lagrange multipliers turn out to be zeros. Only the ones corresponding to the support vectors turn out to be positive. Why is that?


title: LSE estimator for the mean of a multivariate Gaussian?

text: For simultaneous estimation of several parameters, the combination of LSE estimators for each parameter is usually not admissible under squared error for the vector of the parameters. For example, when estimating the mean vector $\theta$ of a multivariate Gaussian distribution $ N({\boldsymbol \theta}, \sigma^2 I)$ with known $\sigma^2$, the combination of LSEs for each dimension of $\theta$ is $X$, and James–Stein's estimator dominates $X$ when the dimension of $\theta$ is greater than $2$, under the squared error. Is James–Stein's estimator admissible? If not,   * what is the LSE in multivariate?   * what is an admissible estimator under the squared error for the mean of a multivariate Gaussian distribution? Thanks and regards!


title: Conceptual distinction between heteroscedasticity and non-stationarity

text: I'm having trouble distinguishing between the concepts of scedasticity and stationarity. As I understand them, heteroscedasticity is differing variabilities in sub-populations and non-stationarity is a changing mean/variance over time. If this is a correct (albeit simplistic) understanding, is non-stationarity simply a specific case of heteroscedasticity across time?


title: How to find marginal distribution from joint distribution with multi-variable dependence?

text: One of the problems in my textbook is posed as follows. A two-dimensional stochastic continuous vector has the following density function: $$ f_{X,Y}(x,y)= \begin{cases} 15xy^2 & \text{if 0 < x < 1 and 0 < y < x}\\\ 0 & \text{otherwise}\\\ \end{cases} $$ Show that the marginal density functions $f_X$ and $f_Y$ are: $$ f_{X}(x)= \begin{cases} 5x^4 & \text{if 0 < x < 1}\\\ 0 & \text{otherwise}\\\ \end{cases} $$ $$ f_{Y}(y)= \begin{cases} \scriptsize{\frac{15}{2}}\normalsize y^2(1-y^2) & \text{if 0 < y < 1}\\\ 0 & \text{otherwise}\\\ \end{cases} $$ I understand how the density function $f_X$ is calculated, by integrating $f_{X,Y}$ from $0$ to $x$ with respect to $y$. I'm however totally lost on $f_Y$, where is the $(1-y^2)$ coming from? If I integrate from $0$ to $1$ with respect to $x$ then I only get $\scriptsize{\frac{15}{2}}\normalsize y^2$, and why is the range $0 < y < 1$? I've graphed the support for $X,Y$, all values where $f_{X,Y}>0$ are colored blue: ![The support for $X,Y$](http://i.stack.imgur.com/CG8sX.png)


title: LibSVM weights of support vectors

text: I am using `LibSVM` classifier in my Java code and I am getting correct results as I verified that with weka GUI, however, when I want to get the weights of the support vectors, I get an empty string. This is a snapshot of my code:               LibSVM svm = new LibSVM();     svm.buildClassifier(train_data);     Evaluation evaluation = null;     evaluation = new Evaluation (train_data);     evaluation.crossValidateModel(svm, train_data, 10, new Random(1));     System.out.println("weights = "+ svm.getWeights());      So, I am wondering if there is any missing step that I should do in using the `LibSVM` so that I get the weights? Thanks for your help.


title: Over-estimated response

text: I have a large survey data set, about 15 years of data. One binary variable (yes/no) in which I'm interested was over-estimated due to a data collection error for about 1 1/2 years. I'm not sure how to handle these data. Could it treated as missing data? Seems like MNAR in that case. If so, and all the other missing data is MAR, how do you run an imputation model in SAS or IVEware? Suppose one could assume the missing values were simply a linear between the two blocks of non-missing data.


title: Difference in two Normal Distributions?

text: I am unsure what I should do if I have two Normally Distributed variables with known parameters and I want to find the probability that one of these variables is greater than the other. Should I use the distribution of differences of the two normal variables, as is shown here?


title: Significance and effect size in treatment design when treatment-response in measured

text: I have a number of plants and all of them undergo a 'treatment'. Then I measure the amount of plants that react to the treatment and the amount that do not. Is there a way to calculate a test statistic (to test for the significant effect of treatment)? If yes, is there a way to infer an effect size measure? E.g.: R code:               n_total <- 120     n_response <- 100     n_no_response <- 20      Or would it be appropriate to use McNemar test in this case in the following way and the use effect size measure as described here:               tab <- cbind(c(n_total, 0), c(n_response, n_no_response))     mcnemar.test(tab)


title: Statistical modelling of outcome of an inquiry

text: Suppose a population of dentists are sent an inquiry. A number of inquiries are filled in and sent back. Each question on the inquiry represents a certain number of points. The total number of points represents the precentage that a dentist adheres to some safety regulation. I would like to test whether there is some significant relation between the outcome of the question "how long are you employed as a dentist" versus the total number of points. Ie., I would like to test whether older/younger dentists work more or less by the rules. I am not quite sure how I would go about testing this. I thought about checking whether the percentage over the whole sample is normally distributed, and then testing whether the percentage within an agegroup is likely to come from the same distribution, but I'm not sure if this is the best way to go about it. Any help would be greatly appreciated.


title: Generating correlated numbers with predefined correlations

text: I'm trying to generate correlated data (preferably multinormal) with predefined correlations (e.g. 0.35 or 0.9). Any idea how I can do it? I'm using R and I did find a way to generate this (using mvrnorm), but you need to supply a covariance matrix. I have a covariance matrix with correlations around 0.9; however, I don't know how I can modify its entries to change the correlation. If I can do that, I'll be able to generate correlated data with the correlations I need. Regards,


title: Multiple comparisons process using Kruskal-Wallis, G-Test and corrections

text: I want to know if changing a button color has a positive effect on different measures (see the 4 measures below), using the following process:   * There's a default button, which is the fixed control group with no changes. This means I'm interested in one-tailed p-values for pair-wise comparisons.   * 3 different button colors, which are the 3 treatment variations.   * I'm using 4 measures to see if the 3 treatment buttons performed better than the control button: 1) Percentage click on the button per treatment 2) Percentage of users signing up per treatment 3) Percentage of paying users per treatment 4) Viral Coefficient per treatment (measures how many new users 1 user invites back to the website)   * Thus a total of 12 pair-wise tests will be performed. **Steps:**   1. For measure 1, I'm doing a one-way ANOVA G-Test (using 4 variables: control + 3 treatments) and multiply that p-value by 4 (using Bonferroni, as I'm doing 4 measurements), then check if there's a significant difference. If there is, I'm doing 3 pair-wise post-hoc G-tests (control <-> treatment 1, control <-> treatment2, etc) and divide that p-value by 2 (to make it one-tailed), then multiply that by 12 (Bonferroni) to find which has a significant effect.   2. For measures 2 and 3, the same procedure as in step 1 is used.   3. For measure 4, I'm doing a one-way ANOVA Kruskal-Wallis (using 4 variables: control + 3 treatments) and multiply that p-value by 4 (using Bonferroni again, as I'm doing 4 measurements), then check if there's a significant difference. If there is, I'm doing 3 pair-wise post-hoc Mann-Whitney tests (again control <-> treatment 1, control <-> treatment2, etc) and divide that p-value by 2 (to make it one-tailed), then multiply that by 12 (Bonferroni) to find which has a significant effect. **Questions:**   * If the process I'm using is the right one for the problem:     1. Is the Bonferroni correction of multiplying the p-value by 4 for each one-way ANOVA correct? My rationale for this is that because there are 4 one-way ANOVA's, you need 4 times the proof before deciding one is signifianct.     2. As for the Bonferroni correction for each post-hoc pair-wise test, is it correct to multiply the resulting p-value by 12?     3. Is it correct to divide the p-value from the post-hoc pair-wise G-test by 2? My rationale for this is that the G-Test returns a two-tailed p-value, but I'm interested in one tail only: is treatment performing better than the control?. Same question goes for diving the pair-wise Mann-Whitney two-tailed p-value by 2.     4. I know the Bonferroni correction is a very conservative approach. So multiplying all post-hoc pair-wise p values by 12 may be too much. The problem is, part of the samples or not independent, but I have no idea how much. This can change for each future experiments. I feel Bonferroni in this case is still a good adjustment, or is there something less conservative? Sidak won't work in this case, as the samples are not 100% independent.   * If what I'm doing is a bad idea in general, would it be better to just run new individual experiments for each measurement separately? This will take 4x as much time for same sample sizes, so if prefer something like the above.


title: Regression with missing covariate data

text: My question is related to regression with missing covariate data. I have a sample of count data representing number of explosion incidents at each year. First $M$ observations were made by considering one criteria and the last $K$ were made considering two criteria. E.g. for the first $M$ years explosions were recorded if there were some human injuries and then, due to the change of regulations, explosion were recorded if they caused human injuries OR the damage were higher than, say, $D$. I introduce two binary covariates: $X_{1}$ for human injury criteria and $X_{2}$ for damage criteria. So, I get that the first covariate $X_{1}$ is always observed and equal to 1 for the first $M$ years but not for the next $K$ years. The second covariate $X_{2}$ for the first year is not observed at all (so it might be 0 or 1), while for the last $K$ years both covariates are not observed but I know that at least one of them is equal to 1.   At first I though that I just put a Bernoulli distribution on these two covariates and integrate out them. Then I perform Bayesian analysis. So, at first I state the model (for simplicity, assume we have just one covariate) $$Y|\theta_{1},\theta_{2}\sim \mathrm{Poisson} \left (\mathrm{exp}\left [ \theta_{1}X_{1}+\theta_{2} \right ] \right )$$ $$X|p\sim \mathrm{Bernoulli}(p)$$ Then perform integration over $X_{1}$: $$\pi(Y|p,\theta)=\frac{1}{Y!}\left ( \mathrm{exp}\left ( -e^{\theta_{1}+\theta_{2}}+Y\left ( \theta_{1}+\theta_{2} \right ) \right )p + \mathrm{exp}\left ( -e^{\theta_{2}}+Y \theta_{2} \right )(1-p) \right )$$ Having this, I form likelihood, state priors and obtain estimates. As I understand, such approach is a so called "Missing completely at random or MCAR". However I understand that in my problem the missingness model is different from MCAR, but I do not know what model would be more appropriate. Any advise would be appreciated.


title: Given a two populations, each split by "categories", how do I test the hypothesis that the mean is larger in the second population?

text: The means for each category varies, but my hypothesis is that the second population has a higher mean within each category. The variance of each category is unknown. The example I'm working with is a set of schools, all of which have been given new computers. The question is does, the introduction of computers improve the students' grades, and by how much? EDIT: The data is average grade by school, each week.


title: Statistics for rounded data

text: I wonder about the existence of "standard" statistical procedures for rounded (log)normal data. Indeed, in my work I often encounter rounded data which potentially cause some problems: "awful" qqplots, zeros for lognormal data... Do there exist in the literature some theoretical studies of usual statistical procedures (such as goodness-of-fit, t-test, ...) for rounded data ? EDIT: I have just discovered the MCMCtobit() function in the MCMCpack R package: _Gaussian Linear Regression with a Censored Dependent Variable; the dependent variable may be censored from below, from above, or both._ In my case the dependent variable is censored by its rounded value, hence the censoring is data-dependent. Is it correct to use MCMCtobit() in this context ? EDIT2 : Now I have discovered the grouped R package. I think this is exactly what I was looking for. Thanks for the comments given below!


title: How to test (in Stata) whether the gender distribution of employees to jobs differs across two companies?

text: I have data on several companies where some are headed by a male CEO while others by a female CEO. As you can imagine, the jobs within these companies have different gender compositions. What I am trying to do is compare the gender distribution of employees across jobs for companies led by a male versus a female CEO. For example - Co A has a FEMALE CEO where Job1 is 40% female, Job2 is 30% female, Job3 is 15% female. Co B has a MALE CEO where Job1 is 60% female, Job2 is 40% female, Job3 is 20% female. Clearly Co B has a greater proportion of female employees, but my question is whether the distribution of female employees across jobs is significantly different for Co A and Co B. What is the appropriate statistical test and command in Stata for testing this? Assume these jobs are identical for simplicity.


title: Is it necessary to exclude all nonsignificant parameters to choose the best model?

text: I'm running LMM models and could I ask if I can just report the model after comparing random intercept with random intercept and slope model without excluding nonsignificant factor?


title: How to test if two regression models are significantly different in R

text: Generally speaking I want to find out if the regression models I built are significantly different from each other. I have monthly energy consumption data (the independent variable) from a couple of houses and Heating Degree Days (`HDD`) covering a period of about 3 years. The data has been divided into groups with Heating Season and House. For each house I have 3 (three heating seasons) regression models, and I would like to test if they are significantly different from each other. Maybe to make it less complicated **I would like to compare model of heating season`HS 1` with `HS 2`, and `HS 2` with `HS 3`**. On different forums I have read that I can use ANCOVA to compare models but not being advanced in statistics, I do not really know how to do it in practice. But I made myself familiar with `R` to a degree so I am trying to solve my problem with this software. I would appreciate your suggestions how to sort out this issue using `R`. Here is an extract of data for one house.               HDD    Energy   HeatingSeason    House     17.5       18            HS 1        1     84.4    206.5            HS 1        1     218.2   677.7            HS 1        1     374.4   824.4            HS 1        1     263.47  602.4            HS 1        1     174.27  514.4            HS 1        1     188.75  436.1            HS 1        1     58.59    44.3            HS 1        1     45.23       0            HS 2        1     26.14       0            HS 2        1     6.03        0            HS 2        1     6.95        0            HS 2        1     13.77       0            HS 2        1     51.82    45.1            HS 2        1     100.57  320.6            HS 2        1     208.78  634.5            HS 2        1     216.02  692.8            HS 2        1     275.34  643.2            HS 2        1     157.03  197.1            HS 2        1     165.42  246.7            HS 2        1     73.06    69.3            HS 3        1     46.21     3.9            HS 3        1     1.11        0            HS 3        1     17.15       0            HS 3        1     60.98    28.7            HS 3        1     164.28  197.1            HS 3        1     225.94  456.8            HS 3        1     315.3   687.8            HS 3        1     407.82  800.5            HS 3        1     342.88  771.2            HS 3        1     346.14  645.1            HS 3        1     191.03  223.7            HS 3        1


title: How do I compare stimulation levels between different stimulants with time points?

text: I have collected data that represent stimulation levels in three different time points ($T_0$, $T_1$ and $T_2$) using three different stimulants ($S_1$, $S_2$ and $S_3$) and a control ($C$). I have tested all three stimulants with three different subjects ($J_1$, $J_2$ and $J_3$) and for each subject and stimulant/ control I have three numbers that represent stimulation levels in $T_0$, $T_1$ and $T_2$ time points. How can I test the hypothesis that a certain stimulant shows significant higher level than the control? How can I test the hypothesis that all three stimulants shows a significant/ insignificant difference? The data is not normally distributed so I would have to use non-parametric tests. I am lost to whether it should be paired/unpaired or to use a regression or something else. Can anybody help?


title: Observed and estimated proportions differ with binomial GLM - which to present?

text: Using the `predict()` command in R, I have been able to calculate estimates (`av2`) and se-intervals (upper, lower)for my data below. The calculated proportion (`av2`) is different from the observed proportion (`av1`) - so 1. I'm not sure why this is? When writing up the experiment which should I quote and graph? Data and working are below.               A<-c(10,10,10,10,10,10,19,19,19,19,19,19)     B<-c("0","1","2","0","1","2","0","1","2","0","1","2")     C<-c("-ve","-ve","-ve","+ve","+ve","+ve","-ve","-ve","-ve","+ve","+ve","+ve")     Dead<-c(1,1,27,0,6,18,2,10,23,0,14,21)     Alive<-c(29,32,2,22,19,4,28,22,3,20,11,0)     Total<-Dead+Alive     gaf<-data.frame(A,B,C,Dead,Alive,Total)          mod2<-glm(cbind(Dead,Alive)~A+B+C, family=binomial)          p<-predict(mod2,newdata=gaf,se.fit=TRUE)     up<-with(p,fit+se.fit)     low<-with(p,fit-se.fit)     invLink<-family(mod2)$linkinv     av2<-with(p,invLink(fit))     upr<-invLink(up)     lwr<-invLink(low)      av1<-Dead/Total          round(av1, digits=2)     # 0.03 0.03 0.93 0.00 0.24 0.82 0.07 0.31 0.88 0.00 0.56 1.00 # observed           round(av2, digits=2)      # 0.01 0.11 0.83 0.02 0.19 0.90 0.03 0.33 0.95 0.07 0.49 0.97 # estimate           round(upr, digits=2)      # 0.02 0.15 0.88 0.03 0.26 0.94 0.06 0.40 0.97 0.12 0.57 0.98 # upper and lower se          round(lwr, digits=2)     # 0.00 0.08 0.76 0.01 0.14 0.86 0.02 0.26 0.92 0.04 0.41 0.96


title: Should I use missing data imputation with a model that already allows incomplete data?

text: I'm just starting to learn about missing data imputation methods, and I'm confused. In every introduction I've read, the author starts by describing listwise deletion and says that it's a bad idea because it reduces your N. They then explain methods you can use to fill in missing values so you can do your analysis with all your participants. That makes it sound like the purpose of imputation is to avoid throwing out other data points that you did observe. Question 1: Is imputation useful only because it lets you use the values you observed? Or do the imputed values themselves improve the analysis? For example, let's say I'm analyzing a repeated-measures dataset with a linear mixed model, and some participants are missing some timepoints. (Let's also say that the data are missing at random). Linear mixed models already work with incomplete data. Is doing imputation on the missing data points still appropriate, even though their missingness isn't causing me to throw out any other data? Question 2: What happens if I use full information maximum likelihood to impute values for the missing data points? It seems like that does provide some additional information for my model -- I get to include the portion of the missing data points that can be inferred from other variables in my dataset, even if those variables aren't part of the model I'm testing. But is that actually legitimate?


title: Family-wise confidence intervals

text: I have a bunch of variables organized into 10 different levels of a grouping factor. I'm doing some ANCOVA on particular variables and also plotting the data using boxplots. I'd like to add 84% confidence intervals to all the groups (since non-overlapping 84% CIs indicate a significant difference at alpha .05 - at least for two groups). I can do all this quite easily in R. My question is - should I be applying a "family-wise" 84% CIs to all the groups? In other words, just as one would devalue an alpha level by the number of groups to obtain a family-wise alpha, should I inflate the CI a reciprocal amount to achieve a family-wise interval? This seems reasonable to me, but I haven't seen this discussed in the literature. If alpha were CI were interchangeable for two or more groups the the family- wise 84% CI would be 99.5%, but i've read that alpha and CI are only interchangable for 1-sample situations. If this is the case, how would I go about calculating the family-wise confidence intervals for 10 (or any number) groups? Any advice would be welcome. best, Steve


title: K-Fold Cross validation and F1 Measure Score for Document Retrieval using TF-IDF weighting and some customised weighting schemes

text: I am developing a search engine system based on the vector space model, and I am confused on what approach I should take to evaluate the system. My case is this:   1. I have a set of indexed documents in pairs (di,ti) stored in a database. And I can calculate the similarities between a query and the set of indexed documents in the system using tf-idf weighting scheme and the cosine similarity measure and rank the document in descending order of the similarity score. How do I evaluate the system using k-fold cross validation and the F1 measure? In this case what would my Training Set, Validation Set and Test Set be ? I am confused which dataset to split into these three. What data do I use for these three sets in terms of query terms, terms, documents and similarities ? As I understand it, I will have to split the dataset into: **Test Set** : A set of query terms **Training Set:** A set of documents/term pairs **Validation Set:** A set of similarity/document pairs ? How can I combine the F1 Score with the K-Fold cross validation ? Sorry if my question is ambiguous, but I don't completely understand these evaluation methods and I am new to all of these. I implemented the system with plain python and django for interacting with the database, I am **not** using any libraries such as numpy, scikit learn. **Update 1:** I will try to explain my system further since I have not included all the details and that's why I am confused between the two. I have another case where I have a set of tags for each of the document, so there is a set of tags assigned to those documents, by some users. So my searching algorithm works for two different cases:   1. The first is as described above where the tf-idf weighting scheme is used and the cosine similarity measure to calculate the similarity between some queries and the documents.   2. The 2nd approach is: I calculate some weights based on the tagging behaviour on the documents, and I am adding these weights to the tf-idf weights i.e if the the term and tag is the same then I add their weights together, if the document doesn't have a term which is the same with the tag then the addition would be 0 + the tag/document weight (since there is no such term in the document term/document relation, and thus its term/document weight would be zero) Now I want to evaluate the 1st and 2nd case and compare the rankings, so for the 2nd case I am assuming I can use k-fold ? What approach would you suggest to evaluate the ranking of the results for such a scenario?


title: What is a statistical journal with quick turnaround?

text: What are good statistical journals with quick turnaround (fast review cycle), suitable for short notes in mathematical statistics and preferably with open access. An example is Statistics & Probability Letters, however, that journal only has sponsored open access.


title: Bayesian probability > 1 -- is it possible?

text: Running a crude model using Bayesian inference, I get some results > 1 (ie, more than 100% "certain") for some combinations of "evidence". For instance, for one bit of evidence the conditional probability of the null hypothesis is 0.85 while the marginal probability is 0.77. If the prior probability is 0.9, the computed posterior probability is 1.008. Which maybe could be ascribed to rounding error, except that the next bit of evidence raises the posterior probability to 1.34. It stands to reason that, for a problem with two hypotheses, one would have a conditional probability less than the marginal probability and the other would be greater. So the resulting P(E|H) / P(E) multiplier would be > 1. So it's hard to see how such > 1 results can be avoided in the general case. Is this just the way Bayesian inference works, or do I likely have an error somewhere in my calcs? ## Data               One "evidence":     Total 6134 samples     Total actual positives 2845     Total actual negatives 3289     Test was true 1623 times     Test was true 465 times when the "gold standard" was positive     Test was true 1158 times when the "gold standard" was negative          conditional probability of positive hypothesis = 465/2845 = 0.1634     conditional probability of negative hypothesis = 1158/3289 = 0.3521     marginal probability of test being true = 1623/6134 = 0.2646     Bayes multiplier for the negative hypothesis = 0.3521 / 0.2646 = 1.3307      As can be seen, the multiplier is significantly > 1, and with several such tests back-to-back it seems hard to avoid probabilities > 1. (Of course, I suppose one can argue that the tests aren't truly independent, and that puts the fly in the ointment.) ## Fudging So does anyone have any suggestions as to how to "fudge" non-independent measurements to improve an estimate? For the two most egregious cases I can come up with a fair estimate of how connected the measurements are, but I don't have a feel for how to factor that knowledge in.


title: 2D analog of standard deviation?

text: Consider the following experiment: a group of people is given a list of cities, and asked to mark the corresponding locations on an (otherwise unlabeled) map of the world. For each city, you will get a scattering of points roughly centered at the respective city. Some cities, say Istanbul, will exhibit less scattering than others, say Moscow. Let's assume that for a given city, we get a set of 2D samples $\\{(x_i, y_i)\\}$, representing the $(x, y)$ position of the city (e.g. in a local coordinate system) on the map assigned by test subject $i$. I would like to express the amount of "dispersion" of the points in this set as a single number in the appropriate units (km). For a 1D problem, I would choose the standard deviation, but is there a 2D analog that could reasonably be chosen for the situation as described above?


title: How to re-scale a five point likert scale into three catagories

text: Can anyone help me in re-scaling a five point likert scale into three categories using the mean scores. I have used a 5 point likert scale:                (1) Not at all, (2) low, (3) moderate, (4) high and (5) very high        Now I want to categorize these into `Low`, `Moderate` and `High`, but I am not sure how to use use the result. I planned try it like this:   1. less than 2.5 = `Low`   2. between 2.5 and 3.5 = `Moderate**   3. greater than 3.5 = `High`


title: What are the statistical assumptions, prior to running hierarchical regression to test interaction effects?

text: What are the assumptions that need to be met when using hierarchical regression, and subsequently simple slopes analysis to test and probe interaction effects and how can I test for them in SPSS? I have a continuous DV and IV, with a categorical moderator (3 categories). Thanks in advance.


title: Why doesn't backpropagation work when you initialize the weights the same value?

text: Why doesn't backpropagation work when you initialize all the weight the same value (say 0.5), but works fine when given random numbers? Shouldn't the algorithm calculate the error and work from there, despite the fact that the weights are initially the same?


title: Standard error of the median

text: Is the following formula right if I want to measure the standard error of the median in case of a small sample with non normal distribution (I'm using python)?                sigma=np.std(data)      n=len(data)      sigma_median=1.253*sigma/np.sqrt(n)


title: How to calculate class-specific accuracy on PCA?

text: I have implemented PCA at Matlab. I just did dimensionality reduction (I don't use but I have class labels in my data set). My aim is to calculate overall accuracy and class-specific accuracy. I know that I can calculate overall (training data set and test data set) accuracy with PoV (Proportion of Variance Explained)--if I am wrong you can explain it, you are welcome--but I don't know how to calculate class-specific accuracy. It sounds to me that each class accuracy is the same as overall accuracy.


title: R: mlogit error on data - "system is exactly singular"

text: So I have data from a randomized blind trial of 1mg of nicotine gum on dual n-back working memory scores; I analyzed them as usual with a t-test and found a small increase in means but a large increase in standard deviations on a f-test! Strange. I also have data for each day on mood/productivity that day on a 1-5 scale. I wondered: is nicotine following an inverse U-curve, where it causes higher scores on the worser days (1-3) and lower scores on the better days (3-5)? I look around and it seems I want a multinomial logistic regression comparing the placebo & active days. I enter the data & load `mlogit`:               nicotine <- read.table(stdin(),header=TRUE)     day      active mp score     20120824 1      3  35.2     20120827 0      5  37.2     20120828 0      3  37.6     20120830 1      3  37.75     20120831 1      2  37.75     20120902 0      2  36.0     20120905 0      5  36.0     20120906 1      5  37.25     20120910 0      5  49.2     20120911 1      3  36.8     20120912 0      3  44.6     20120913 0      5  38.4     20120915 0      5  43.8     20120916 0      2  39.6     20120918 0      3  49.6     20120919 0      4  38.4     20120923 0      5  36.2     20120924 0      5  45.4     20120925 1      3  43.8     20120926 0      4  36.4     20120929 1      3  43.8     20120930 1      3  36.0     20121001 1      3  46.0     20121002 0      4  45.0     20121008 0      2  34.6     20121009 1      3  45.2     20121012 0      5  37.8     20121013 0      4  37.2     20121016 0      4  40.2     20121020 1      3  39.0     20121021 0      3  41.2     20121022 0      3  42.2     20121024 0      5  40.4     20121029 1      2  41.4     20121031 1      3  38.4     20121101 1      5  43.8     20121102 0      3  48.2     20121103 1      5  40.6          library(mlogit)     Nicotine <- mlogit.data(nicotine,shape="wide", choice="mp")     mlogit(score ~ (active + mp)^2, Nicotine)     Error in solve.default(H, g[!fixed]) :        Lapack routine dgesv: system is exactly singular     Calls: mlogit ... mlogit.optim -> as.vector -> solve -> solve.default      The error also happens even with the simplest call I can think of:               mlogit(score ~ active, Nicotine)     Error in solve.default(H, g[!fixed]) :        Lapack routine dgesv: system is exactly singular     Calls: mlogit ... mlogit.optim -> as.vector -> solve -> solve.default      Reading the documentation for `mlogit` didn't much help, and look at the other questions having the same error, they're different enough I can't tell whether they apply or not. Thank you for your assistance.


title: Regsubsets with leaps fails

text: I am trying to fit a problem with `regsubsets` with leaps in R. My problem is particularly strongly collinear, which is why I chose to use it in the first place. The number of variables is about 200 and I have about 2 million independent observations. All the variables have a strong correlation structure with each other. On running regsubsets with `really.big = TRUE`, and `nvmax = 5` and `nbest = 1`, I get the following: > Error in leaps.setup - 31 linear dependencies found and it crashes. All I am looking to do is a simple forward stepwise, say order the variables in the order of correlation and run nested regressions. Is that too much for the software to handle? I think the problem is well posed in that sense.


title: How to test for causation in a static panel data model with small t?

text: I have a static panel data model with small T (T=5) that makes it impossible for me to use granger causality as it requires a long time span. So my question:   * Is there any alternative solution to test for causation even in a small T context? Any hint will be highely appreciated!


title: An easy explanation for the parallel coordinates plot

text: I have read and seen a lot of Parallel coordinates plots. Can someone answer the following set of questions:   1. What are parallel coordinates plots (PCP) in simple words, so that a layman can understand?   2. A mathematical explanation with some intuition if possible   3. When are PCP useful and when to use them?   4. When are PCP _not_ useful and when they should be avoided?   5. Possible advantages and disadvantages of PCP


title: What is the difference between summary and loadings for princomp?

text: Example code:               (pc.cr <- princomp(USArrests))  # inappropriate     summary(pc.cr)     loadings(pc.cr)  ## note that blank entries are small but not zero      I am getting different outputs from each, and I am not sure I understand what the difference is. Here is the output:               > summary(pc.cr)     Importance of components:                                Comp.1      Comp.2      Comp.3       Comp.4     Standard deviation     82.8908472 14.06956001 6.424204055 2.4578367034     Proportion of Variance  0.9655342  0.02781734 0.005799535 0.0008489079     Cumulative Proportion   0.9655342  0.99335156 0.999151092 1.0000000000               > loadings(pc.cr)  ## note that blank entries are small but not zero          ...                         Comp.1 Comp.2 Comp.3 Comp.4     SS loadings      1.00   1.00   1.00   1.00     Proportion Var   0.25   0.25   0.25   0.25     Cumulative Var   0.25   0.50   0.75   1.00      p.s: how can I get access to the table created by summary(pc.cr)?? (I can't seem to find it in str) Thanks.


title: KPSS test for large samples

text: I am trying to interpret some KPSS results for a large sample (n=1800 - http://data.is/TB2z5b). In general the KPSS test can be used for samples with time series with more than 100 observations? Is there a reference I can consult?


title: Interpreting decision trees outputs in R?

text: I have created 2 decision trees, I just want to assure if I am making correct interpretation about it. Here is my first tree: Note: Right click on image and select view image to view it clearly! ![tree-1](http://i.stack.imgur.com/OhasY.jpg) Here is interpretation: Let us assume a value of ipttl be 40. The tree shall start by checking the given value of ipttl for this condition at root ie. ipttl<44.5, now as 40<44.5 so the control moves to left branch from root. Now it will check the ipttl value for condition ipttl>=43.5, so as 40>=43.5 condition is false so the controls moves to right branch now it will be checked for condition ipttl>=42.5 so as 40>=42.5 condition is false so control shall move to right branch and it means this ipttl value of 40 belongs to category 8. Is my interpretation correct? I managed to create another decision tree in R. But the tree seems to be quite populated at bottom and I am unable to judge it. moreover why is it having same value of date_time everywhere in tree. ![tree-2](http://i.stack.imgur.com/PmPXB.jpg) My R-output is as follows:               > tree     n=5544 (756 observations deleted due to missingness)          node), split, n, loss, yval, (yprob)           * denotes terminal node           1) root 5544 4829 6 (0.13 0.13 0.12 0.12 0.13 0.13 0.13 0.13)          2) date_time< 1.351177e+09 2727 2025 1 (0.26 0.26 0.24 0.25 0 0 0 0)            4) date_time< 1.351176e+09 1398  696 1 (0.5 0.5 0 0 0 0 0 0)              8) date_time< 1.351176e+09 702    0 1 (1 0 0 0 0 0 0 0) *            9) date_time>=1.351176e+09 696    0 2 (0 1 0 0 0 0 0 0) *          5) date_time>=1.351176e+09 1329  660 5 (0 0 0.5 0.5 0 0 0 0)             10) date_time< 1.351176e+09 660    0 3 (0 0 1 0 0 0 0 0) *           11) date_time>=1.351176e+09 669    0 5 (0 0 0 1 0 0 0 0) *        3) date_time>=1.351177e+09 2817 2102 6 (0 0 0 0 0.25 0.25 0.25 0.25)            6) date_time< 1.351177e+09 1417  702 6 (0 0 0 0 0.5 0.5 0 0)             12) date_time< 1.351177e+09 715    0 6 (0 0 0 0 1 0 0 0) *           13) date_time>=1.351177e+09 702    0 7 (0 0 0 0 0 1 0 0) *          7) date_time>=1.351177e+09 1400  697 8 (0 0 0 0 0 0 0.5 0.5)             14) date_time< 1.351178e+09 703    0 8 (0 0 0 0 0 0 1 0) *           15) date_time>=1.351178e+09 697    0 9 (0 0 0 0 0 0 0 1) *      The main question is how can I know if the tree is able to classify my target variable (category) properly as it is having only one value of date_time in tree plot? Regards,


title: Compare means of 2 variables (no individual data available)

text: **Short question:** I want to find out whether the means of 2 variables differ significantly. I know how large the dataset is (N~2000), but I only have the mean scores averaged over all participants. I do not have information about standard deviation. Is there a test that can do that? **Details:** 2 variables x1 and x2 were measured at times t1 and t2. Both increase. I need to compare _relative_ change over time. The relative change for x1 is x1t2 / x1t1. Since the huge majority of subjects at t1 have values of 0 on my x, calculating relative increase is impossible for these (and I would lose 90% of the data). So what I do instead is averaging x1 and x2 scores over the whole population at both measurement points, and then calculate the relative increase of these mean scores. Example: Assume the mean of x1t1 to be 0.5 and the mean of x1t2 to be 2.0. The relative increase is 4 (or 500%). Assume the relative increase of x2 to be 5 (or 600%). Now I want to test whether 4 and 5 differ from each other significantly.


title: What does "degrees of freedom" mean in the context of cointegration?

text: I recently read a paper where the author says they used cointegration "with two degrees of freedom". What does degrees of freedom mean in this context?


title: Information gain as a feature selection for 3-class classification problem

text: I am facing a sentiment analysis task where I am using Naive Bayes to classify documents as Positive, Negative or Neutral. I have thought of using Information Gain as my filter for feature selection. Taking into account that I need to classify into 3 classes, would it be ok to use IG without further modifications (since I think IGain works best with 2 classes)?


title: The best way to solve particular classification problem?

text: I got training set (time series) of size approximately 2 million precedents {x,y}. Each x is a vector of size 20 and each y is a binary vector of size 10 like {1,0,0,1,1,0,1,1,1,0}. For a new input x I would like to obtain vector of probabilities {P(y_i = 1 | x)}. So, what would you recommend?


title: Nagelkerke pseudo-R2 with positive log likelihoods

text: I'm trying to calculate a pseudo-R2 for linear mixed models using Nagelkerke's method . My understanding is that Nagelkerke's pseudo-R2=1-EXP[(-2/n)(l(B)-l(0))], where l(B) and l(0) are the maximized log-likelihoods for the fitted model of interest and the null (intercept only) model, respectively. The pseudo-R2 is then rescaled to a maximum of 1 by dividing by R2(max), where R2(max)=1-EXP[(2/n)*l(0)]. My problem is that in some cases the log-likelihood for the model of interest (and sometimes also the null model) is positive. Consequently, R2>R2(max), leading to rescaled R2 values that are either negative or >1. Am I doing something wrong here? If not, is there a different method for calculating reasonable pseudo-R2 for models with positive log likelihoods? Thank you, in advance, for any advice. Laura Perry


title: Computing average of data in multiple (different) columns in SPSS

text: I have a SPSS data file, which i am trying to reduce. However the data which belongs together is spread within multiple columns. I.e. I have one row per subject, but each subject has done multiple different conditions as well as procedures. The conditions are in columns `condition_1`, `condition_2` etc. The procedures are in rows `procedure_1`, `procedure_2` etc. The values (reaction times) are in columns `reaction_time_1`, `reaction_time_2`. Now I want to compute the average of all reactions times for each person where the condition and procedure is the same. So if `condition_x` is the same as `condition_y` and `procedure_x` is the same as `procedure_y` I want to get the average over all of those items `reaction_time_x` and `reaction_time_y`. How can I do this with SPSS? I easily know how to compute averages in general, but in this case I have to select only a subset of the rows based on the conditions, and I do not know how to do this.


title: Code for fitting multivariate t distribution (ECME)

text: I am looking to fit some data with central multivariate t distributions. I understand an ECME algorithm is the most promising way to do this. I wondered if anyone knew of any existing code for fitting t distributions with ECME that would get me started (in any language)?


title: Quantifying applicability domain for predictive models?

text: I have a big dataset and I want to build a classification model (svm, rf, ann etc.). Then I split the original dataset into training set and test set. I build the model using training set. After it was done, I use the model to predict the test set. Here, the training set and test set should be similar in data structure, that means they are located in the same data domain, right? If the test set is not located the data domain of the training set, then the prediction result for the test set should be worse than those in the train set domain. it is right to express this? So, how can I quantify the applicability domain for a predictive model? that means, when I have the new dataset, how can I check if it located in the applicability domain of the model? some links: http://en.wikipedia.org/wiki/Applicability_domain http://www.moleculardescriptors.eu/tutorials/T7_moleculardescriptors_ad.pdf Any points are appreciated.


title: Model prediction for every possible outcome

text: Let's say I'm running a model to look at questions asked on a form. My model looks like: Price = B0 + B1(Gender) + B2(Age) + B3(Credit Type) + B4(Home) ... I'm trying to calculate expected price for every possible scenario.               So price for gender = 0, age = 18, home = 1, etc.     Price for gender = 1, age = 18, home = 1, etc.     ... and so forth      In R, is there a way to get the predicted price (output) for all possible combinations of the independent variables.               library(car)     mod = lm(sold ~ age + gender + marital + educ + cars + license + credit + type + home + id)     summary(mod)      EDIT: After reading some of the comments below, I realized the stupidity of what I was doing. But I'm still wondering if there's a way to find the expected price. A general measure of what would the price be given the values/variables in this model. So lets that that from the model, is there a way to glean the expected price. Let's say that it's $10 given the model. EDIT 2: Actually, just to figure it out in R, can anyone point to how I can perform the original task. The predicted value for all variations of the predictors.


title: Performance metrics to evaluate unsupervised learning

text: With respect to the unsupervised learning (like clustering), are there any metrics to evaluate performance?


title: Characterizing rank along "lower envelope" of points

text: Suppose I have data that looks like this. ![enter image description here](http://i.stack.imgur.com/vxlIc.png) I want to quantify how well the "lower envelope" of the points are described by a monotonic function. The measure will be used to compare sets of points and should prefer sets with more closely monotonic lower envelopes. Emphasis should be placed on individual outliers of low y and high x that break from monotonicity (red points). By "lower envelope" I roughly mean the points circled in yellow. ![lower envelope](http://i.stack.imgur.com/7voEL.png) How can I turn what I want to characterize into a number? **My thoughts on how I might do this:** For defining the lower envelope of the points, I could use conditional kernel density estimation. The lower envelope would then be defined as some bottom % of y for given x. Following this procedure I could then either characterize the points that are within that % cutoff or characterize the density distribution itself. Characterizing the points seems straightforward. I was thinking of Spearman's rho. This approach has the advantage of being simple however the data is such that some regions of x values may be heavily sampled and others may be more sparse. I'm not sure if that would cause a problem. Characterizing the conditional density estimate might be preferable but I have no idea how to get a sense of rank in y given density estimates along x. One way would be to pick a % y for the conditional density estimate along x and calculate rank from points along x at that % y. Is there a better way to do this?


title: Learning from one positive

text: Say that in a binary classification problem you have several negatives and only one positive. What types of models are good to learn from this data, and predict the label for a new instance? Anything available in Python? From what I read here, scikit-learn doesn't support it.


title: Is it possible for my mixed model results to be the same as those for a repeated measures ANCOVA? Why don't I get the same results? SPSS

text: **The experiment:** Two groups of children perform tasks with three levels of difficulty at two different time points and brain activity is measured. All of the children perform all of the tasks and do so at both timepoints. **Dependent variable:** Brain Activity (BA)   **Repeated measures:** task (three levels), time (2 levels)   **Independent:** Group   **Covariate:** Age I have been using a repeated measures ANCOVA to analyse the data but recently have attempted a mixed model for increased flexibility in contrasts / post hoc options. A 3 way interaction that is meaningful to me, `group x task x time`, is significant with the ANCOVA method but not with the mixed model, even when modeling the same main effects and interactions. The syntax (SPSS) for each model is below. I am concerned that I am using the mixed model incorrectly. Does anyone know what I could do differently with my mixed model to make it more similar to what I have done with the ANCOVA?               DATASET ACTIVATE DataSet1.      GLM Time1Task1 Time1Task2 Time1Task3 Time2Task1 Time2Task2 Time2Task3 BY Group WITH Age        /WSFACTOR=fixation 2 Polynomial target 3 Polynomial        /METHOD=SSTYPE(3)        /CRITERIA=ALPHA(.05)        /WSDESIGN= task time task*time        /DESIGN= Group Age.          MIXED EEG BY Group Time Task WITH Age        /CRITERIA=CIN(95) MXITER(100) MXSTEP(10) SCORING(1) SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)        /FIXED=Group*Time*Task Task Time Time*Age Group*Task Group*Time Task*Age Time*Task Time*Task*Age | SSTYPE(3)        /METHOD=REML        /RANDOM=INTERCEPT Age | SUBJECT(id) COVTYPE(UN)        /REPEATED=Target*Time | SUBJECT(id) COVTYPE(UN).      In the ANCOVA, the tests of within subjects effects table ends up having 9 effects. I have specified the same 9 effects in the mixed model. I realize this is a rather ridiculously large number of interactions included in the mixed model, but I wanted to mirror what the ANCOVA gives me (it is not possible to get an output for `task x time x age` in the GLM unless I include all of the above interactions - the model is only customizable within "within subject" factors and within "between subject" factors and then it automatically includes all interactions between "within" and "between" subject factors. I have played around with:   * not including random effects (age, subject ID)   * including less effects in the mixed model while still including the interaction of interest, task x time x age   * different covariance matrices And none of these things mean that the `task x time x age` interaction is significant. I am concerned that there is something fundamentally wrong about the way I have depicted repeated measures testing with a covariate in the mixed model rather than these minor adjustments.


title: Conditional logistic regression vs GLMM in R

text: I have paired data (GWAS case/control study) and I have heard using conditional logistic regression or generalized linear mixed models (GLMM) is appropriate. Which should I use in this case? Why would you use one over the other. More importantly can you guys point me towards resources for doing these methods in `R`? I'm finding a lot of material for `SAS`, which I do not prefer. I can provide more details if necessary.


title: Interpreting regression coefficients of log(y+1) transformed responses

text: I have measurements $y_1$,...,$y_i$,...,$y_n$ taken from a set of replicates in a factorial designed experiment. In order to use a linear regression I define my response $z_i = log(y_i + 1)$. The log is used in order to make the normality assumption hold and the model fit, and the + 1 is used since some $y_i$'s are 0. If the model is ${z = a + \beta X}$, then the interpretation of $\beta$ is ${\beta = {\sum_{i=1}^nz_i}/n = {\sum_{i=1}^nlog(y_i + 1)}/n}$ so ${{\log(\pi_{i=1}^n(y_i + 1))} = \beta n}$, which gives that the geometric mean $\sqrt[n]{\pi_{i=1}^n(y_i + 1)} = e ^\beta$. My question is whether there is some back transformation or some other way to get an interpretation of the geometric mean of $y_i$'s rather than of $(y_i + 1)$'s as a function of $\beta$. The reason I'm asking is that for cases where $y_i$'s are small (close to 0) the value of $\beta$ underestimates the magnitude of the effect.


title: R Using Random Forest or SVM

text: # Problem: **The following code is to predict the last row's price using all previous rows as a training set** ** The prediction is very nice when the last row price smaller than 10, but becomes unreasonable when the test price gets bigger. ** Same problem happens using SVM package               file <- read.csv('test.csv',head=F)          N<-nrow(file)     train <- file[1:N-1,]     test <- file[N,]          forest <- randomForest(train[,1] ~., data = train)     rf.predict <- predict(forest, newdata=test)          rf.predict  #10.08446      data:                      V1     V2     V3     1     3.0    1.0    2.0     2    -6.0   -6.2   -6.2     3   -30.2  -31.2  -30.7     4    -2.1   -1.3   -1.8     5   -11.4   -9.5  -10.7     6     8.4    8.8    8.6     7     6.3    6.5    6.4     8    -0.2    3.2    1.5     9    -8.8  -10.2   -9.3     10  -10.2  -11.0  -10.7     11  -16.3  -15.6  -16.0     12   -2.5   -2.5   -2.5     13  -17.2  -16.7  -16.8     14   10.4   10.0   10.2     15   -3.0   -5.3   -4.0     16    2.1    3.5    2.5     17   -2.0   -1.4   -1.5     18   11.7   12.2   11.8     19   11.2   11.2   11.2     20   -1.0   -0.4   -0.9     21    3.8    1.9    3.2     22   -9.3  -10.7   -9.7     23  -17.7  -15.7  -17.5     24 1000.5 1000.0 1001.0


title: The effect of skewed continuous predictors in a binary logistic regression model

text: I am analyzing data with a binary outcome and a variety of continuous and categorical (including dichotomous) predictor variables. My approach is to perform a binary logistic regression and to treat any predictor with more than 20 unique values as continuous. Several arguments against categorization, especially well-documented on Frank Harrell's site are a good reason to not categorize. However, at a recent meeting where I discussed my analyses approach, a faculty suggested that I would get more accurate risk estimates for the data if I categorized the variables which have a skewed distribution and outliers. Their logic was that the tail of the skewed distribution and the outliers in that tail will have a detrimental effect on the risk estimates generated by the logistic regression and that categorization will address this by erasing the effect of the tail and the skew. I have several predictor variables that definitely have skewed distributions and some outliers. Is the claim that a variable with a skewed distribution (with outliers) is more likely to produce inaccurate risk estimates compared to the categorized version of the same variable, true? How do skews and outliers in the tail affect logistic regression estimates?


title: conditional expectations value

text: I need to calculate the following integral $$\int_{\mu+c}^{\infty} y\cdot \frac{1}{\sigma\sqrt{2\pi}}e^{(y-\mu-w)^2/2\sigma^2}dy$$ So essentially $y\sim N (\mu+w, \sigma^2)$ and im trying to calculate its expected value from $\mu+c$ to $\infty$ I know I will have to use integration by parts, but don't know how to proceed. More importantly however, when I calculate the integral, will I get an expression which includes $w$? I really need the $w$ in there.


title: A non-negative definate matrix has a non-negative generalized inverse

text: I'm having trouble proving a N.N.D matrix has a N.N.D G-Inverse. So far I have: If we assume x = Az where x >= 0 and A is a nnd matrix. So if Y is a G-inverse than: x = Az = YAz = Yx >= 0 . Thus Lx = LAz = Az >= 0 . So L is any G-inverse of A is also nnd.


title: Univariate fixed effect Vs Multivariate model -Negative Covariance, positive parameter estimate, but why?

text: I am trying to compare the results of two models.   The first model looks at y with x as a fixed effect. The second looks at the covariance between x and y. Both models have repeated measures for x and y, so individual is included as a random effect. I am running the models in MCMCglmm. I am running them as Gaussian for ease in this example and the residuals suggest this is an ok approximation. data is called d1               library(MCMCglmm)     head(d1)          d1$Indiv <-as.factor(d1$Indiv)          (d1$Y, d1$X) # suggests a positive relationship between X and Y      these are all starting priors and not optimised but look ok on first checks               prior1 = list(R = list(V = diag(1),nu=0.002),      G =  list(G1 = list(V=diag(1), nu=1, alpha.mu=c(0), alpha.V=diag(1)*1000)))          M1 <- MCMCglmm(Y~X,     random= ~Indiv,     rcov=~units,     data = d1,     prior = prior1,      nitt= 100000, thin =200, burnin = 100)      this gives a positive parameter estimate for X, suggesting a positive slope.               prior3= list(R = list(V = diag(2),nu=0.002),      G =  list(G1 = list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*1000)))          M2 <- MCMCglmm(cbind(Y,X)~1,     random= ~us(trait):Indiv, ## could run with just Indiv                  rcov=~us (trait):units,     data = d1,     family=c(“gaussian”,”gaussian”),     prior = prior3,     nitt= 100000, thin =200, burnin = 100)          posterior.mode(M2$VCV)      this gives a negative covariance/correlation then I thought that maybe if the variance was difference for X this could be causing problems so I standardised X               d1$stX <-( d1$X -mean(d1$X))/sd(d1$X)     M1.1 <- MCMCglmm(Y~stX,     random= ~Indiv,     rcov=~units,     data = d1,     prior = prior)      this gives a positive parameter estimate for stX, suggesting a positive slope. So I thought the difference between the models could be because there is an individual effect on X and Y individually in the covariance model. So I mean that the fixed effect model allows Y to have an individual effect but not a unique one for X. So I fitted a random slope, which in theory should account for some differences in the variance of X (I say this as a random slope alters the X intercpt)               M1.2 <- MCMCglmm(Y~stX,     random= ~us(1 + stX):Indiv,     rcov=~units,     data = d1,     prior = prior)      this gives a positive parameter estimate for stX, suggesting a positive slope. I have looked at other things like the correlation between blups from a model of: Y~1 , random ~ Indiv and X Y~1 , random ~ Indiv - still get a positive relationship. I am bemused. I have tried simulating data with a set covariance structure and trying to generate large differences in the effect if individual on X and Y to see if I can recreate this effect, but I struggling to create a situation with a strong covariance and large individual effects on only one trait, suggesting this is not the cause of these discrepancies. I wonder if the autocorrelation in X could be a problem? For Yt, X = Bt-1 + Bt-2 So Yt =5, X = B4 + B3 So Yt =6, X = B5 + B4 This was potentially an error in the design and the calculation of X – I can change it very easily it but I am intrigued- could this temporal correlation reverse the direction of a covariance? I am really keen to better understand how the model works to be sure I understand if something like this happens again Data below – sorry I cropped it but you need a lot to run the model correctly Many thanks Sam Indiv Y X 1896 0 0 1685 2 1 2021 0 1 2054 0 1 2093 1 1 2219 3 1 201 2 2 1201 0 2 1340 2 2 1363 2 2 1579 0 2 1596 2 2 1667 3 2 1677 3 2 1990 2 2 2064 0 2 2165 2 2 2178 2 2 2207 3 2 2219 3 2 2240 3 2 2240 3 2 2626 3 2 1 3 3 970 3 3 1141 1 3 1149 1 3 1149 2 3 1201 1 3 1201 0 3 1202 1 3 1243 2 3 1252 3 3 1252 3 3 1336 1 3 1336 2 3 1348 2 3 1448 1 3 1490 0 3 1504 2 3 1580 2 3 1621 2 3 1637 1 3 1654 2 3 1670 2 3 1692 2 3 1692 2 3 1817 3 3 2004 2 3 2011 1 3 2021 0 3 2048 3 3 2055 3 3 2089 1 3 2092 2 3 2110 1 3 2118 2 3 2121 3 3 2178 1 3 2191 1 3 2231 1 3 2527 0 3 258 3 3.5 1 1 4 1 2 4 167 3 4 212 3 4 970 2 4 1052 0 4 1141 1 4 1202 1 4 1205 2 4 1246 1 4 1445 2 4 1450 1 4 1467 2 4 1488 1 4 1568 3 4 1580 2 4 1580 0 4 1620 2 4 1629 2 4 1630 2 4 1635 1 4 1643 2 4 1645 3 4 1649 0 4 1681 3 4 1698 2 4 1703 3 4 1703 3 4 1716 3 4 1795 2 4 1795 2 4 1795 2 4 1797 1 4 1801 2 4 1814 3 4 1823 2 4 1885 1 4 1900 2 4 1913 2 4 1979 1 4 1979 2 4 1982 3 4 1986 2 4 1986 2 4 1986 0 4 2004 2 4 2011 1 4 2021 0 4 2058 2 4 2058 2 4 2096 1 4 2191 2 4 2217 1 4 2227 3 4 2269 3 4 2351 2 4 2385 2 4 2544 2 4 2545 1 4 2623 2 4 2629 2 4 1492 2 4.5 2303 2 4.5 1 2 5 167 2 5 1202 1 5 1208 2 5 1231 2 5 1237 2 5 1240 1 5 1345 2 5 1448 2 5 1492 1 5 1492 0 5 1568 3 5 1570 2 5 1630 2 5 1645 3 5 1646 1 5 1681 3 5 1823 2 5 1879 2 5 1978 2 5 2011 1 5 2014 2 5 2030 2 5 2063 1.5 5 2095 2 5 2240 2 5 2526 1 5 2633 2 5 1448 1 5.5 1208 2 6 1250 2 6 1338 3 6 1338 3 6 1462 2 6 1784 3 6 1787 3 6 1948 1 6 2014 1 6 2015 2 6 2015 1 6 2023 2 6 2030 2 6 2109 2 6 2218 2 6 2267 2 6 2521 2 6


title: Unknown formula

text: I came across a formula presented in a financial context, and have been trying unsuccessfully to try a work out what it means. $SL = 2 * 100 * \sqrt{V_L + (\alpha + \beta)^k \sigma^2_{n - V_L}}$. I have been scouring the net and can't see anything that comes close, so any help would be appreciated. It's driving me nuts


title: What distribution family would best fit this graph?

text: Sorry if this type of question is not kosher. I'm new around here, so please forgive me. Anyway, I have a dataset that describes the probability that users will like certain articles from my corpus _(with 5 articles chosen randomly and graphed below. Along the x axis are individual users, and along the y-axis their score. The users are sorted by score.)_ I'm looking to fit a distribution to each line, but am unsure of which one to choose that would likely best these lines. I was thinking of using an exponential distribution? But there's a lot of information that exponential probably would not capture here. Does anyone have a good idea? EDIT: I will update my title to include the name of a reasonable distribution so that this question will be more searchable and others would possibly learn from it in the future. ![enter image description here](http://i.stack.imgur.com/9OAix.png)


title: Confusion about error term in ARMA-model when predicting future values

text: Let's say I'm fitting data for example to ARMA(1,1)-model: $x_t = \phi x_{t-1} + \epsilon_t + \theta \epsilon_{t-1}$. Now I estimate the parameters $\phi$ and $\theta$ and solve some values for them, e.g. $x_t = 0.7 x_{t-1} + \epsilon_t + 0.8 \epsilon_{t-1}$ (I just made the numbers up) Let's say this the best model for my problem at hand and I begin forecasting. So I make the forecast $x_{t+1} = 0.7 x_{t} + \epsilon_{t+1} + 0.8 \epsilon_{t}$. Now where I get confused is with the future value of the error term $\epsilon_{t+1}$. How can I solve this value? Should I predict the residuals also or should I use the expected value of $\epsilon_{t+1}$, which should be 0. Hope my question is clear =) Thank you for any help =) EDIT: About the value of $\epsilon_{t+1}$. Should I estimate the value of $\epsilon_{t+1}$ by assuming (as in literature normally is assumed) that the noise process $\epsilon_t$ is normally distributed $\epsilon_t$ ~ $iidN(0,\sigma_{\epsilon}^2)$ and then use estimation techniques (Least squares, Maximum likelihood, Yule-Walker) to estimate the value for noise process variance $\widehat{\sigma}_{\epsilon}^2$ and then just evaluate value for $\epsilon_{t+1}$ ~ $iidN(0,\widehat{\sigma}_{\epsilon}^2)$ from the estimated Gaussian distribution?


title: Prediction Intervals in Ridge Regression?

text: For some reason, I can't seem to find the formula for prediction intervals in ridge regression anywhere. I know that the coefficient estimates are biased, but are the predictions (dependent variable) biased as well?


title: Datamining and time series forecasting

text: Could we say time series forecasting is a part of data-mining or it's just a data-mining tool?


title: James-Stein shrinkage 'in the wild'?

text: I am taken by the idea of James-Stein shrinkage (i.e. that a nonlinear function of a single observation of a vector of possibly independent normals can be a better estimator of the means of the random variables, where 'better' is measured by squared error). However, I have never seen it in applied work. Clearly I am not well enough read. Are there any classic examples of where James-Stein has improved estimation in an applied setting? If not, is this kind of shrinkage just an intellectual curiosity?


title: Statistically significant vs. independent/dependent

text: What is the difference between having something statistically significant (such as a difference between two samples) and stating if a group of numbers are independent or dependent.


title: Measuring entropy for event sequences?

text: In `TraMineR` the function `seqient()` for calculating the entropy of a sequence object is not available for event sequences (as defined with `seqecreate()`). Does this mean that measures of "variation" does not make logical sense for event sequences? Are there alternatives to measures such as entropy, complexity, and turbulence for event sequences? A workaround would naturally be to use the `TSE_to_STS()` function in `TraMineRextras` to convert the event sequence to an STS sequence object, after which functions such as `seqient()` becomes available. However, would the statistics coming out of this be reliable?


title: Generalize the usage of moments in method of moments?

text: In Method of Moments for estimation, if there are $k$ parameters to estimate, we usually consider $i$-th moments, $i=1,...,k$, so that we have k equations for k unknowns.   1. I wonder if it is wise to consider more moments of different orders, i.e. $i$-th moments, $i=1,...,n>k$, so that we are to solve a over-determined linear system? Why?   2. Also would it be better if we choose moments of other orders, instead of $i=1,...,k$? Thanks and regards!


title: Significant box's test equality covariance matrices, How bad it is?

text: Having equal number of cases across two groups (60 total cases), if one wants to perform a split-plot Anova and finds out that the Box's Test of Equality of Covariance Metrices is less than 0.001 then what does it mean? I've read that it means that I have violated the assumption of multivariate normal distribution?!! To be honest, I have no idea. I would be grateful to know that by using Pillai's trace rather Wilks's Lambda would the problem be solved? That's all? Any threat to the validity of my research results? (The Pillai's trace returns sig < 0.05.)


title: Differences on exploratory factor analysis, confirmatory factor analysis and principal component analysis

text: Before it is pointed, I am aware that a very similar question was already asked. Still, I am in doubt regarding the concept. More specifically, it is mentioned by the most voted answer that: > In terms of a **simple rule of thumb** , I'd suggest that you: > >   1. Run factor analysis if you assume or wish to test a theoretical model > of latent factors causing observed variables. > >   2. Run principal components analysis If you want to simply reduce your > correlated observed variables to a smaller set of important independent > composite variables. > > **Question 1:** I am having difficulties on understanding based on the results I obtained from R where exactly I am inputing my _theoretical model of latent factors_. I am using the functions from statsmethods. On both **factanal()** and **princomp()** the inputs were the same: A table where each row represented one data point and the columns consisted of different attributes I was interested on reducing. Thus, this add to my confusion on where is this pre assumed model play its role. I noticed that for factor analysis function I used parallel analysis also suggested by the site using the nScree() function to determine the number of factors and I specified if I wanted a varimax (orthogonal) or promax (oblique) rotation. Is that what is it mean by the model? Being able to choose the amount of factors and the type of rotation? The results being provided as visual graphs for both PCA and EFA also doesn't seem to highlight this difference which adds to my confusion. Where does this distinction can be observed on them? ![PCA](http://i.stack.imgur.com/SGK56.jpg) PCA ![EFA](http://i.stack.imgur.com/Eg4MN.jpg) EFA **Question 2:** \-- Answered I bought a book to study about this from Richard L. Gorsuch. On this book there is something that the author caught attention on the difference between PCA (Principal Component Analysis) and EFA (Exploratory Factor Analysis): It is mentioned that PCA is for **population** while EFA is for **sample**. Is that true? I didn't see that being mentioned on any discussion I read so far. Is it irrelevant? **Question 3:** I noticed that all those methods seems to impose the normal distribution constraint. I also read that for larger sets this constraint can be ignored. Is that true or PCA, EFA and CFA are sensible to distribution constraint violations? **Question 4:** Where from the results of PCA and EFA should I note that one is talking about latent factors (EFA) and the other is just clustering on components (factors) the variables? The outputs from R looks the same to me. Is it just the way I perceive what the factors being shown as output? I noted that both show me the table where I can see which I can observe which of my variables are expressed the most of my factors. **What is the difference on the interpretation I should have on which variable belongs to which factor** in respect to PCA and EFA? EFA is saying those with higher expression seems to be more explained by that latent factor while PCA is trying to say that factor is holding those variables from what is it observed? **Question 5** Finally the last question is regarding CFA (Confirmatory Factor Analysis). On the same function website the following image is being shown: ![Confirmatory Factor Analysis](http://i.stack.imgur.com/NNrvh.gif) I read that CFA is usually followed after EFA for hypothesis testing. In that sense, EFA tells you which are the latent factors (which are the output factors) and then you use CFA assuming those factors you observed from EFA for hypothesis testing? **Question 6** For EFA one of the available rotations on the literature is direct oblimium. I heard that it can accounts for both promax and varimax so 'it takes the best of two words'. Is that true? I am also trying to find a function that employs them on R, since the one suggested on the site does not. I would be happy to get any suggestion on this one. * * * I hope it is noted that this question is way more specific on the doubts regarding EFA and PCA and also adds to CFA so not to get closed for being repeated on the subject. If at least one of the questions is answered I am more than happy too as to clarify the confusion in my head. Thank you.


title: How many people need to be screened if the disease has an R0 of 2, and starts with 1 infection in 1000?

text: I need to calculate how many people I would need to screen to catch a disease in its early stages. If the disease has an R0 of 2, and starts with 1 infection in 1000, and I screen one person in every 1000 - how many people will the disease have reached by the time I screen an infected individual? How can I do this on excel?


title: Is the expected value a valid decision-making method in a very short term?

text: This might be related to game theory more than statistics, but I decided to ask this question here. Let's assume you're offered a lottery. There are a hundred balls in a bowl: 99 white balls and one red ball. You are allowed to pick one ball at random (obviously without seeing what color it is). If the ball you picked is red, you win \$100. However, if you pick a white ball, you lose \$1. The expected value of the lottery is .99* (-\$1) + .01 * \$100 = \$.01. This suggests that you should play the lottery, since it allows you to win one cent on average. The expected value in this case means that if you play the lottery an infinite amount of times, you'll average a gain of one cent per game. However, what if you're only allowed to play the lottery once or twice or three times? Is the expected value still a viable decision-making device in such a scenario? Thank you for your responses.


title: Interpretation of logistic regression intercept with one dummy coded categorical variable

text: Feel free to critique my overall approach instead of answering my question directly. I want to look at bivariate relationships among a binary outcome and multiple predictor variables before conducting multiple regression. EDIT: [The data are multiply imputed] and the categorical predictors have been dummy coded. Therefore, for each nominal categorical variable in the original data, there are k-1 dummy variables, with the omitted category serving as the reference group. It seemed to me that there was no meaningful "bivariate" relationship between the outcome and a single [dummy] variable from a collection of related variables (this would change the comparison category to "everything else," which could not be the case in a multiple regression). For this reason, I did not use bivariate correlations. Instead, I did a series of "bivariate" logistic regressions including in a single regression all k-1 dummies of a variable. So, like: outcome = race_black race_other (race_white omitted). In any case, the reference category cannot be examined directly since it was omitted in the multiple imputation process and has missing values. I figured if any coefficient of a binary dummy were significant, I would include the whole variable (group of k-1 dummies) in the omnibus regression. After laboriously compiling a table of the results of these "bivariate" regressions without including the intercept, it occurred to me that (MAYBE?) the intercept is the "dummy" for the omitted category. If it is significant, is it in fact indicating that the omitted group is different from the mean on the outcome? If that's true, should I include the variable in my omnibus regression even if none of the explicit category variables are significant predictors? I'm worried that I am thinking wrongly about the comparison groups and meanings of significant coefficients.


title: Multiple probabilities to single

text: This is a pretty simple question. If I win 7 dollars with a probability of 0.33, and 3 dollars with a probability of 0.35, can I substitue this multiscenario event to a single event of winning 7(33/68) + 3(35/68) with a probability of 0.68. My goal is, given a money game, to quickly compare which side is more advantageous (without having to calculate the expected value of the game). When it is a simple problem, it's all fine. Like:      1\. I receive 30 dollars if $x$ happens and I pay 10 dollars if $x$ doesn't happen, I can easily compare 30/40 with the probability of $x$ not happening, if it's bigger, then it's worth it, else it isn't. But multiscenario problems are hard for me to deal with. For example:      2\. I receive 25 dollars if $y$ happens and 12 dollars if $z$ happens, and I pay 32 dollars if $k$ happens and 19 dollars if $p$ happens , where $y,z,k$, and $p$ are independent events, with $P(y) + P(z) = 1 - P(k) - P(p)$. What is a quick way (without evaluating the expected value) to tell if it's worth it?   I was thinking of reducing each of the 2 multiscenario events in 2 single- scenario events, and deal like in the case of 1.