title: How to find classifier success performance?

text: **EDIT** I am playing with pattern recognition techniques and just to get a grip of it for simplicity I have tried to develop a classifier which categorizes strings into 3 classes based on the probabilistic frequency count under labels A,B,C; each label indicating the dominance of the letters which means if a string of fixd length L=8 contains maximum A then it is classified under A and so on. SO,these numerical results form the features. Is there a code which plots and gives the success rate of classifiers when using the same sample and training set which are data points. The data points are the features which indicate the frequency count of letters in a string. I am interested to use k-NN, Bayes Classifier and Piecewise Component Analysis (PCA). I am aware of the cross validate function but unable to use it for the purpose. There are 3 classes each containing 100 rows of single column data and an unknown sample of same size.   1. The issue is just to confirm that would all the other generic classes also have this same probabilistic numbers as their features? If so, then how to work with k-NN since it computes the distance between the strings (in htis case I guess it would be simply asscii values or euclidiean) and my classifier computes the frequency count.   2. What would be the features for k-NN,Bayes,and PCA?    3. For k-NN would k=number of classes -1 ?   4. I have plotted the ROC for my classsifier and as the design goes that a string would certainly be classified under one of the classes.So,each curve of the ROC indicates a class,hence there are 3 such curves.Is this approach ok?   5. How to proceed with a comparative ROC for all the 4 classifiers?


title: Fitting all statistical distributions in R and choose best fit

text: Is there any function in R that fit all statistical distributions, and choose best fit based on log-likelihood and Kolmogorov-Smirnov D (KSD) statistic? Like a software named, EasyFit


title: Regression $R^2$ and correlations

text: I understand that for simple linear regression, the sample correlation coefficient is the square root of the $R^2$. But that's just for a simple (i.e., single variable) regression $Y=\beta_0+\beta_1X+\varepsilon$. How about multiple regression, e.g., $Y=\beta_0+\beta_1X_1 + \beta_2X_2+\varepsilon$? Is there any relationship between the correlations $corr(Y, X_1)$, $corr(Y, X_2)$ and the regression $R^2$? Thanks!


title: Difference between superpopulation and infinite population

text: What is the difference between Superpopulation and Infinite population? Please explain this with examples. Thanks!


title: If $X \sim \mathrm{Gamma}(\alpha,1)$ and $Y \sim \mathrm{Gamma}(\beta,1)$, is $X/(X+Y) \sim \mathrm{Beta}(\alpha,\beta)$?

text: According to Wikipedia : $\mathrm{Beta}(\alpha,\beta) = \mathrm{Gamma}(\alpha,\theta) / (\mathrm{Gamma}(\alpha,\theta) + \mathrm{Gamma}(\beta,\theta))$ However, when I try to simulate this in R :               > m <- 2^11     > a <- rgamma(m,1,1) / (rgamma(m,1,1) + rgamma(m,1,1))     > quantile(a,1:9/10)     0.05576023 0.12110211 0.19886341 0.29088744 0.41851181      0.56698135 0.80151216 1.20748646 2.13406961     > qbeta(1:9/10,1,1)     0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9      Those two distributions don't look the same on the lower and upper percentiles to me, even when I increase the simulation loop index `m`. Simulation works for higher $\alpha$ and $\beta$ though. I understand than since $\mathrm{Gamma}(1,1)$ is identical to $\mathrm{Exponential}(1)$, dividing one exponential number by the sum of two exponential numbers can produce a result greater than 1, whereas the Beta distribution supports only real numbers on $(0,1)$, so I don't understand how the formula above could be correct for any $\alpha$ and any $\beta$ (even small ones) ? Did I make a confusion between the Gamma function and the Gamma distribution, or is it something else?


title: How to visualize two bar charts with very different scales without looking redundant

text: More an aesthetics question to do with presentation of statistical data, say you have 2 sets of data,                    speed     weight     a      2.2      500     b      4.7      222     c      7.3      999     d      3.1     1000      So if you plot the speeds and weights on the same bar chart, you will get squishy tiny bars for all the speed values (making them unreadable). The data should be viewed together though, because there is some relationship between speed and weight that you're trying to show. What I came up with was adjacent bar charts: ![enter image description here](http://i.stack.imgur.com/5wiME.png) But it kind of looks redundant, with the repeated axes there.


title: How to measure the number of people in a picture of a crowd?

text: **Background** : Israel (and the middle east in general) is filled with protests. I am curious, when given a picture, to estimate how many people are in that picture (often a picture of a large crowd). What modelling can offer some solution to this problem? (and, of course, can it be done with any open source package. Say, R?)


title: Confidence interval for a regression parameter via prediction

text: Consider a simple Poisson-regression - GLM - model. There $\exp\left(\beta\right)$s are used as Incidence Rate Ratios (IRR), but their calculation is sometimes not completely straightforward, for instance, if we have interactions (when one has to take the product of such terms). It is especially true if we need a confidence interval for an IRR: it is easy to obtain a CI for a $\beta$, but in that case we need a CI for the _exponentiated_ of the _sum_ of such parameters. Motivated by this, I had an idea that might be completely wrong, but I thought it worth asking. What if I calculate IRRs and their confidence intervals by using predictions? Let me show an example (in R). Take the well-known dataset `warpbreaks`:                head( warpbreaks )      mod <- glm( breaks ~ wool * tension, warpbreaks, family = poisson )      To obtain IRRs:                exp( coefficients( mod ) )        (Intercept)          woolB       tensionM       tensionH woolB:tensionM woolB:tensionH          44.5555556      0.6334165      0.5386534      0.5511222      1.8930264      1.2072719       It can be seen what I meant by ''not straightforward'': to obtain IRRs for each covariate combination we need               irrmod <- c( 1, exp( coef( mod ) )[ 2:3 ],     exp( coef( mod )[ 2 ] + coef( mod )[ 3 ] + coef( mod )[ 5 ] ), exp( coef( mod )[ 4 ] ),                  exp( coef( mod )[ 2 ] + coef( mod )[ 4 ] + coef( mod )[ 6 ] ) )     cbind( irrmod )                 irrmod              1.0000000     woolB    0.6334165     tensionM 0.5386534     woolB    0.6458853     tensionH 0.5511222     woolB    0.4214464      And we haven't even talked about CIs for these... So, here is the idea: first, we obtain predictions for each possible combination of explanatory variables.               preddata <- expand.grid( wool = levels( warpbreaks$wool ),         tension = levels( warpbreaks$tension ) )     preds <- predict( mod, newdata = preddata, type = "response", se.fit = TRUE )     irrpred <- preds$fit            1        2        3        4        5        6      44.55556 28.22222 24.00000 28.77778 24.55556 18.77778      The problem is that these are measured on the response scale, and the (possible?) solution: divide all of them with the first! As these are all categorical variables _I think_ this results in IRRs. Indeed, in the concrete example:               expbeta0 <- irrpred[ 1 ]     irrpred <- irrpred / expbeta0     cbind( preddata, irrpred, irrmod )       wool tension   irrpred    irrmod     1    A       L 1.0000000 1.0000000     2    B       L 0.6334165 0.6334165     3    A       M 0.5386534 0.5386534     4    B       M 0.6458853 0.6458853     5    A       H 0.5511222 0.5511222     6    B       H 0.4214464 0.4214464      so these are in fact equal. My **first question** : is this method theoretically sound, i.e. is it mathematically guaranteed that the two approaches are equivalent in every situation? I'm quite sure they are, just to be safe. But now comes the part in which I'm absolutely not sure: what about CIs? Taking this approach further, it is logical to calculate the following:               upr <- preds$fit + ( qnorm( 0.975 ) * preds$se.fit )     lwr <- preds$fit - ( qnorm( 0.975 ) * preds$se.fit )     upr <- upr / expbeta0     lwr <- lwr / expbeta0      I have to admit that I'm not too familiar with how `predict.glm` calculates standard errors, so I might be completely wrong, but can this method be used to obtain CIs for the IRRs...? An obvious problem is that this way we get a ''true CI'' even for the reference category (which, of course, should not have a CI, being fixed at $1$). The other problem is that the results are not equal to that of `confint` even in the categories where they are comparable:               cbind( preddata, lwr, exp( confint( mod ) )[ , 1 ], upr, exp( confint( mod )[ , 2 ] ) )     Waiting for profiling to be done...     Waiting for profiling to be done...       wool tension       lwr exp(confint(mod))[, 1]       upr exp(confint(mod)[, 2])     1    A       L 0.9021242             40.3356421 1.0978758             49.0597741     2    B       L 0.5555195              0.5407439 0.7113134              0.7405983     3    A       M 0.4668193              0.4558586 0.6104874              0.6347433     4    B       M 0.5672254              0.4670131 0.7245452              0.6486862     5    A       H 0.4784616              1.4908755 0.6237828              2.4069226     6    B       H 0.3579064              0.9356292 0.4849864              1.5571553      Note that only the second and third line is interesting now, as only they pertain to the same combination. (And the fifth row from `upr/low` vs. the fourth row of `confint`.) While they are _quite close_ , they do not match exactly... So my **second question** is: what about the CIs for the IRRs? Can this method, with `predict`, be used to calculate them?


title: SVM parameter selection with NM simplex (or other algorithms)

text: I'm having some trouble getting the NM Simplex to find a good minimum for selecting hyperparameters of a rbf SVC. Not only am I tuning the 2 SVC parameters (C and gamma) I also have five class weights that I'm trying to tune. Grid searches seem fairly infeasible in this case. The NM Simplex doesn't seem test inputs far away enough from the initial guess. I am transforming the inputs (e.g. the NM Simplex is searching in range i while 2^i is inputted as the C parameter for the SVC). Even so, the NM Simplex always seems to end up back at the initial guess. Any suggestions or any other minimization algorithms that might be worthwhile? Thanks


title: Rank ordering and/or classification problem

text: I have about 10 variables about products. I need to rank order the products by quality. Now, the quality is very subjective thing, so I got all 300 products ranked by domain experts. Now I have 10 variables, and ranking assigned by an expert. About 1% of products should be excelent, about 10% should be very good, 30% should be good, 30% should be average and the rest is below average. What sort of modeling technique would you use to make sue of this data and create system which will automatically be able to classigy product into those categories and mimic "expert judgement". Thanks


title: Standard error of parameter estimates in regularized regression

text: In a regularized linear regression model (e.g., ridge regression, lasso, etc.), what is the best way to obtain standard errors for parameter estimates? If cross-validation is used, is it statistically sound to estimate the standard error for each parameter as the standard deviation across folds?


title: Probabilistic model for individual income change

text: For an agent-based simulation, I need to model the change in the income of an individual over time. I just got an answer about Probability distribution of income , so I can determine the income of individuals at t=0, but how would their income change over time? Is this is a Markov series, and if so, what would be reasonable parameters for that model?


title: Open source data mining package comparsion

text: Is there sources where I can see performance of different ML models applied to same data? Currently I want to understand performance percision/recall difference between mahout online logistic regression vs Mallet MaxEnt, LibLINIEAR.


title: Logistic regression model does not converge

text: I've got some data about airline flights (in a data frame called `flights`) and I would like to see if the flight time has any effect on the probability of a significantly delayed arrival (meaning 10 or more minutes). I figured I'd use logistic regression, with the flight time as the predictor and whether or not each flight was significantly delayed (a bunch of Bernoullis) as the response. I used the following code...               flights$BigDelay <- flights$ArrDelay >= 10     delay.model <- glm(BigDelay ~ ArrDelay, data=flights, family=binomial(link="logit"))     summary(delay.model)      ...but got the following output.               > flights$BigDelay <- flights$ArrDelay >= 10     > delay.model <- glm(BigDelay ~ ArrDelay, data=flights, family=binomial(link="logit"))     Warning messages:     1: In glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  :       algorithm did not converge     2: In glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  :       fitted probabilities numerically 0 or 1 occurred     > summary(delay.model)          Call:     glm(formula = BigDelay ~ ArrDelay, family = binomial(link = "logit"),         data = flights)          Deviance Residuals:            Min          1Q      Median          3Q         Max     -3.843e-04  -2.107e-08  -2.107e-08   2.107e-08   3.814e-04          Coefficients:                 Estimate Std. Error z value Pr(>|z|)     (Intercept)  -312.14     170.26  -1.833   0.0668 .     ArrDelay       32.86      17.92   1.833   0.0668 .     ---     Signif. codes:  0 â***â 0.001 â**â 0.01 â*â 0.05 â.â 0.1 â â 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 2.8375e+06  on 2291292  degrees of freedom     Residual deviance: 9.1675e-03  on 2291291  degrees of freedom     AIC: 4.0092          Number of Fisher Scoring iterations: 25      What does it mean that the algorithm did not converge? I thought it be because the `BigDelay` values were `TRUE` and `FALSE` instead of `0` and `1`, but I got the same error after I converted everything. Any ideas?


title: How to make Random Forests more interpretable?

text: Are there any methods that one could utilize to make Random Forest more interpretable? Random Forest performs much better than CART but it is a lot less interpretable.


title: How to sample from the prior predictive distribution in jags?

text: Is there an example available of how I can sample from the prior predictive distribution (without data) in jags? I would like to get a better sense for the contribution of the prior in a multilevel model using jags. I could probably set up a simulation manually using R, but since I already have the distributional relationships set up in jags, it would be preferable if I could do it without setting up duplicate simulation code in R.


title: Normalizing depth data

text: I have some Kinect data of somebody standing (reasonably) still and performing sets of punches. I am given it in the format of an x,y,z co-ordinate for each joint of which they are 20, so I have 60 data points per frame. I'm trying to perform a classification task on the punches however I'm having some problems normalising my data. As you can see from the graph there are sections with much higher 'amplitude' than the others, my belief is that this is due to how close that person was to the kinect sensor when the readings were taken. (The graph is actually the first principal coefficient obtained by PCA for each frame, multiple sequences of the same punch are strung together in this graph.) Looking back at the data files it looks like those that are 'out' have a z co- ordinate (depth from sensor) of ~2.7m where as the others tent to hover around 3.3-3.6m. **How can I perform a normalization with the depth values to make them closer to each other for each sequence?** I've already tried differentiation to get the velocity, although it helps to normalise the output actually ends up too similar and makes it very hard to classify. ![enter image description here](http://i.stack.imgur.com/85mgV.jpg) Edit: I should mention I am already using a normalization method by subtracting the hip position from each joint in an attempt to make the co- ordinates relative.


title: Matrix singularity issues with gaussian variogram

text: I have defined my gaussian variogram like this $r(h) = \text{nugget} + \text{partial_sill}\cdot(1 - \exp(-\frac{3h^2}{\text{range}^2}))$ I set nugget = 0.1343 partial_sill = 0.3125 range = 19.8642 I tried to create a variogram for a spatial grid of size 5x5, meaning for each of the location in the 5x5 grid, I calculated the value of the variogram with distance equal to its distance to the other locations in the grid. So I had a matrix of size 25x25. However, this matrix is not positive definite and is singular. How can I remedy this issue with gaussian variogram?


title: What is randomness?

text: In probability and statistics, the concept of "random" and "randomness" are frequently used. Often the concept of a random variable is used to model events that occur due to chance. My question regards the term "random". What is random? Does randomness really exist? I am curious what people that have a lot of experience in working with random events think and believe about randomness.


title: MLE of the mixture parameter in mixing two normal densities

text: Imagine that we have mixture of two normal distributions with mixture parameter $\theta$: $$p(y_i|\theta) = \theta\phi(y_i;\mu_1, \sigma_1^2) + (1 - \theta)\phi(y_i; \mu_2, \sigma_2^2)$$ Assume that the ONLY unknown parameter is $\theta$. I want to get the MLE estimate of $\theta$. Here is what I do: $$L(\theta | data, \mu_1, \sigma_1^2, \mu_2, \sigma_2^2) = \prod_{i = 1}^{n} p(y_i|\theta) = \prod_{i = 1}^{n} \Big(\theta\phi(y_i;\mu_1, \sigma_1^2) + (1 - \theta)\phi(y_i; \mu_2, \sigma_2^2)\Big) \Rightarrow$$ $$L(\theta | data, \mu_1, \sigma_1^2, \mu_2, \sigma_2^2) \propto \prod_{i = 1}^{n}\Big[ \theta \Big(\phi_{i1} - \phi_{i2}\Big)\Big] = \theta^n \prod_{i = 1}^{n}\Big(\phi_{i1} - \phi_{i2}\Big) = \theta^nG(Y, param)$$ where $\phi_{i1}$ refers to the first normal density with $y_i$ plugged in and $\phi_{i2}$ refers to the second normal density. Note that there is no $\theta$ involved in G(Y, param). Now, in order to get the MLE, I take log and then derivative: $$l(\theta | data, \mu_1, \sigma_1^2, \mu_2, \sigma_2^2) = log(L(\theta | data, \mu_1, \sigma_1^2, \mu_2, \sigma_2^2)) \propto n log(\theta) + log(G) \Rightarrow $$ $$l'(\theta | data, \mu_1, \sigma_1^2, \mu_2, \sigma_2^2) = \frac{n}{\theta}$$ If my calculations above is right, then it's a very weird derivative and I don't know how to get the MLE estimate. I appreciate if you could guide me on that.


title: Power of a Multiple Linear Regression

text: I have a set of models that are the result of a multiple linear regression. I would like to calculate the power for each of these models. I found this tutorial on calculating the power using G*Power. Since the data has already been collected, I am using the Post Hoc type of Power Analysis, and I was hoping to use the "Determine" button to automatically calculate the effect size. However, I see that one need to provide a **partial** R2, but I would think that for a multiple linear regression (which I think would have multiple partial R2 values), it would make more sense to provide the R2 for the resulting model. Am I misunderstanding the difference between R2 and partial R2? If so, how do I obtain the partial R2 to use in G*Power?


title: The sample size applied to a non-normal distribution

text: I have a single variable that represents my population values (sample of data):               [1]  94.51  59.81  63.84  94.51  94.51  94.51  94.51  94.51  94.51  94.51     [11]  59.81  94.51  94.51  94.51  47.90  29.16  50.36  23.51  44.41  33.14     [21]  47.90  29.16  47.90  29.16  47.90  29.16  47.90  29.16  47.90  29.16     ...     [331]  23.44  24.52  12.37  29.12  24.52  12.37  29.12  24.52  12.37  29.12     [341]  24.52  12.37  29.12  24.52  12.37  29.12  24.52  12.37  29.12  24.52     [351]  12.37  29.12  24.52  12.37  45.25  25.78  49.84  29.12  24.52  12.37     [361]  29.12  24.52  12.37  29.12  24.52  12.37               > summary(group$V1)         Min. 1st Qu.  Median    Mean 3rd Qu.    Max.          6.11   35.94   59.13   62.31   86.10  111.50      > mean(group$V1)      [1] 62.30546     > sd(group$V1)      [1] 29.55491      The corresponding histogram is: ![Histogram of BitScores of the Population](http://i.stack.imgur.com/Cp6dW.png) And the Shapiro test of normality:               Shapiro-Wilk normality test                   data:  group$V1         W = 0.9466, p-value = 3.161e-10      With the last information my conclusion is that the population is not distributed normally. The objetive is extract a sample from these population, but I have problems to apply a method to determine the sample size, because in some methods the assumption is based on the normality of population. (According with these reference) The sample is required to comparate this group with a random group with the same sample size, and the single variable to evaluate is the Bitscore. Some references, suggestions, approaches? Thanks in advance.


title: Feature Selection Packages in R

text: I am very new to R. I am learning machine learning right now. Very sorry, if this question appears to be very basic. I am trying to find a good feature selection package in R. I went through Boruta package. It is a good package but I read that it is only useful for classification. I want to do implement feature selection in R for regression tasks. I went through the caret package documentation but for my level, it is very difficult to understand. Can any one please point me to a good tutorial or list any good packages or most frequently used packages in R for feature selection. Any help would be appreciated. Thanks in advance.


title: Determining if time has influence on correlation between x and y

text: Is there a sort of statistics model I could use to determine if time has an influence on the correlation between two sets of data? For example say I have a column $X$ and a column $Y$, both just columns of numbers, and then also a time column. Thus each row has time, $X$ and $Y$. Is there some sort of statistical analysis I could do that tells me if time has an influence on the relationship between $X$ and $Y$? I ran across the Pearson Correlation Coefficient to measure correlation between two variables. I'm completely guessing here, but for example could I repeatedly use the Pearson Correlation for different segments of time, and then simply determine if there is a change in the Pearson Correlation Coefficients? Is there a better way?


title: Is it correct to say the "smooth curve fitted by LOESS"?

text: I use the package `scatter.smooth` in R, and I want to have a figure text describing the line in the scatterplot for my report. Is is okay to say "the line is a Smooth Curve Fitted by LOESS ( _a reference to the package_ )" or need I explain it more? (It's a software engineering master thesis) Otherwise, how can I describe what the line is in one short sentence?


title: Would a Bayesian admit that there is one fixed parameter value?

text: In Bayesian data analysis, parameters are treated as random variables. This stems from the Bayesian subjective conceptualization of probability. But do Bayesians theoretically acknowledge that there is one true fixed parameter value out in the 'real world?' It seems like the obvious answer is 'yes', because then trying to estimate the parameter would almost be nonsensical. An academic citation for this answer would be greatly appreciated.


title: How to calculate the prediction interval for my 24-hour forecast?

text: I'm doing time series analysis and I would like to calculate the prediction interval. I have predicted the next 24 values for my time series, but now I would also like to calculate the prediction interval for them (say 95% interval). How to do this? Can someone show me the steps? I found a pdf on the subject, but it only showed how to make the prediction interval for one step- forecast...here is the pdf if you want to see... http://www.learningace.com/doc/2087382/d6c17c68ca3c0a03278d4a97276a0223/forecast- conf-int-11 Thank you for any help =) Update to question: When answering my question, could you please explain what to do and not redirect me to a R-package or something. I want to know what I'm doing =) Thank you


title: What happens to the covariance matrix when the errors are independent?

text: I was wondering, what happens to the covariance matrix of the errors, when I assume that all the errors are stochastically independent? Is the covariance matrix still: $$\sigma^2I = \begin{bmatrix}\sigma^2 & ... & 0\\\\. & . &.\\\ 0 & ... & \sigma^2\end{bmatrix}$$ I really appreciate your answer!!!


title: In a "bursty" dataset, how do you filter for the few important values that make up the bulk of the information?

text: Note sure if there is an existing stats concept for this but I have a dataset that consists of mostly small data points with a few large ones. e.g. 1 2 1 3 1 2 87 3 2 1 1 1 1 3 1 2 1 1 1 99 How can I filter this data set to only get the values that disproportionately make up the bulk of the information? I am currently filtering by data points that exist a few standard deviations out but this doesn't tell me what % of the total I am getting. (e.g. if I go 2 standard deviations out am I getting 70% of the info? if I go 5 is it 95%? I only know what % of the number of data points it represents not the percentage of the data) EDIT: I want to remove as many data points as possible without removing the important data points. So if I have a mean of 5 and a std-deviation of 20, I filter out data points that are less than 45 (20 + 20 + 5). This removes say 95% of the data points but then the dataset can look like: 50 46 90 80 44 99999 57 87 88. The Pareto Principle here is applied recursively with this 99999. In this scenario I'd like to only keep the 99999 since it accounts for 99% of the data but I don't know this from just using a std-deviation rule of thumb. For example, many people will agree that 1% of people can hold 99% of the wealth. If you slice into that data further you find that 1% of that 1% holds 99% of that wealth. Meaning that 0.01% of people hold 98% of the wealth. This second piece of information is surprising since it shows the "big guys" of the "big guys". This might even go further with the "big guys" of the "big guys" of the "big guys" (big-guys^3) Maybe 1 person holds 95% of all the wealth. How can I analyze my data for this? Given a pie or bar-chart it would be obviously at a glance.


title: Comparison of p values for Levene mean test and Levene median test?

text: I am doing Levene's mean test and Levene's median test (Brown-Forsythe). I want to compare the p-values of these two tests to see which is better. I get large p-values for both tests which are 0.562 (Levene mean) and 0.611 (Levene median) for normal distribution.   * Which test shows the better type I error rate?    * does Levene's mean test perform best when the data follows a normal distribution?


title: Normalization of categorical factor variables

text: For supervised machine learning for prediction, if I had some feature variables that are real, and also some features that are categorical--which have been coded using dummy variables (010, 001 etc)--I have normalized the real variables so that each of them sums to one. I am wondering what kind of a preprocessing I should do for the rest of the categorical features, before I run cross-validation routines and regression methods.


title: Parameter estimation for linear system with correlated noise

text: Let's say that I have a linear system (e.g. an electrical circuit) and I am trying to estimate the values of parameters in the system (e.g. the resistance, inductances, and capacitances). I do this by measuring the response of the system to some input signal (e.g. a step function). If the noise in the measurement were uncorrelated and Gaussian, then one obvious approach to parameter estimation would be to use the maximum-likelihood estimator, which would reduce to a (non-linear) least square fit. Simple enough. However, if the measurement is very sensitive, so that the dominant source of noise comes from compoments internal to the system (e.g. Johnson noise in a resistor), the noise at the output will be correlated. This means that maximum likelihood estimator for the parameters will not reduce to a simple least squares procedure. I'm sure that there is a set of statistical practices and literature devoted to parameter estimation for systems like this, but I'm not sure where to start looking. Can anyone point me in the right direction? **Update** I just want to clarify that although the system is linear in the sense that a linear superposition of two signals is also a valid signal, the output of the system as a function of time is certainly _not_ linear. In general it would be some sum of exponentials, where the coefficients in front of each exponential are function of the various time constants.


title: Building a time series that includes multiple observations for each date

text: I'm trying to apply a time series to quarterly sampled data (animal biomass) over a 10 year period with 3 reps per quarter. So 40 dates but 120 total observations. I have read up to SARIMA'a in Shumway and Stoffer's Time Series Analysis and it's Applications as well as skimmed Woodward, et. al.'s Applied Time Series Analysis, and my understanding is each model is based on a single observation at each point in the time series. QUESTION: How can I include the variation in each observation in my model? I could build a series on the mean, but I would loose the variation at each observation and I think that is critical to my understanding of what is happening.


title: Choosing optimal alpha in elastic net logistic regression

text: I'm performing an elastic net logistic regression on a health care dataset using the glmnet package in R by selecting lambda values over a grid of alphas from 0 to 1. My abbreviated code is below:               > alphaslist<-seq(0,1,by=0.1)     > elasticnet<-lapply(alphalist, function(a){cv.glmnet(x, y, alpha=a, family="binomial", lambda.min.ratio=.001)})     > for (i in 1:11) {print(min(elasticnet[[i]]$cvm))}      which outputs the mean cross validated error for each value of alpha from 0.0 to 1.0 by 0.1:               [1] 0.2080167     [1] 0.1947478     [1] 0.1949832     [1] 0.1946211     [1] 0.1947906     [1] 0.1953286     [1] 0.194827     [1] 0.1944735     [1] 0.1942612     [1] 0.1944079     [1] 0.1948874      Based on what I've read in the literature, the optimal choice of alpha is where the cv error is minimized. But there's a lot of variation in the errors over the range of alphas. I'm seeing several local minimums, with a global minimum error of 0.1942612 for alpha=0.8. Is it safe to go with alpha=0.8? Or given the variation, should I rerun cv.glmnet with more cross validation folds (say 20 instead of 10) or perhaps more increments between alpha=0.0 and 1.0 to get a clear picture of the cv error path?


title: Probability that someone will like image

text: I've got following problem:   \- We have set of N people   \- We have set of K images   \- Each person rates some number of images. A person might like or not like an image (these are the only two possibilites) . \- The problem is how to calculate likelihood that some person likes a particular image. I'll give example presenting my intuition.   N = 4   K = 5   \+ means that person likes image   \- means that person doesn't like image   0 means that person hasn't been asked about the image, and that value should be predicted               x 1 2 3 4 5         1 + - 0 0 +        2 + - + 0 +       3 - - + + 0       4 - 0 - - -      Person 1 will probably like image 3 because, person 2 has similar preferences and person 2 likes image 3.   Person 4 will probably not like image 2 because no one else likes it and in addition person 4 does not like most images. Is there any well known method, which can be used to calculate such likelihood?


title: Is it okay to take the average of several normalized time series?

text: I have 3 time series (10 years) of economic data for 3 countries. For another analysis that I am doing, I should try and reduce these into preferably 1 time series. For this analysis I don't care so much about the absolute values, rather about the variation. The three time series are for the same thing, but the magnitudes vary between the countries, so obviously I need to normalize. But should I normalize all three time series and take the average, or should I take the average and normalize this? I'm leaning towards the first option (first normalize the 3 time series and then take the average) but I am not sure this is statistically correct or okay…?


title: How do you draw structural equation/MPLUS models?

text: I am looking for a software tool (preferably open source) to draw structural equation/mixture models efficiently and prettily. After looking into xfig and graphviz I now stick to the general vector graphics package inkscape because it seems most flexible. I would like to poll the stat.stackexchange community: How do you draw your structural equation/mixture models? What software do you use?


title: Geometric interpretation of the maximum likelihood estimation

text: I was reading the book _The Identification Problem In Econometrics_ by Franklin M. Fisher, and was confused by the part that he demonstrates the identification by visualizing the likelihood function. The problem could be simplified as: > For a regression $Y=a+Xb+u$, where $u \sim i.i.d. N(0,\sigma^2I)$, $a$ and > $b$ are the parameters. Suppose $Y$ has a coefficient $c$ which equals to > unity. Then the likelihood function in the space of $c, a,b$ would have **a > ridge along the ray corresponding to the vector of true parameters and its > scalar multiples**. When considering only the place given by $c=1$, the > likelihood function would have a unique maximum at the point where the ray > intersected that plane. My questions are:   1. How should one understand and reason about the ridge and the ray mentioned in the demonstration.   2. Since the ray are the true parameters and scalars, why is the ray not on the plane given by $c=1$ since the true value of parameter $c$ is 1.


title: Sample from continuous uniform distribution with open interval

text: If I want to sample from a continuous uniform distribution with interval $(a,b]$, how can I do it in `R`? Or is it just the same as sampling from $(a,b)$ in `R` using               runif(1,min = a, max = b)      ?


title: Is there a term for generating time-line based data from individual points? Also how would I do this?

text: Is there a term that describes what I'm trying to do below? Also, how would you do this using something like JMP or Excel? (or do I need to code this in something like perl?) Given this sort of data:               ID| opened     | closed     | quantity     --------------------------------------     1 | 2010-01-01 | 2010-01-03 | 1     2 | 2010-01-02 |            | 2     3 | 2010-01-02 | 2010-10-05 | 3      I'd like to get this data and then graph with x being a time line and Y being total quantity open:               on date    | total quantity open     ---------------------------------     2010-01-01 | 1     2010-01-02 | 6     2010-01-03 | 5     2010-01-04 | 5     2010-01-05 | 2


title: How to interpret a control chart containing a majority of zero values? (Part 2)

text: (I know it is not quite appropriate to quote it as "part 2", but since the question has been dormant for quite a while, I hope by doing this will rise peoples' attention again, you may have a look of part I here.) I have come across an article online talking about the case similar to mine, that most of the time the case count is zero, make sometimes when the case number increases to one, it already shoots above the control level and consider the case as "out of control". Since c-chart will be easier for my bosses to read and interpret, I wonder if the method is sound, or did anyone have some more official reference on this method? (I have googled quite a while but I can find nothing) The article can be found here. An to further my question, I want to ask one more thing: for the assumption of c-chart that the case count needs to follow a Poisson distribution, is it applicable to all lambda (i.e. mean of case count)? Thanks again.


title: Lasso fitting by coordinate descent: open-source implementations?

text: What open-source implementations -- in any language -- exist out there that can compute lasso regularisation paths for linear regression by coordinate descent? So far I am aware of:   * glmnet   * scikits.learn Anything else out there?


title: Output of linear SVM model in Matlab using SVM-light

text: I am using SVM-light with Matlab, for linear SVM. I would like to understand the output model, but I cannot find any documentation or help about it. Here is the output:                       sv_num: 639        upper_bound: 547                  b: 1.4023           totwords: 576             totdoc: 2000          loo_error: -1         loo_recall: -1      loo_precision: -1           xa_error: 14.9500          xa_recall: 82.6000       xa_precision: 86.8559            maxdiff: 9.9611e-004         r_delta_sq: 575.0052        r_delta_avg: 23.9792       model_length: 0.4847               loss: 357.5542              vcdim: 136.1071              alpha: [2000x1 double]        lin_weights: []              index: [2000x1 double]             supvec: [639x576 double]        kernel_parm: [1x1 struct]     example_length: 23.9792                  a: [2000x1 double]      So far, I understand there are `sv_num=639` support vectors. However, there are 638 positive index in the `index` vector ; ditto for the positive values of the `alpha` vector. The first value of the `alpha` vector is zero, and the non-zero values are then between 2 and 639. Same remark about the `supvec` vector. So it seems there are indeed 638 support vectors. Actually, I would like to get the hyperplane `W` and the parameter `b` myself to see if I could understand the output model correctly. Since the kernel is linear, I thought of a sum over some coefficients times some support vectors. For information, `a(i)` "should" store the values `alpha_i` times `y_i` for training samples `x_i`. And there are 638 non-zero values. Has anyone used svmlight before? Can you understand the output model? Or is there a documentation anywhere? * * * I manage to get one thing right:               MyCoeff = model.a ./ y;     MyIndex = 1+max(model.index,0); % between 1 and 639      The support vector corresponding to any non-zero `MyCoeff(i)` is:   * `x(i,:)`   * `model.supvec(MyIndex(i),:)` And `MyIndex` is 1 iff `MyCoeff` is zero. And `supvec(1,:)` is the zero vector. * * * Here are four pieces of code I have made:               MyCoeff = model.a ./ y;     W1 = zeros(1, h*w);     for i=1:length(MyCoeff)         W1 = W1 + MyCoeff(i)*x(i,:);     end;          W2 = zeros(1, h*w);     for i=1:model.sv_num         W2 = W2 + model.alpha(i)*model.supvec(i,:);     end;          eq((MyCoeff~=0),(model.index>0))     MesIndices = 1+max(model.index,0); % entre 1 et 639     W3 = zeros(1, h*w);     for i=1:length(MyCoeff)         W3 = W3 + MyCoeff(i)*model.supvec(MesIndices(i),:);     end;          W4 = zeros(1, h*w);     for i=1:length(MyCoeff)         W4 = W4 +  model.alpha(MesIndices(i))*x(i,:);     end;      Here are the outputs: ![http://imgur.com/FErWW.png](http://i.stack.imgur.com/A51VP.png) I have understood the `index` link between the support vectors `supvec` and the matrix `x`. But I have not found the link between `alpha` and `a`. * * * I got it: apart from the difference of indices, `model.a` and `model.alpha` are the same.               for i=1:length(MesIndices)         W1 = W1 + model.a(i)*x(i,:);%*y(i);         W2 = W2 + model.alpha(MesIndices(i))*model.supvec(MesIndices(i),:);%*y(i);         W3 = W3 + model.a(i)*model.supvec(MesIndices(i),:);%*y(i);         W4 = W4 +  model.alpha(MesIndices(i))*x(i,:);%*y(i);     end;      The question now is whether I should multiply by `y(i)`. I believe I should NOT since `model.a` is already `y` times the weight of the corresponding support vector, and $y\in\\{-1,1\\}$. So the good result should be the ORANGE face, which you can get by using one of the four lines above. Too bad I could not find any documentation. :')


title: question on presenting multiclass classification results

text: I'm trying to figure out how to present results (e.g. precision, recall, accuracy, ROC curve) for a class project. I have a multi-class classification problem(specific classification algorithm not really important) of classes **A, B, C and D**. The classification algorithm was **trained** on **A, B and C**. During the test phase, classification algorithm would determine unit test data to be **class A** if **likelihood_of_test_data >= likelihood_threshold_of_A** which is obtained in training phase. Similarly classes B and C are also recognized in test phase. If **likelihood_of_test_data < likelihood threshold of A, B and C**, it would be **recognized** as belonging to **class D**. How can I draw the ROC curve in this case? Do I have to calculate precision and recall for each of the classes of A, B, C and D separately? Or I calculate just the accuracy?


title: How to draw ROC curve with three response variable?

text: I have a toy example reproduced below in which the response variable has three possible classes. I am trying to create an ROC but not sure how to deal with it when there are three classes. Any help will be appreciated. Thanks               library(ipred)     control = rpart.control(maxdepth = 20, minsplit = 20, cp = 0.01, maxsurrogate=2, surrogatestyle = 0, xval=25)     n <- 500; p <- 10     f <- function(x,a,b,d) return( a*(x-b)^2+d )     x1 <- runif(n/2,0,4)     y1 <- f(x1,-1,2,1.7)+runif(n/2,-1,1)     x2 <- runif(n/2,2,6)     y2 <- f(x2,1,4,-1.7)+runif(n/2,-1,1)     y <- c(rep(-1,floor(n/3)),rep(0,ceiling(n/3)), rep(1,ceiling(n/3)))     dat <- data.frame(y=factor(y),x1=c(x1,x2),x2=c(y1,y2), matrix(rnorm(n*(p-2)),ncol=(p-2)))     names(dat)<-c("y",paste("x",1:p,sep=""))     dat          plot(dat$x1,dat$x2,pch=c(1:2)[y], col=c(1,8)[y],           xlab=names(dat)[2],ylab=names(dat)[3])     indtrain<-sample(1:n,300,replace=FALSE)     train<-dat[indtrain,]; dim(train)      test<-dat[setdiff(1:n,indtrain),]; dim(test)      test          mod <- bagging(y~.,  data=train, control=control, coob=TRUE, nbagg=25, keepX = TRUE)     mod     pred<-predict(mod, newdata=test[,-1],type="prob", aggregation= "average"); pred      For two class case, I use to do the following but it is no longer valid for three classes.               yhat <- pred[,2]     y = test[, -1]     plot.roc(y, yhat)


title: Understanding a Quantile-Quantile Plot

text: I must have generated at least 5 Q-Q plots until now when trying to fit my data into a known distribution but I just noticed something that I could not understand. In the figure shown below, from what I've read from the wiki, X-axis is supposed to read "Negative Binomial Theoretical Quantiles" and Y-axis is supposed to read "Data quantiles". Agreed that this makes perfect sense. But when I looked at the figure, the X and Y axis go beyond 100 but how can there be quantiles beyond 100? What do they mean if they exist? Or is this graph produced by the qqplot of R totally different? Can someone help me understand this? The way I was generating this data was using the following script:               library(MASS)          # Define the data     data <- c(67, 81, 93, 65, 18, 44, 31, 103, 64, 19, 27, 57, 63, 25, 22, 150,               31, 58, 93, 6, 86, 43, 17, 9, 78, 23, 75, 28, 37, 23, 108, 14, 137,               69, 58, 81, 62, 25, 54, 57, 65, 72, 17, 22, 170, 95, 38, 33, 34, 68,                38, 117, 28, 17, 19, 25, 24, 15, 103, 31, 33, 77, 38, 8, 48, 32, 48,                26, 63, 16, 70, 87, 31, 36, 31, 38, 91, 117, 16, 40, 7, 26, 15, 89,                67, 7, 39, 33, 58)          # Fit the data to a model     params = fitdistr(data, "Negative Binomial")          #using the answer from params create a set of theoretical values     plot(qnbinom(ppoints(data), size=2.3539444, mu=50.7752809), sort(data))     abline(0,1)      ![alt text](http://i.stack.imgur.com/QeB2Y.png)


title: Can the items measure one underlying construct?

text: For an assignment I am asked to test whether different items (all dichotomous) can actually measure one concept. The items: _Please tell me whether or not you think it should be possible for a pregnant woman to obtain a legal abortion . . . A. If the woman's own health is seriously endangered by the pregnancy? B. If there is a strong chance of serious defect in the baby? C. If she became pregnant as the result of a rape? D. If the family has a very low income and cannot afford any more children? E. If she is not married and does not want to marry the man?_ The results of my **latent-class analysis** are that:   * The model with two and three classes do not have an acceptable model fit: > >     2-class: X² = 492.485; L² = 337.1439; df = 20; p = .000; AIC = 297.1439; > BIC = 195.7409 >     3-class: X² = 170.096; L² = 26.4617; df = 14; p = .0226; AIC = -1.5383; > BIC = -72.5204 >     * There is a clear pattern in the estimated classes. One group agreeing on all items X1, one agreeing on only the medical/ethical reasons (health, defect and rape) , one never agreeing X3 _(exept sometimes for the health item, I think because this is a balance between killing the baby or killing the mother, so both answers are not pro-life)._ ![enter image description here](http://i.stack.imgur.com/3s6WK.png) X. are the three classes, the first horizontal line of numbers are the proportions of the sample in each class, and then the probabilities of those classes to say yes (=1) or no (=2) to the different item. I was wondering if the simple answer is: the items can not measure one construct because the model does not have an acceptable fit. Or is the structure of the latent classes more informative to answer this question? I know that putting constraints on the model could lead to an acceptable model fit, so in that way my criterium to base my statement on seems flawed. (something like: they can not measure one construct, but when we put constraint on the model, they can)


title: Multi-stage selection model with panel data in R

text: My data consists of individual level observations nested within countries over time. I would like to use multilevel models along with some sort of selection model. I have three related questions. 1) Are there any issues or concerns with using a Heckman selection model twice? I have a model with two selection stages. While I can't conceive of any reason why this shouldn't work, I also haven't seen any examples of a heckman selection model being used twice. Assume I also have the necessary exclusion restrictions. 2) My primary tool is R. There exists a package, sampleSelection, but my understanding is that algorithm it uses may not be appropriate for panel data / multilevel models. In stata, my understanding is the aforementioned is possible, but I'm not sure a package in R exists for doing heckman selection with multilevel models. 3) Are there better ways for accomplishing what I am trying to do? Perhaps using instruments alone would be better? If so, are there any R packages you could point me to? If there are no R packages, I am still interested in this answer since I would be happy to code something if necessary. Thank you.


title: Asymmetry between high order and low order interaction terms

text: Suppose we have a distribution for $\mathbf{x}\in \\{-1,1\\}^n$ and are interested in representing $p(\mathbf{x})$ as a linear exponential family with sufficient statistics of the form $1,x_1,x_2,\ldots,x_1 x_2,\ldots,x_1 x_2 \cdots x_n$ Suppose we know that ~~entropy~~ some measure of complexity of target distribution is bounded. We can view the parameter vector as a probability distribution $D$. Can we say anything about how the mass in $D$ is distributed among terms of various orders for some interesting complexity measure? Are there complexity measures that force mass onto low-order terms? Restricting attention to graphical models means discarding 3rd and higher order terms. I'm interested in some explanation why this works so well This was motivated by a similar question on CSTheory **Edit 11/29** Below is the probability simplex of of distributions over $\\{-1,1\\}$ valued $x_1,x_2$, in "probability" coordinates and in "log-linear" coordinates. In second picture, x,y axes correspond to coefficients in front of $x_1$ and $x_2$, z axis is the coefficient of $x_1 x_2$ term. You can see first and second order dimensions are symmetric, so bounding entropy (or any other symmetric function of distribution) is not enough to force mass onto low-order coefficients.... ![](http://yaroslavvb.com/upload/save/stats.SE.simplex.png) Mathematica source


title: What to do in case of nonsensical (negative survival time) predictions?

text: I have a rather small dataset of clinical variables and survival times. I've tried running linear regression on these in a cross-validated fashion. I note that in some cases, the regression will give a negative prediction. As these are survival times, this is nonsense. I am wondering if there is a generally accepted fix for this type of situation (set negative predictions to zero?) I had considered using a categorical predictor of some kind, but my dataset is too small to populate a meaningful number of categories. Thanks for any pointers or suggestions.


title: Intercept only logistic GEE model

text: I have repeated measures binomial data for several subjects across several experimental sessions. I am interested in testing whether or not across all experimental sessions the group proportion differs from 0.5. Does it make sense to fit a GEE glm with logistic link function and run an intercept only model and test if this intercept differs from 0? For example, the R code I ran is               gee.logit.0 <- geeglm(cs ~ 1, id = interaction(session,subject),std.err='fij',family = binomial(link = "logit"),data=data)      Which resulted in               Coefficients:             Estimate Std.err Wald Pr(>|W|)         (Intercept)   0.4212  0.0557 57.3  3.8e-14 ***      May I thus claim that the proportion I am interested is significantly different from 0.5 (since exp(0.4212)/(1+exp(0.4212))= 0.604) with piece of mind that I corrected for the within-subject correlation? Thank you for your help!


title: Why Phillips-Perron test returns "singularities in regression"?

text: I'm doing an unit root test using **Phillips-Perron** in the _tseries_ library (pp.test). I tried this code:               > pp.test(c(1:1000))      and the result is:               Error in pp.test(c(1:1000)) : Singularities in regression      doing a research I found the lines on the pp.test function with this error:               if (res$rank < 3)        stop ("singularities in regression")      _where **res** is the linear model (lm)._ Why i get this kind of error? If I do the same thing with another unit root test (in the same library) like ADF I don't get errors. Thank you


title: Correlations for nominal data

text: I have a gender variable in my questionnaire, I have done it where 1= Male and 2= Female. I would like to correlate this variable with their answers to compare the difference between male opinions against female opinions. could someone help me with which correlation to run please and how to understand it.


title: Averaging LASSO coefficients for repeated random partitioning of data

text: Is it reasonable to average LASSO coefficients from repeated reshuffling of training/test sets? Suppose I randomly divide my data into testing & training sets, then within the training set use 10-fold cross-validation to choose an optimal $\lambda$, then refit on full training data and record the model coefficients. Now, suppose that I repeat this process $k$ number of times. Each iteration will choose slightly different coefficients. One might think to average these together. However, each iteration may not choose the same set of non-zero coefficients, therefore the average of all coefficients may contain many more non-zero coefficients than any single solution. I found this similar brief discussion here, but I do not want to be confused by the additional discussion of multiple imputation: Combining LASSO coefficients across imputed datasets Please note that no answer was ever accepted for this question.


title: wilcox.test usage in R

text: I would like to know if there is a difference in the values taken from measurements of a given characteristic sampled at primary `forest=1` and secondary `forest=0` I'm using two vectors, one with the data itself (`DATA$Measurements`), and the other one with the category of sampling point (`DATA$Point`). How should I perform the `wilcox.test` fucntion? I'm doing this way:               wilcox.test(DATA$Measurements~DATA$Point)      But I notice that are other forms of performing this test:               wilcox.test(DATA$Measurements,DATA$Point)     wilcox.test(DATA$Measurements,DATA$Point,paired=TRUE)      Which one is more suited for testing my question of interest, and what is the difference between these options? (Is there a difference between points for this given characteristic)?


title: Random forest for binomial data

text: I have data where the response is of the form (k,n) - k successes, n trials for a given subject. The covariates are just some subject-specific predictors. Is it possible to use randomForest for this? I have thought about replicating the data so that for each subject we have n lines, k of which have the y-value '1' and n-k of which have the y-value '0' - is there an alternative? This makes cross-validation pretty messy, because I want to be able to predict p=k/n for new subjects, and I think the cross-validation would have to be 'locked' so that each training set contains all the rows for a given subject, or none of them.


title: Solving regularized least squares problems using Matlab optimization toolbox

text: I am trying to solve a least squares problem where the objective function has a least squares term along with L1 and L2 norm regularization. I am unable to find which matlab function provides the ability to perform such an optimization in addition to specifying constraints. I looked into the MATLAB optimization toolbox which also does not provide too much freedom to specify own objective functions (although I hope I am wrong in this case!) in the case of leastSquares optimization problems. If anyone knows how to model such an optimization problem in matlab please do help me out. The optimization problem is as follows: ![enter image description here](http://i.stack.imgur.com/W8YP6.png) Here **aj** is an m X 1 vector and **A** is an mXn matrix , **wj** is an n x 1 vector. I don't know how to incorporate the two additional regularization terms into the optimization problems as they only accept 2 matrices and perform the least squares operation on them. I would like to know how to implement this objective function and pass it into matlab optimization toolbox leastSquares function. Any help would be appreciated.


title: How to interpret probability plot?

text: I was reading about how to find the distribution of data from this question I generated the following probability plot from the scipy package.scipy prob plot Can someone please help me figure out how to go about interpreting such plots? The default distribution parameter is norm which is the normal distribution. The data is the number of packets sent from machine X to machine Y. I understand the following: the image does not follow normal distribution, but what else can I say about the data for the horizontal lines in terms of some well-known distribution if any exists. ![enter image description here](http://i.stack.imgur.com/1EdBr.png)


title: Comparing sample to a known population

text: Suppose you have recently run a survey on a known population (N=1000). You surveyed a subset of that population randomly (n=250). Now, suppose that you have a number of variables that are known. Like geography, expenditures, number of sales. Are there statstical tests that takes into account the known variables from the population to identify the accuracy, or precision of your sample. In other words, is there a way, beyond MoE, to denote the accuracy of your test based on these known parameters?


title: How to come up with research questions in statistical theory?

text: (I think this question is appropriate as a community wiki? If not, please close it, no hard feelings.) Does anybody have links to webpages/PDF files of advice for grad students in statistical theory (or even better if econometric theory) trying to come up with good research questions? I have been able to find such advice for medical researchers, computer scientists, applied economists, etc., but not theoretical statisticians (or econometricians). (Other than asking professors in my department in person.) There is a corresponding CW on MathOverflow.


title: To what extent is the distinction between correlation and causation relevant to Google?

text: ### Context A popular question on this site is " What are common statistical sins?". One of the sins mentioned is assuming that "correlation implies causation..." link Then, in the comments with 5 upvotes it is suggested that: "Google makes $65B a year not caring about the difference." At the risk of over-analysing a light quip, I thought this might be a useful discussion point for fleshing out the distinction between correlation and causation and the practical relevance of the distinction; and perhaps it could highlight something about the relationship between machine learning and the distinction between correlation and causation. I'm assuming the comment is addressing technologies that underlie the generation of search engine results and advertising display related technologies. ### Question   * To what extent is the distinction between correlation and causation relevant to Google's income generation, perhaps focusing particularly on the generation of income through advertising display related technologies and quality search results?


title: What is the mathematical difference between Cpk and Ppk?

text: I am trying to understand the mathematical difference between Cpk and Ppk used in Statistical Process Control (SPC). I have gone through the web but all I found a confusing and more confusing theoretical implication of the two. Before I can try to induce their application, I will need to find out What mathematical difference are? And even that I am not able to find a good source. I have googled and wiki-ed. My understanding is Ppk uses all data in calculating its standard deviation and Cpk uses all data but minus 1 in calculating its standard deviation. Is that it? For example, if I have a set of data "5" , "4" , "5.5" , "4.5" , "3.5" what is their Cpk and Ppk?


title: Weird formula for linear regression

text: I'll try to make the matter as clear as possible given the circumstances. My boss asked me to look at an old report a former employee wrote around a couple of months ago. Apparently the report contains some very useful information but, unfortunately, no one is able to understand it because it is written in "mathematics". The main node of the report is a linear regression formula for the predictions of future data for non-linear function. Consider to have the following: ![enter image description here](http://i.stack.imgur.com/hISMx.png) The Y-axis report a percentage (1 equals 100%). Now imagine to only have the part of the data before the red line. My former colleague figured out that the formula to estimate the second part of the data was: $$R(t) = 1 - d t^{1/\alpha}$$ While I did understand most of the formula I cannot understand what $d$ and $\alpha$ represent. By any chance, are they a common sign for something I am missing? Please, ask if something is not clear.


title: Area under the "pdf" in kernel density estimation in R

text: I am trying to use the 'density' function in R to do kernel density estimates. I am having some difficulty interpreting the results and comparing various datasets as it seems the area under the curve is not necessarily 1. For any probability density function (pdf) $\phi(x)$, we need to have the area $\int_{-\infty}^\infty \phi(x) dx = 1$. I am assuming that the kernel density estimate reports the pdf. I am using integrate.xy from sfsmisc to estimate the area under the curve.               > # generate some data     > xx<-rnorm(10000)     > # get density     > xy <- density(xx)     > # plot it     > plot(xy)      ![plot of the density](http://i.stack.imgur.com/9zT9s.png)               > # load the library     > library(sfsmisc)     > integrate.xy(xy$x,xy$y)     [1] 1.000978     > # fair enough, area close to 1     > # use another bw     > xy <- density(xx,bw=.001)     > plot(xy)      ![density with bw= .001](http://i.stack.imgur.com/agW3K.jpg)               > integrate.xy(xy$x,xy$y)     [1] 6.518703     > xy <- density(xx,bw=1)     > integrate.xy(xy$x,xy$y)     [1] 1.000977     > plot(xy)      ![density with bw = 1](http://i.stack.imgur.com/6QQ6s.jpg)               > xy <- density(xx,bw=1e-6)     > integrate.xy(xy$x,xy$y)     [1] 6507.451     > plot(xy)      ![density with bw=1e-6](http://i.stack.imgur.com/SZjuJ.jpg) Shouldn't the area under the curve always be 1? It seems small bandwidths are a problem, but sometimes you want to show the details etc. in the tails and small bandwidths are needed. **Update/Answer:** It seems that the answer below about the overestimation in convex regions is correct as increasing the number of integration points seems to lessen the problem (I didn't try to use more than $2^20$ points.)               > xy <- density(xx,n=2^15,bw=.001)     > plot(xy)      ![density with higher number of points to sample at](http://i.stack.imgur.com/f4rob.png)               > integrate.xy(xy$x,xy$y)     [1] 1.000015     > xy <- density(xx,n=2^20,bw=1e-6)     > integrate.xy(xy$x,xy$y)     [1] 2.812398      * * *


title: How to use a parametric statistical test with R?

text: I would like to test the null hypothesis that all coefficients of my categorical variable with 4 Levels are really zero using Wald test & R. I'm new to the R family. And just have some problems to implement that question. My model is `lm(Y ~ X + Z)` where `X` is continuous and `Z` qualitative/categorical with 4 levels. Example: So I run `model<-lm(prestige~height(cm)+weight(kg)+education)` where `education` is `yes` or `no`. As I saw, I would continue to with `coef(model)` and afterwarts test the $H_0$ as stated above. Then I would use `linearHypothesis(model,c("height(cm)=weight(kg),"education(yes)=education(no)+1")` And I would get the following output: `Linear hypothesis test` `Hypothesis:` `height-weight = 0` `education(yes)-education(no) = 1` So my question is how to compute Z or in der example education with more than 2 Levels...


title: Non linear regression in R Exponential curve fitting

text: I am trying to fit an exponential model to some data. The data is:               Wavelength  aCDOM     350.01  0.80605     350.22  0.78302     350.43  0.78302     350.64  0.78302     350.85  0.78302     351.06  0.78302     351.27  0.78302     351.48  0.78302     351.68  0.75999     351.89  0.75999     352.1   0.75999     352.31  0.75999     352.52  0.75999     352.73  0.73696     352.94  0.73696     353.15  0.73696     353.36  0.73696     353.57  0.73696     353.78  0.73696     353.99  0.73696     354.2   0.73696     354.41  0.73696     354.62  0.73696     354.83  0.73696     355.04  0.73696     355.25  0.71393     355.46  0.71393     355.67  0.71393      I know that the best model fit this type of data is an exponential function in the for `y ~ a * exp(-b * c)` where `a` is the absorbance estimate at the reference wavelength, `b` is the spectral slope which is the value I am looking for, and `c` is the wavelength minus the reference wavelength (440nm). I have tried using the `nls` function by estimating the starting parameters following the method outlined in Fox (2002) but I keep getting the following error in `nlsModel(formula, mf, start, wts)`:               singular gradient matrix at initial parameter estimates      The formula I am using is `model1 <- nls(St0104 ~ a * exp(-b * (Wavelength-440)), start = list(a=0.1, b=0.0012), trace=T)`, I'm not even if I have calculated the starting parameters correctly. I have tried instead to plot the log of `Y` and do linear as well as polynomial regression, but I know this is not correct for this data and the residuals are horrible. I'm pretty new to R but I have tried everything I can but this is driving me nuts can anybody help please?


title: Calculation of seasonality indexes for complex seasonality

text: I want to forecast retail items (by week) using exponential smoothing. I'm stuck right now in how to calculate, store, and apply the sesonality indexes. The problem is that all examples I've found deal with a sort of simple seasonality. In my case I have the following problems: 1\. Seasons don't occur on the same week every year: they are movable. Mardi-gras, lent, easter, and a few others. 2\. There are seasons that change depending on the year. For example, there is a national holiday season. Depending on whether the holiday is close to the weekend, customers will or will not leave town. So it's kind of like having two seasons: one where customers leave town, and one where they don't leave town. 3\. Sometimes two (or 3) seasons occur at the same time. For example, we had "Mardi-Gras" season occurring at the same time as Valentine's season.   4\. Sometimes seasons change in duration. For example, "Halloween season" started earlier this year. Christmas is also another example, where it seems like every year we start earlier to carry the products. It seems to me that I need to find a way to set some sort of "seasonal profiles" which then, depending on the particular scenario are somehow added to obtain the correct seasonal index. Does that make sense? Does anybody know where I can find practical information on how to do this? Thanks, Edgard


title: Different Scoring with libsvm

text: I have been using libsvm for binary classification, e.g. given data $(x_1,\dots,x_n)$ and labels $(y_1,\dots,y_n)$ where $y_i\in\\{0,1\\}$. Now I am wondering if the prediction of my model can be tuned in such a way that the output is in $[0,1]$ e.g. the classification is more specific. Does anyone know if libsvm supports something like this? Cheers


title: UMVUE for normal distribution $\sigma$

text: Let $X_1,X_2,...,X_n$ be a random sample from a normal distribution with mean $\mu$ and variance $\sigma^2$.   I showed that $(\bar X,S^2)$ is jointly sufficient for estimating ($\mu$,$\sigma^2$) where $\bar X$ is the sample mean and $S^2$ is the sample variance. Then assuming that$(\bar X,S^2)$ is also complete I have to show that $$\sqrt{ n-1\over 2}{\Gamma ({ n-1\over 2})\over\Gamma (\frac n2)} S$$ is a Uniformly Minimum Variance Unbiased Estimator for $\sigma$. I think I have to use Lehman Scheffe theorem as $(\bar X,S^2)$ is jointly sufficient and complete for $\sigma$. But how can I find a function which is unbiased for $\sigma$ that contains both $(\bar X,S^2)$. I don't understand how to work when there's a **joint** sufficiency and completeness.


title: How to get the overall effect for linear mixed model in lme4 in r?

text: I would like to get the significant value and effect size of the independent variable in overall, rather than the normal output from lme4 in R. It is just like the thing people report when running ANOVA. Do you have any idea how I can get this?


title: How to measure/argue the goodness of fit of a trendline to a power law?

text: I have some data to which I am trying to fit a trendline. I believe the data to follow a power law, and so have plotted the data on log-log axes looking for a straight line. This has resulted in an (almost) straight line and so in Excel I have added a trendline for a power law. Being a stats newb, my question is, what is now the best way for me to go from "well the line _looks_ like it fits pretty well" to "numeric property $x$ _proves_ that this graph is fitted appropriately by a power law"? In Excel I can get an r-squared value, though given my limited knowledge of statistics, I don't even know whether this is actually appropriate under my specific circumstances. I have included an image below showing the plot of the data I am working with in Excel. I have a little bit of experience with R, so if my analysis is being limited by my tools, I am open to suggestions on how to go about improving it using R. ![alt text](http://i.stack.imgur.com/NhPf2.jpg)


title: Predicting on data consisting of many independent short time series?

text: I have a dataset which contains many (hundreds) of short (3-30 obs) timeseries of different lengths. Each series is currently represented as a number indicating how long it took to the next event. To give you an idea:               0, 30, 34, 39, 45, 47, 51, 51     0, 5, 5, 7, 40      I'm sure that the values do not depend on time of day/week/month/... and that the timeseries do not influence each other. I do expect that the previous times between events have predictive value for how long it will take to see the next event. The source of my data are many (those hundreds) of devices being monitored and all these devices are identical. What i'm trying to figure out is how to build a single model to make predictions on how long it will take for the next event to happen in a device. My problem is that there are models enough when having one time series or even multiple time series if they are aligned in time (also not the case with my data). **So what model to use that allows leveraging as much of the data as possible?** If necessary, i could limit my analysis to a subset of the data having equal amounts of data points in their time series.


title: How to find the $(a_n,b_n)$ for extreme value theory

text: In the solution to this question (Extreme Value Theory - Show: Normal to Gumbel), the OP asked for the sequence $(a_n, b_n)$ such that $\Phi(a_nx+b_n)$ converges to the Gumbel CDF. Not only did I not able to understand the derivation in the accepted answer (I don't see the derivation for the last paragraph), I am also curious to know how one would go about deriving the sequence for distributions other than the standard normal. For example, the sequences in examples 7.5 and 7.6 in this document (page 6) seems to be pulled out of thin air.


title: Are these distributions independent?

text: I'm trying to compare the medians of two non-normal distributions: $A$ and $B$.   * $A$ is the distribution of completion durations for one type of task   * $B$ is the distribution of completion durations for another type of task However, since these distributions are measured from workers that can complete tasks from both distributions,   * is it incorrect to assume that these distributions are independent? Ideally, completing task $A$ should not have any influence on the time it takes to complete task $B$. However, I can imagine where a worker might have difficulty completing an instance of task $A$, become disgruntled, and therefore take a longer time to complete task $B$. Any help will be appreciated. Thanks. **edits** :   * I'm analysing data that has already been collected (rather than designing the study).   * There are actually more than two distributions in my data, I just wanted to focus on two first.


title: Do some of you use Google Docs spreadsheet to conduct and share your statistical work with others?

text: I know most of you probably feel that Google Docs is still a primitive tool. It is no Matlab or R and not even Excel. Yet, I am baffled at the power of this web based software that just uses the operating capability of a browser (and is compatible with many browsers that work very differently). Mike Lawrence, active in this forum, has shared a spreadsheet with us using Google Docs doing some pretty fancy stuff with it. I personally have replicated a pretty thorough hypothesis testing framework (including numerous parametric and nonparametric tests) originally done in Excel in Google Docs. I am interested if any of you have given Google Docs a try and have pushed it to its limits in interesting applications. I am also interested to hear about the bugs or flaws you have encountered with Google Docs I am designating this question "for community wiki" denoting that there are no best answers for this. It is more of a survey than anything.


title: Literature on measuring speed and memory usage of software

text: In one of my experiments for measuring the performance of a specific type of software I measure the execution time and memory usage. Now I know that measuring these things can be extremely problematic since the operating system interferes with how the program is run and other programs/processes can be run during the execution of the benchmark. I remember the following best practices:   * Run multiple times   * Take the smallest time it took to perform the benchmark, not the average.   * Take the memory usage with a HUGE grain of salt and only talk about it when the memory used is large (if its in kilobytes or a few megabytes you can't say much about it) and when the difference is significant (at least a couple of tens of megabytes). While its fun that I can remember these best practices and that I can argue for the correctness of each of them I'd rather tie these best practices to some literature. However I have not yet found any that speaks specifically about speed/memory benchmarks. Does anybody know of papers or other material I can cite to substantiates these best practices? _(I'm not sure if this question belongs on cross validated. But since I asked a related questionhere I thought it would be best to start here. If this question fits better on academia, cs theory or stack overflow, please help me move it to the correct sub-site)_


title: Correct estimation of arguments for glmmLasso function

text: I am using `glmmLasso` for variable selection. In my case, `n` is slightly less than `p` and `p` are bioclimatic variables for different time periods, so are highly correlated. How do I choose the right values for the arguments: `lambda` and `control`. I have tried with different values of `lambda`, `maxIter` (even to 10,000) and `control` (`start` and `steps`). But, the algorithm never converges. The p.values of all the selected variables are 0 and that makes me wonder if the non-convergence is the cause. What factor could help achieve convergence in a reasonable number of iterations? Among the variables selected, since p.values are all 0 (the rest of the variables are NA), can I estimate the relative variable importance based on the `StdErr`? Also, the `lambda` resulting in the lowest `BIC` values selects too many variables as being significant. Would it be OK in this case,(as I am only doing a rough variable selection for modelling) to not worry about the BIC, but choose the lambda that gives me a reasonable number of variables that also make sense given my data. What is the `start` arguments in `control`, anyway? All zeros are not acceptable, and so are vectors below a certain length. What do the values and vector length depend on? Apologies, if that shows ignorance of the mathematics behind mixed-models... For details on the data, please refer to Variable selection using mixed-models (lme4) I went through the `GMMBoost`, also by Groll, but did not find something to guide me in this case. Has anyone used glmmLasso for analysis and faced similar situations? Glad to hear any suggestions


title: learning hyper parameters: are we allowed to touch the prior parameters after observing the data?

text: There are many algorithms/applications that aim to learn the hyper parameters i.e. the parameters of a prior distribution from the observed data. A typical algorithm works in an iterative function where we optimise some cost function for a given set of hyper parameter values and then the next step is to update the hyper parameters using the current estimates of the other model parameters. Perhaps the EM algorithm works this way? My question is more philosophical. I thought that the priors should be fixed and should not depend on the data at all as they reflect our belief before we observe any data. So is tuning the hyper parameters not slightly unprincipled?


title: CART & stationarity assumptions

text: I have a question regarding if there are any stationarity requirements for a regression tree. If I were to regress 2 share prices on each other (normal linear regression), it would be classified as spurious regression due to the non-stationarity of the data ((stochastic) trend in the share prices). So in order to correct that we would regress the log returns (log returns hopefully yielding a stationary series) on each other - please correct if wrong. But with a regression tree there are no distributional assumptions - nor t-tests, etc. So should I be doing a regression tree based on the log returns of the share prices or on the share prices themselves? Which would be more appropriate I initial thoughts are that using stationary series in the data would be better but seek more of an understanding why Thank you for any help


title: Regression where predicted values always under-predicts the dependent variable

text: I am new to this field and was wondering if there is any regression technique available in literature that always under-predicts the dependent variable (y). I am actually planning to implement two step modeling where the first step will involve regression of y over x's and in the second level, I will use the residuals from first step to perform a different analysis. The second step requires the response variable (the residual of model in first step) to be always positive. Thanks in advance for your help! Edit: I am trying to implement Bass Diffusion Model (BDM) on Sales data, but first I want to take out the effect of some marketing variables using multiple regression. Implementing BDM in second step will require the residuals of the first step to be positive.


title: What if time variable is not significant in longitudinal analysis, can we remove it in the model?

text: In my longitudinal data, I firstly build a model with two fixed effects, session.week and sync. The former one is just the time variable. I actually have two subquestions : _**(1) What if time variable is not significant? can I remove it?_** When I run the following code with R                 quality.Model.1<- lmer(quality ~  session.week + sync  + (1|group.name/student.id), data = dfSubset)       pvals.fnc(quality.Model.1)      What I get is the following                 $fixed                     Estimate MCMCmean HPD95lower HPD95upper  pMCMC Pr(>|t|)       (Intercept)     3.3166   2.9966     0.9498     4.6601 0.0154   0.0000       session.week   -0.0458  -0.0265    -0.2369     0.1867 0.7934   0.6352       sync            1.3079   1.6866     0.1026     3.7183 0.0244   0.0206            $random                        Groups        Name Std.Dev. MCMCmedian MCMCmean HPD95lower HPD95upper       1 student.id:group.name (Intercept)   0.3303     0.1417   0.1535     0.0000     0.3912       2            group.name (Intercept)   0.0000     0.1975   0.3523     0.0000     1.2247       3              Residual               0.6365     0.6798   0.6859     0.5408     0.8552      As you see, the time variable is not significant, but sync is significant (p<.05) If I run a summary of the above model, we got:                 AIC   BIC logLik deviance REMLdev      113.3 124.4 -50.67    96.63   101.3      Actually although this is a longitudinal study, I don't actually expect the time would change my independent variable. The non-significance result actually confirms my expectation. Since the time variable (session.week) is not significant, can I remove it in my model, so it becomes:                quality.Model.2<- lmer(quality ~  sync  + (1|group.name/student.id), data = dfSubset)      pvals.fnc(quality.Model.2)      The result of this models is :                 $fixed                     Estimate MCMCmean HPD95lower HPD95upper  pMCMC Pr(>|t|)       (Intercept)      3.091    2.897     1.5119      3.931 0.0032   0.0000                sync    1.413    1.691     0.3116      3.320 0.0092   0.0074            $random                        Groups        Name Std.Dev. MCMCmedian MCMCmean HPD95lower HPD95upper       1 student.id:group.name (Intercept)   0.3371     0.1480   0.1583     0.0000     0.3903       2            group.name (Intercept)   0.0000     0.1787   0.2955     0.0000     1.0023       3              Residual               0.6284     0.6721   0.6784     0.5373     0.8353      And if we summary the result, we got:                 AIC BIC logLik deviance REMLdev       108.7 118 -49.35    96.91   98.71      We don't lose the significance of sync, but the AIC has reduced. Can I conclude that my second model (without the time variable) is better? And can I just remove the time variable safely? It just made me feel wired to remove the time variable in a longitudinal study. _**(2) What if the the significance of sync is canceled when testing a model with interaction terms?_** Suppose we don't give up time variable, and I built a mixed model with an interaction term as follows:                 quality.Model.3<- lmer(quality ~  session.week * sync  + (1|group.name/student.id), data = dfSubset)       pvals.fnc(quality.Model.3)      and we got the following :                 $fixed                                  Estimate MCMCmean HPD95lower HPD95upper  pMCMC Pr(>|t|)       (Intercept)                  4.6887   4.3486      1.058      7.477 0.0200   0.0004       session.week                -0.3859  -0.3402     -0.953      0.313 0.2842   0.1744       sync                        -0.6694  -0.2381     -4.267      3.874 0.8706   0.6835       session.week:sync            0.5110   0.4610     -0.395      1.368 0.2916   0.2018            $random                        Groups        Name Std.Dev. MCMCmedian MCMCmean HPD95lower HPD95upper       1 student.id:group.name (Intercept)   0.3488     0.1465   0.1576       0.00     0.3941       2            group.name (Intercept)   0.0000     0.1861   0.3387       0.00     1.2614       3              Residual               0.6260     0.6808   0.6869       0.53     0.8444      **Nothing is significant at all, including the interaction term!!!** , after introducing the interaction term. If we do summary on it:                    AIC   BIC logLik deviance REMLdev          113.7 126.7 -49.85    94.92   99.71      The AIC is larger than both quality.Model.1 and quality.Model.2. _**(3)What to choose_** So, which model should I choose if I want to see the effect of sync on quality? If I remove the time variable as fixed effects, then it's done perfectly. If I don't remove the time variable, then the main effects model suggests the significance of sync, but the interaction model cancels out everything.


title: Removing attributes with few observations in R

text: I have roughly 15 variables / attributes characterizing 6k customers in my data set. As they are categorical I have transformed them into 1 attribute for each possible value (1-out-of-K coding). An example could be Region with values "A", "B" and "C", which is transformed into 3 variables: `Region_A`, `Region_B` and `Region_C`. The same goes for other variables such as the `Sales Channel`. After this transformation I now have around 70 attributes. I would like to examine if there are any significant 2-way interactions between the different variables with regards to a response variable (concerning `customer quality`) using logistic regression. For instance, it is interesting to see if there is an interaction between `Region_A` and `Sales Channel 1`. However, there are very many possible interactions and therefore I would like to start by removing some variables, which have very few observations connected to them. An example could be that only 3 customers come from `Region_A`. More specifically, I would start by removing all attributes that have 5 observations or less connected to them (out of 6k observations). However, I cannot find out how to do that. Thus I have the following questions:   1. Does my thinking make sense? Or should I approach the issue in another way?   2. How do I remove all attributes in a dataset which has fewer than 5 observations connected to them? The values of the variables are always 0 or 1 as the customer is either from `Region A` (=1) or not from `Region A` (=0).   3. After removing these variables there should be fewer interactions. However, it would still be quite a large amount. I would therefore also like to only examine interactions with 5 observations or more. I am thinking this could be done using a formula in the logistic regression, but can you help me how I would find the right variables for the formula?


title: Ideas for ratio estimation (X and Y from different surveys)?

text: We are trying to estimate the rate of event X per hour exposed to certain activities (e.g. commuting) by combining statistics from two multi-year cross- sectional surveys, one for the denominator and the other for the numerator. However we don’t really know how to do this and are reaching out for any ideas. The first question is, is it possible to estimate MEAN(X/Y) and VAR(X/Y), where X indicates cases and Y indicates hours of exposure? The typical ways to estimate ratio (e.g. delta method) may not work here, because X and Y are from two different surveys. We noticed that an old thread in other site (http://www.stata.com/statalist/archive/2011-01/msg00723.html) suggests theoretically it is possible to estimate MEAN( X /Y) even though X and Y are from different surveys. Would doing so be exceedingly complicated (e.g. re- weighting the X, re-weighting Y…etc)? So we came up with plan B: what about just estimate MEAN(X)/MEAN(Y), and hence VAR(MEAN(X)/MEAN(Y))? Certainly, by doing so we are not trying to calculate the average of individual rates, but simply the ratio between the total events divided by total exposure hours, a commonly-used method in epidemiological studies in estimating risk at the population level. Then main question is, under plan B how do we calculate VAR(MEAN(X)/MEAN(Y))? Also, is estimating MEAN(X)/MEAN(Y) as straightforward as we thought, i.e. just substitute MEAN(X) with the sample mean, and do the same thing for MEAN(Y)? Thank you kindly


title: Main effects and interaction in multivariate meta-analysis (network meta-analysis) in R

text: I would like to conduct a multivariate meta-analysis (or network meta- analysis) to investigate the effect of two factors (A & B, 1 & 2). For every combination of the two factors I have a number of studies that report results that describe the difference in effect size as well as the variance thereof (A1 vs B1, A1 vs B2, A2 vs B1, A2 vs B2). Every effect size comes from an separate study. I would like to investigate the main-effects as well as the interaction of my two factors. Here I provide a dummy example of the data in R that I create initially for trying out the gemtc package:               set.seed(100)     treat1 <- c(rep("A1", 10),rep("A1", 8),rep("A2", 15),rep("A2", 5))     treat2 <- c(rep("B1", 10),rep("B2", 8),rep("B1", 15),rep("B2", 5))     TE <- c(rnorm(10,5,1),rnorm(8,2,1),rnorm(15,3,1),rnorm(5,1,1))     seTE <- c(rnorm(10,1,1),rnorm(8,2,1),rnorm(15,0.5,1),rnorm(5,0.2,1))     studlab <- paste(rep("study_", 38), as.character(c(1:38)), sep="")     my_data <- data.frame(treat1, treat2, TE, seTE, studlab)     my_data      I would be thankful for a pointer to an appropriate R package (gemtc, netmeta, mvmeta ?) that could handle such data and how to model not only main effects but also interactions of my factors. Some hints how to get started would be fantastic...


title: How do I find out specificity/sensitivity/predictive value with only prevalence and # of positives

text: I understand that               specificity = True Positive / (True Positive + False negative)     sensitivity = True negative / (False positive + True negative)          positive pv = True positive / (True positive + False positive)     negative pv = True negative / (true negative + False negative)      However, assume I know only which percent of population is resulted as positive in a certain test. also, I do know the prevalence. How can I calculate specificity/sensitivity/predictive values in this case? the question can be reduced into:               fill blanks given x, y, z, w               test    diagnosis    sum     positive    ?         ?         x     negative    ?         ?         y     sum         z         w


title: Clinical trials significance

text: I have been a studying high energy physics (HEP) for the last few years but I recently started working on a project in medical imaging. I have been a little surprised (not entirely I was aware that 95% was commonly used) to find major studies reporting using 95% cl. In HEP the convention is not worry too much until 3 sd if not 5. I appreciate in the 'real world' such idealistic situations cannot always be created but even so 95% doesn't seem such a high level of confidence. What is the rationale behind this? Is it simply a pragmatic one in the interests of completing in a timely manner? I have already come across trial or two which appear, to say the least, couterintuitive in their findings. Can anyone recommend a good book to jump into the mathematics of this? I have borrowed a fairly qualitative introduction but things like power calculations and kaplan meier curves are essentially totally new to me.


title: Dichotomizing when variable is skewed

text: I have a meeting frequency variable and I want to determine whether it has an effect on hypothesized relationship, by 1) controlling for it and 2)examining whether it is a moderator. Unfortunately, it wasn't a very sensitive measure and distribution is skewed. The breakdown is: 10% daily, 5% twice a week, 52% weekly, 27% fortnightly and 6% monthly. Someone suggested that because of the distribution (52% in a single category), I wouldn't get anything meaningful out of using it as a ordinal/likert variable. They said, instead, I should dichotomize it as weekly or more (67%) or less than weekly (33%). I'm unclear on whether or not I should do this. The articles I have read on dichotomizing variables (e.g. MacCallum et al., 2002) say that it is only okay to dichotomize when "the distribution of a count variable is extremely highly skewed". Meeting frequency is skewed (standardised skewness = ~-2), but it is not **extremely highly** skewed. Thoughts? I'm thinking at the moment, I'm better off treating it like a likert variable, and then in the limitations section of discussion saying that the variable had low variance and was not very sensitive and this limits the validity of the findings. Any assistance on this would be great. **The article I read about this is:** On the practice of dichotomization of quantitative variables. MacCallum, Robert C.; Zhang, Shaobo; Preacher, Kristopher J.; Rucker, Derek D. Psychological Methods, Vol 7(1), Mar 2002, 19-40


title: How to calculate power (or sample size) for a multiple comparison experiment?

text: I collected data on 20 groups (with 30 elements each). A multiple comparison procedure (pairwise t-test with Holm correction) shows that in general there are three sets of groups: the high with 4 groups, the low with 2 groups, and the middle with the remaining 14 groups. Each set is not significantly different for the groups within but it is significantly different from the groups in the other sets. (This is a simplification, because there are some other significant and non-significant results for the extremes of each set, but I am making a simplification of the results so I can write a concise summary of the experiment both to you and to the readers of the paper.) If this result is going to be used for decision making, for example treating members of groups the middle set as equivalent, one must be sure that the results are "real" and not just due to the small sample size. Thus I need to calculate some measure of power (power = 1- the probability of accepting H0 when it is false) or some measure of sample size to show that either a new experiment with a larger sample size is needed, or that indeed the differences are "probably true". But statistical power of WHAT?   1. It is not of the whole 20 groups ANOVA, since that analysis rejected the null.   2. Should I run the ANOVA of the 14 groups in the middle set and calculate the power of that? But that seems will overestimate the power (or underestimate the needed sample size) since the extreme groups in the middle set are "almost" different.   3. Should I calculate the power for the least significant pairwise t-test in the middle group (with a Bonferroni corrected alpha)? But that will terribly underestimate the power since the two most similar groups are very likely "really" not different. Any ideas? Any references I can follow? What I know so far:   1. The R package pwr calculates the power or sample size for t-test, one way ANOVA, and other tests.   2. On the relative sample size required for multiple comparisons, by Witte, Elston AND Cardon discusses the use of the Bonferroni corrected alpha values in the calculations of sample size for multiple comparisons.  EDIT - Aug 2013 There has been some upvote movement in this question, so I decided to add some more information, or better clarification regarding this topic. I did not quite agree with the two answers posted. I do not think it is a data-mining/clustering problem. But probably I did not phrase the question correctly. That paper is published so I can not only point to it, here, but also discuss what I needed. In the paper I (and colleagues) discuss the differences between productivity and citations among different computer science subareas, based on a random sample of 30 researchers in each area. The paper includes a compact letter display that shows the significant **differences** between any two of the 20 CS subareas. But I wanted to show significant **equivalences** between the areas. That is when it is very likely that two areas have the same productivity or the same citations per paper, given the 30 sample points for each area. I know of equivalence tests (or Two One Sided Tests - TOST) - there have been some discussions in CV on that, but nowhere did I see multiple equivalence tests! My idea was to use power was that the definition of power = 1- the probability of accepting H0 when it is false is exactly what I need to state that two areas have the same productivity - I make the statement that they have the same productivity (H0) and that statement is true with "power" confidence level! I still do not know how to do that, and the paper has no statement of **probable equivalence** between some CS areas, which is in fact the more interesting result! I would again appreciate any comments or help.


title: How to add outliers to an existing data?

text: I want to test few similarity measures for outlier detection. I've got some data from UCI repository, for example: Breast-Cancer. Is there a smart way to add artificial outliers to an existing data? Thank you.


title: Simple introduction to linear models in R

text: It was hard for me to understand linear models in R. There are a lot of documents for the case, but many of them are technical manuals rather than teaching the concept. I found this article really simple and instructive, I hope it would be useful for the other people who have the same problem. Do you have any better suggestion?


title: Sliding threshold for hypergeometric test

text: I have a list of 20,000 genes with an associated p-value (the total number of balls in the urn). I also have a list of specific genes that are "targets" of interest (these are the white balls in the urn). I want to do a hypergeometric test to see if there's an over-representation of targets among the genes I call significant (i.e., overabundance of white balls from the balls I draw from the urn). Easy so far with `phyper` in R. But I can change the number of significant genes (number of balls I draw) by varying the p-value threshold I'm using to call something significant (.001, .01, .05, etc). What I'm wondering - is there an alternative test that examines an overabundance of "targets" (white balls) at the top of the list of genes ranked by p-value? That is, can I do this test without specifying a p-value cutoff? That is, something like a more generalized hypergeometric test where I don't have to specify the number of balls I draw, that computes a result for all possible number of balls drawn (genes selected)? Thanks.


title: How to give entire word a quality score to account for the phred score of each letter?

text: I have a word representing a DNA sequence and each letter of that word is a nucleotide base (A,T,G,C). Each letter has a phred quality score. The formula used for the phred quality score is $Q=-10\log_{10} P$ where $Q$ is the phred score and $P$ the probability of that letter to be incorrect( for example, there was a problem in DNA sequencing). My question is: **what is the best way to give the entire word a quality score, so that it takes into account the phred score of each letter?** One thing I though about was to take the average and another was to add the phred scores.


title: How  to test  whether subgroup mean differs from overall group that includes the subgroup?

text: How can I test whether the mean (e.g., blood pressure) of a subgroup (e.g., those who died) differs from the whole group (e.g., everyone who had the disease including those that died)? Clearly, the first one is a subgroup of the second one. What hypothesis test should I use?


title: How to obtain covariance matrix for constrained regression fit?

text: **Is there an easy way to obtain the covariance of the parameters from a constrained regression fit?** I am using the PCLS function in the MGCV package in R to fit the constrained regression, however I am open to other approaches. The constraint I am imposing is that the coefficients must be positive.