title: Difference in values of tf-idf matrix using scikit-learn and hand calculation

text: I am playing with `scikit-learn` to find the `tf-idf` values. I have a set of `documents` like:               D1 = "The sky is blue."     D2 = "The sun is bright."     D3 = "The sun in the sky is bright."      I want to create a matrix like this:                  Docs      blue    bright       sky       sun        D1 tf-idf 0.0000000 tf-idf 0.0000000        D2 0.0000000 tf-idf 0.0000000 tf-idf        D3 0.0000000 tf-idf tf-idf tf-idf      So, my code in `Python` is:               import nltk     import string          from sklearn.feature_extraction.text import TfidfVectorizer     from nltk.corpus import stopwords          train_set = ["sky is blue", "sun is bright", "sun in the sky is bright"]     stop_words = stopwords.words('english')          transformer = TfidfVectorizer(stop_words=stop_words)          t1 = transformer.fit_transform(train_set).todense()     print t1      The result matrix I get is:               [[ 0.79596054  0.          0.60534851  0.        ]      [ 0.          0.4472136   0.          0.89442719]      [ 0.          0.57735027  0.57735027  0.57735027]]      If I do a hand calculation then the matrix should be:                           Docs  blue      bright       sky       sun                 D1    0.2385    0.0000000  0.0880    0.0000000                 D2    0.0000000 0.0880     0.0000000 0.0880                 D3    0.0000000 0.058      0.058     0.058       I am calculating like say `blue` as `tf` = `1/2 = 0.5` and `idf` as `log(3/1) = 0.477121255`. Therefore `tf-idf = tf*idf = 0.5*0.477 = 0.2385`. In this way, I am calculating the other `tf-idf` values. Now, I am wondering, why I am getting different results in the matrix of hand calculation and in the matrix of Python? Which gives the correct results? Am I doing something wrong in hand calculation or is there something wrong in my Python code?


title: Feature selection and cross validation

text: I'm working on a project and I would like to know if the following strategy is good/correct. Sorry if this is a basic/stupid idea (I'm new to this). The input is a dataset with 2.500 features and 1.000 instances. I have to apply a feature selection on this set. In the beginning I randomly make a learning set (which contains 70% of the original data) and a test set (with the remaining 30%). Now I split the learning set in 10 folds (to do cross-validation), so we've split the learning set into 10 (more or less equal) parts and now we take 9 parts out of 10 to do a feature ranking (using chi squared), then we take the best 10 features (highest ranking) to do a model evaluation (using SVM) on the remaining part of the 'splitted learning set' (and we do this 10 times, so each part will be used as a 'test set'). Doing this '10 folds' argument gives us 10 feature rankings (with 10 model evaluation/accuracies). Out of these 10 results we can take the best feature ranking (according to it's accuracy based on the model or to take the sum of the feature rankings) and then finally test it on the 30% test set from the beginning. Please give your opinion whether this is good or really stupid...


title: How to calculate precision and recall for multiclass classification using confusion matrix?

text: I wonder how to compute precision and recall using a confusion matrix for a multi-class classification problem. In specific, one observation can only be assigned with most probable class / label. I would like to compute   * Precision = TP / (TP+FP)    * Recall = TP / (TP+FN) for each class, and then compute the micro-averaged F-Measure.


title: Using panel OLS coefficients for forecasting - gretl

text: I conducted a panel OLS regression using gretl, on a balanced panel with around 25 zones and 10 time periods. I would like to use the model to forecast future values of y (dependant variable), given known future changes of x (explanatory variable). I would like to ask: how to deal with (a) the constant, (b) the unit dummies; (c) the time dummies when it comes to forecasting. Do i have to include any/all of the above in computing y-hat, or simply use the coefficient for the x term?


title: What is the expectation of one random variable divided over another (both independent)?

text: Suppose I have $X,Y$, which are independent random variables. Why is it that $E(\frac{X}{Y}) = E(X)E(\frac{1}{Y})$? Also, why is it that $E(X^2Y^2)=E(X^2)E(Y^2)$? How is it that the square of an independent random variable is also independent in relation to $Y$ or $Y^2$? Thanks!


title: Structural equation model with very small samples

text: I have two sets of a longitudinal data that I hypothesize to measure same latent construct.I am trying to test this hypothesis using Structural equation modelling technique. Basically, I am trying to use confirmatory factor analysis in longitudinal data using SEM. My models are relatively simple CFA models. The problem is that my sample size is small in term of SEM standards. We collected data only from 20 individuals for nearly two years (number of time points nearly equals 18) because the data collection is very expensive. I am new to SEM and I was not able to decide if SEM is right choice for me. My question is given such a small sample size, is it a good idea to use SEM based techniques to model the data and test my hypothesis? Based on my researches, I concluded that I should not fit SEM using least squares, weighted least squares or Maximum likelihood since they require large samples which I do not have. It seems that they recommend fitting SEM using Partial least squares or Baysian approaches for small sample size but the problem is my sample size is very small. Are there any articles that discuss fitting SEM using very small sample size? I am looking for suggestions, literatures or recomendations to fit such data with small sample size.


title: Kernel estimation of hazard or density function

text: I am working with R package "muhaz" to get the kernel estimate of the hazard function. If we have the following example which is the sample of uncensored data:                   x =c(0.001409591, 0.008451294, 0.019719245, 0.046626865 ,0.244813744,            0.410032119, 0.489286425, 0.697342699, 0.731541987, 0.954493031)         fit = muhaz(x, rep(1, length(x)), bw.method="local", b.cor="both", n.min.grid=51,         n.est.grid=101, kern="epanechnikov")      I have two questions: 1) How can I write the estimated hazard function as a formula in terms of the result of fit? How can we get the estimated hazard function for any desire value not just for grid points (fit$est.grid points)? 2) Is it possible to get the estimated density function? because as far as I know the result does not provide it. How can we get the estimated density function for any desire value? not just for grid points (fit$est.grid points).


title: On the hardness of data to learn

text: Almost in all texts which are discussing theorems of statistical learning, they assume analyzing arbitrary unknown distribution (the worst case). But in practice different problems (different data) have different levels of hardness, for example linear separable data is easier to learn than the data that are less (or not) separable by hyperplanes. Are there any works on formalizing the hardness of data (similar to the stuff done in complexity theory) Note: I am not sure if the same as analyzing VC-dimension or not, but I think it is not.


title: How to generate random points in the volume of a sphere with uniform nearest neighbour distances

text: With respect to post (1) and post (2), I generated a large number of uniformly distributed points inside the ball of radius $R$ using $\frac{R_s U^{1/3}}{\sqrt{X_1^2 + X_2^2 + X_3^2}} (X_1, X_2, X_3)$, where $U$ is uniformly distributed between 0 and 1, and $X_1, X_2, X_3$ are independent normal random variables with mean 0 and variance 1. The following figure shows a uniform spherical distribution obtained by this method using 10000 independent draws in a sphere of radius 10. ![enter image description here](http://i.stack.imgur.com/kzBDU.png) By computing the nearest neighbour distance $d_i$ of every point, I observed that the diagnostic plot of nearest neighbour distances does not follow a uniform distribution. Does this non-uniform distribution mean that one can cluster the points? Does it mean points dont have spatial randomness? If so, then how can I generate random points with uniform nearest neighbour distances. * * * Temporary images for @Anony-Mousse consideration: ![enter image description here](http://i.stack.imgur.com/txV1K.png) ![enter image description here](http://i.stack.imgur.com/1wLfB.png)


title: Does Bayesian data analysis take into account estimates' bias?

text: For example, the standard deviation is known to have bias that depends on the number of samples observed. If I wanted to do Bayesian inference on the SD of samples from two populations, and have unequal sample sizes from each, would the resulting inference take this bias into account? (Yes, I'm aware that the bias in the standard deviation doesn't really manifest discernibly unless you have very small N, but I'm more interested in the theoretical question. Plus, I do have hierarchical data with lots of units of observation, but few observations within each unit and a systematic bias for more observations in one particular experimental condition than another, causing me trepidation in my desired analysis of the SD within each level of unit-crossed-with-condition.)


title: Use daily production as a check on whether monthly production goal is likely to be met

text: This is a real small business problem, not a homework problem. A monthly production goal is 60000. We want to be sure that say 95% of months meet this goal. The goal is a minimum, so we want to be sure we aren't under this, but we don't care how much we are over. We want to use daily production as a check on whether we are likely to meet the monthly production goal. We don't really care how low one day's production is as long as other days make up for it, but if we are too low one day it should be a signal that there will be a problem in meeting the monthly production goal. How do we determine what is too low for daily production? The simplest approach is to just divide the monthly production by the number of workdays (say 16) to get daily production goals. This works, in the sense of if you meet each daily production goal you will meet the monthly production goal, but it ignores the natural variability in daily production. It is too severe. We should be able to have daily production that is some amount less than this without having to be worried yet. You can assume a normal distribution of daily production with errors that are independent from day to day, and that we have a history of daily production from which we can estimate the variance in daily production. Thanks for your help.


title: Finding the expected value of two normal random variables

text: Suppose $a_1 = b + c_1$ and $a_2 = 2b + c_2$ where $b, c_1, c_2$ are all $N(0,1)$ Find $E[b|a_1,a_2]$ My attempt: As $E[b] = 0$, I assume $E[b|a_1, a_2] = 0$. Is this a logical assumption?


title: Likelihood Ratio Test and Wald test provide different conclusion for glm in R

text: I'm reproducing an example from Generalized, Linear, and Mixed Models. My MWE is below:               Dilution <- c(1/128, 1/64, 1/32, 1/16, 1/8, 1/4, 1/2, 1, 2, 4)     NoofPlates <- rep(x=5, times=10)     NoPositive <- c(0, 0, 2, 2, 3, 4, 5, 5, 5, 5)     Data <- data.frame(Dilution,  NoofPlates, NoPositive)          fm1 <- glm(formula=NoPositive/NoofPlates~log(Dilution), family=binomial("logit"), data=Data)     summary(object=fm1)      * * * **Output** * * *               Call:     glm(formula = NoPositive/NoofPlates ~ log(Dilution), family = binomial("logit"),          data = Data)          Deviance Residuals:           Min        1Q    Median        3Q       Max       -0.38326  -0.20019   0.00871   0.15607   0.48505            Coefficients:                   Estimate Std. Error z value Pr(>|z|)     (Intercept)      4.174      2.800   1.491    0.136     log(Dilution)    1.623      1.022   1.587    0.112          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 8.24241  on 9  degrees of freedom     Residual deviance: 0.64658  on 8  degrees of freedom     AIC: 6.8563          Number of Fisher Scoring iterations: 6      * * * **Code** * * *               anova(object=fm1, test="Chisq")      * * * **Output** * * *               Analysis of Deviance Table          Model: binomial, link: logit          Response: NoPositive/NoofPlates          Terms added sequentially (first to last)                             Df Deviance Resid. Df Resid. Dev Pr(>Chi)        NULL                              9     8.2424                 log(Dilution)  1   7.5958         8     0.6466  0.00585 **     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      * * * **Code** * * *               library(aod)     wald.test(b=coef(object=fm1), Sigma=vcov(object=fm1), Terms=2)      * * * **Output** * * *               Wald test:     ----------          Chi-squared test:     X2 = 2.5, df = 1, P(> X2) = 0.11      Estimated coefficients are perfectly matching with the results given in the book but SE's are far apart. Based on LRT test the slope is significant but based on Wald and Z-test slope coefficient is insignificant. I wonder if I miss something basic. Thanks in advance for your help.


title: What are the algorithms to detect bots clicking on a shortened URL?

text: I own a url shortening service. I want to deliver only legitimate statistics to my clients. There are possible scenarios that a particular user writes a script to automatically open the shortened URL, thus making the statistics look bad. What are the approaches one can follow to detect if a click is legitimate or not? Please point me to appropriate resources like tools or research papers. The very basic approach that I can think of is to monitor the IP address of the user and block if the number of requests exceed a threshold.


title: What does "units" mean in epiR package?

text: From epiR manual for **_epi.2by2_** function: > units: multiplier for prevalence and incidence estimates. In examples it is set to 100, `units=100`. What does it mean? **Note:** I am working on ~5K case ~5K control data, no stratas, outcome vs exposure - 2by2 tables. **Update:** This PDF might be relevant.


title: Can I reconstruct a normal distribution from sample size, and min and max values? I can use mid-point to proxy the mean

text: I know this might be a little ropey, statistically, but this is my problem. I have a lot of range data, that is to say the minimum, maximum and sample size of a variable. For some of these data I also have a mean, but not many. I want to compare these ranges to each other to quantify the variability of each range, and also to compare the means. I have a good reason to assume that the distribution is symmetrical around the mean, and that the data will have a Gaussian distribution. For this reason I am thinking I can justify using the mid-point of the distribution as a proxy for the mean, when it is absent. What I want to do is reconstruct a distribution for each range, and then use that to provide a standard deviation or standard error for that distribution. The only information I have is the max and min observed from a sample, and the mid-point as a proxy for the mean. In this way I hope to be able to calculate weighted means for each group, and also to work out the coefficient of variation for each group as well, based on the range data I have and my assumptions (of a symmetrical and normal distribution). I plan to use R to do this, so any code help would be appreciated as well.


title: Which group of an indepent variable is more likely to influence the dependent variable and by how much

text: Give the following table                                 Male    |   Female   |   Total     Likes              10     |     3      |     13     Indifferent        20     |     20     |     40     Dislikes           10     |     9      |     19     Total              40           32           72      Here the independent variable is gender which is nominal. The dependent variable which is ordinal is either likers, indifferent, or dislikes. I can use the Chi Square test or Exact Fisher to determine if a H0 such as _men and women are equally likely to be likers_ should be rejected. However, the question I would like to answer after this is: **which gender is significantly more likely than the other to be a liker? And, by how much?** For example: Men are x% significantly more likely than Women to be likers. I thought that the next step would be to perform a test between liker and each group: i.e. males and likers. But I got stumped because I thought I should be comparing males to females with respect to likers. Is there any standard way of approaching this question?


title: Assumptions implied by "pairwise marginal" parameterization of MRF

text: I'm trying to understand the assumptions of different parameterizations in a Markov network. In this case, I'm trying to understand the assumptions (and effects) that result from parameterizing pairwise potentials with pairwise frequencies (i.e., $\psi_{(a, b)} = p(a, b)$), which the authors use in: Mitrofanova, A. et al. Prediction of Protein Functions with Gene Ontology and Interspecies Protein Homology Data. IEEE/ACM Trans. Comput. Biol. and Bioinf., 8, 775–784. Let's say we have a 3 node network (a)-(b)-(c), where the RV are Bernoulli. By chain rule of probability (and independence), in all cases, $p(a, b, c) = p(a) p(b|a) p(c|b)$. First, if we define $p(a, b, c) = \psi_{a, b}(a,b) \psi_{b, c}(b,c)$, where $\psi_{a, b} = p(a, b)$ and $\psi_{b, c} = p(b, c)$, then it appears that $p(b)$ must be equal to 1 (or 0). What would happen in larger, loopy networks? Second, if we use the same number of free parameters, constrained so each local potential sums to one, and learned the parameters (but the potentials need not be marginals), then I calculate that $p(b) = \psi_1(b) \psi_2(b)$, where $\psi_1(b) = \Sigma_a\psi_{b,c}$. As long as I don't assume $\psi_1(b) = \psi_2(b) = p(b)$, then is it a valid, useful parameterization? Are there any assumptions? Third, let's now remove the constraint on the local potentials, and add in a global normalization constant instead. If I compare $p(a) p(b|a) p(c|b)$ to $\frac{1}{Z} p(a,b) p(c,b)$, I can see that $Z = p(b)$. That's good. If, however, I try to calculate $Z$ by summing over all combinations of $a, b, c$, then I get that Z equals $\Sigma_b{p(b)p(b)}$. What happened?


title: weight parameter elimination in neural network

text: I am a little new to neural network. I have two questions: 1. Can we use neural network in the small $n$ large $p$ situation? 2. Any regularization methods designed for parameter elimination? (Shrink some $\omega$ to 0). For example to add $L_1$ regularizer as Lasso? I know some methods for avoiding over-fitting, (weight decay, early stopping, ect.) but none of them shrink parameters to 0.


title: Visualizing two variables which have a very similar values

text: I have two columns of data, and would like to show their differences through the visualization approach. The current issue is that these two columns are in-fact very close to each other. In other words, I would like to zoom in the minor differences that are existed using a graphical way. Are there any suggestions? Thanks. Here is part of the data set. The third column is the difference between the first column and second column.               0.999998295 0.999998331 -3.56E-08     0.99999381  0.999986678 7.13E-06     0.999997682 0.999998494 -8.12E-07     0.999989429 0.999989062 3.67E-07     0.999998877 0.999998505 3.72E-07     0.999995327 0.999996709 -1.38E-06     0.999995779 0.999995146 6.33E-07     0.999997484 0.999996701 7.83E-07     0.999998829 0.999998076 7.53E-07     0.999997523 0.999998836 -1.31E-06     0.999998287 0.999996964 1.32E-06     0.999983543 0.999989135 -5.59E-06     0.999992705 0.999994249 -1.54E-06     0.999995663 0.999998168 -2.50E-06     0.999997263 0.999999195 -1.93E-06     0.999994474 0.999996713 -2.24E-06     0.999993058 0.999994063 -1.00E-06     0.999998295 0.999996262 2.03E-06     0.99999736  0.999998185 -8.25E-07     0.999983507 0.996875105 0.003108402     0.999996339 0.999996597 -2.58E-07


title: Fitting copula for discrete margins

text: How can I fit a copula for a bivariate vector of negative binomial and Bernoulli margins? I would prefer a Frank or Clayton copula.


title: Can I probe cross-level interactions without random slope in a hierarchical linear model?

text: I have a HLM model with significant variance in the level 1 intercepts across groups but no significant variance in the level 1 slopes across groups and find significant cross-level moderation effects. Does it make sense to interpret these or are random slopes a necessary condition for probing cross- level interaction effects?


title: Is data driven identification in a simultaneous equation model possible?

text: Suppose I have to estimate 3 models: $y_1 = y_2\beta_1 + y_3\beta_2 + X\beta_3 + u_1$ $y_2 = y_1\alpha_1 + y_3\alpha_2 + X\alpha_3 + u_2$ $y_3 = y_1\gamma_1 + y_2\gamma_2 + X\gamma_3 + u_3$ I account for the endogeneity of $y_1, y_2$ and $y_3$ with a approach which creates a additional variable which allows me to not exlude any other regressors. Suppose these additional regressors are $c_1$ for $y_1$, $c_2$ for $y_2$ and $c_3$ for $y_3$. Now the models will be: $y_1 = y_2\beta_1 + y_3\beta_2 + X\beta_3 + c_2\beta_4 + c_3\beta_5 + u_1$ $y_2 = y_1\alpha_1 + y_3\alpha_2 + X\alpha_3 + c_1\alpha_4 + c_3\alpha_5 + u_2$ $y_3 = y_1\gamma_1 + y_2\gamma_2 + X\gamma_3 + c_1\gamma_4 + c_2\gamma_5 + u_2$ Given that the influence of the $c_i$ in each equation is significant, identification is achieved. The reason is that in each equation one $c_i$ is missing (for example $c_1$ in the $y_1$ equation). It turns out that I cannot leave $c_1$ out of the equation regarding $y_1$ (same for all the other equations). If I would, the model would be seriously biased. My questions are: When I estimate the full model, i.e., all $c_i$ are included in every equation and some of the regressors contained in $X$ are not signficant and these regressors are all different for each equation, is identification achieved? I would say yes. But this approach does only achieve identification through the insignificance of a regressor. I'm not quite sure if this is the right angle to it.


title: Dataset for article recommendation with demographic data

text: I am experimenting with a recommendation systems for articles (news, magazines etc). Is there any dataset available for the same which includes the demographic information of users as well? I am actually trying to incorporate demographic information as well in recommending, hence looking for such a dataset. I really don't want to simulate such data as I think real data would give much more interesting insights than simulated one. Anyone aware of such? I found one http://contest.plista.com/ but unfortunately its has german dataset and demographic information seems to be missing. Would appreciate any information here.


title: Deriving posterior distribution of multivariate normal

text: I have a problem that I can't seem to solve: The prior: $\beta_h \sim \mathcal{N}(\bar{\beta},D)$. Somehow the posterior is $\beta_h \sim \mathcal{N}(M, \Omega)$, where $M=(D^{-1}+X'X)^{-1}(D^{-1}\bar{\beta}+X'p)$ and $\Omega=(D^{-1}+X'X')^{-1}$, where $X$ is a data matrix and $p$ is a data vector. How can I prove $M$ and $\Omega$ are as such? What would be the likelihood function? I have 7 similar 'derive' questions that make no sense to me.


title: What is the difference between hot deck and single stochastic imputation?

text: From what I understand hot deck imputation involves 2 parts 1) Choosing the donor pool: This based on variables related to the missing variable. Maybe using a regression technique to test for association or just ask a subject matter expert. 2) Choosing the value to be imputed from the donor pool based on characteristics such as distance or through random sampling techniques. For single stochastic imputation the process is 1) Determine the distribution of the dataset after listwise deletion of item nonresponding units 2) Use stochastic techniques to choose a value from the distribution. Both these processes seem very similar to me. What exactly is the difference between the two ? Thanks


title: Eviews vs R (Statistical Software)

text: Are there many features in Eviews that R misses? I have heard that especially when dealing with time series R is less extensive than Eviews. Is this true? Which of the two packages contains most statistical tools (or are they comparable?)?


title: When to correct p-values in multiple comparisons?

text: I'm afraid that related questions didn't answer mine. We evaluate the performances of >2 classifiers (machine learning). Our Null hypothesis is that performances do not differ. We perform parametric (ANOVA) and non-parametric (Friedman) tests to evaluate this hypothesis. If they're significant, we want to find out which classifiers differ in a post-hoc quest. My question is twofold: 1) Is a correction of p-values after multiple comparisons testing necessary at all? The German Wikipedia site on "Alphafehler Kumulierung" says that the problem only occurs if multiple hypotheses are tested on the same data. When comparing classifiers (1,2),(1,3),(2,3), data only partially overlaps. Is it still required to correct the p-values? 2) P-value correction is often used after pairwise testing with a t-test. Is it also necessary when doing specialised post-hoc tests, such as Nemenyi's (non-parametric) or Tukey's HSD test? This answer says "no" for Tukey's HSD: Does the Tukey HSD test correct for multiple comparisons?. Is there a rule or do I have to look this up for every potential post-hoc test? Thanks!


title: Normalisation for regression

text: First off, I know little about statistics, so some of this question may seem naive. I'm trying to perform linear regression to model the relationship between x and y where: -x is a company's daily stock volume on a date -y is variable that is taken from the same date, however is something unrelated to stock volume. It is the volume of activity on that wikipedia for that company. I assume that the variables need normalising. Specifically, x needs to be normalised as overall index volume fluctuates. My first thoughts were to divide the daily volume by the total index volume. I'll do the same with the y variable. I just wondered if this seems sensible? Thanks EDIT: I've just noticed a typo in the question, the Y variable description has changed.


title: Comparing the means of two time series

text: I have data sets of the returns of two indexes in the same market (two different sets of stocks constituting each index), with 496 observations for each. I want to compare if the means are statistically different. I believe the variances are different, so I think I have to check if the variances are statistically different first. How would I do these things?


title: Zero regression coefficient when correlations are not zero

text: I don't really have a motivation for this - but I was thinking about this and couldn't work it out. Suppose I have a random variables $X$ and $Y$ which are correlated. Is it possible that the partial correlation between $X$ and $X\cdot Y$ is zero after taking into account Y? In other words, would a regression of $X$ on $Y$ and $X\cdot Y$ possibly result in a zero coefficient on $X\cdot Y$?


title: When using the beta distribution as a prior distribution for binomial, why won't the distribution results match with the calculated probability?

text: Let's say I have 1 success in 4 bernoulli trials, and I wish to plot the distribution of the parameter $p$ of the corresponding binomial distribution. I'm using R. The probability of seeing 1 sucess and 3 failures in 4 tests for $p=0.25$ is, for these parameters:               > n <- 4     > p <- 0.25     > dbinom(1, n, p)     [1] 0.421875      To get the distribution for the parameter, I use a beta distribution $Beta(k+1, n-k+1)$. But when I try to calculate the value for $p=0.25$, I get a different result:               > k <- 1     > dbeta(p, k+1, n-k+1)     > [1] 2.109375      I tried to divide that by the function beta $B(k+1, n-k+1)$, but that did not work either:               > dbeta(p, k+1, n-k+1)/beta(k+1, n-k+1)     > [1] 42.1875      That's exactly 100 times the result I expected. I am certainly missing something - maybe something basic, I'm no expert. Another thing that intrigues me: `dbeta` region sums to 1 without dividing by the beta function.


title: Discriminant analysis with random effects

text: Is it possible to do discriminant analysis with random effects? Is there an R package for this? Context: I have habitat use data for two species of frogs from radio telemetry, but nested within 'species' the data is for individuals with highly autocorrelated data. For example, I have ~280 relocations for species 1, with a mean of 20 relocations per individual, and ~ 210 relocations for species 2, with a mean of 13 relocations per animal. Relocations per animal are highly autocorrelated. At each relocation I collected structural habitat data, and I would like to identify which parameters are most strongly identified with each species. My preference is to have this autocorrelation accounted for in the analysis, but my alternatives with DA, I think, are to:   1. Use mean values for each individual, however then I would be breaking assumptions of sample size as I would have fewer samples than variables (therefore not a viable alternative); or   2. Cut relocations from individuals with lots of data points.  I've done some searching through CranR, docs, etc, but haven't had much success. The MASS package does discriminant analysis and random effects, but RE in DA? I'd be grateful to be pointed in the right direction (or told that I'm completely out to lunch!)


title: What is the definition of one-sample Kolmogorov–Smirnov test?

text: I have a time series data set of photon arrival times from a detector and need to know whether the arrival time is uniform.It is a continuous distribution? I have calculated the maximum $D$ between the normalized CDF of the photon arrival times. Then what should I do? There is a "$Pr(k\le x)$" in the Wiki link. What does it stand for in my question? Could somebody tell me the basic definition of ks-test? In fact, I can calculate ks-test probability via scipy.stats.kstest. Does anybody know the meaning of `args` for a uniform distribution? And the two output values?               >>> stats.kstest(sample, 'uniform',args=(1,2,3,4,5,6))      (1.0, 0.0)      >>> stats.kstest(sample, 'uniform',args=(0.1,0.2,0.3,0.4,0.5,0.6))      (1.0, 0.0)      >>> stats.kstest(sample, 'uniform',args=(0.1,0.2))      (1.0, 0.0)      >>> stats.kstest(sample, 'uniform',args=(1,2))      (0.98999999999999999, 0.0)      >>> stats.kstest(sample, 'uniform',args=(1.1,2.1))      (0.98499999999999999, 0.0)       **about uniform** I think I should compare two cdf:the normalised real data and the cdf of a uniform distribution that the photons arrive uniformly.   Please take a look at my plot.The green points are from real data.There are hundreds of green points.Every green point means an arrival of a new photon.x-axis is time and y-axis is normalised,in fact,it is the percent(cumulative arrival photons/total photon number).The red straight line stands for the cdf of uniform distribution. stats.uniform picks random values every time.Is stats.uniform appropriate?   ![enter image description here](http://i.stack.imgur.com/yvUBz.png) **about ks-test** I took an image of the 2rd version Numerical Recipes in C.   What is the relation between 14.3.9 and wikipedia's Pr(k<=x)?Just an approximation,right?   The significance should be decided only by two values max(D) and sample size,right? Pr(k<=x) is the cdf of D?How to define x?   I can not get a consistent result with stats.kstest().   You mean if I use the real time data which is not normalised,I should use two- sample ks-test,right?   ![NumerialRecipes2rdCh14](http://i.stack.imgur.com/hGMz5.jpg)


I have conducted a survey which asks respondents to indicate the informational sources regarding a product. There are 8 information sources (e.g., internet, newspaper, brochure and so on) and people can choose more than one source. I want to use information source a predictor variable in a multiple regression. However, because people can choose more than one information source, it does not seem possible to to treat it as a single categorical variable. **How can a categorical variable where respondents can choose more than one response be used as a predictor in multiple regression?**


title: Comparing model fit with heteroscedastic data

text: I am developing a physiological test using R that requires some parameters optimised. In comparing the new method against the existing method, the values of individual readings correlate in a linear way but are very heteroscedastic - the variance approximately proportional to the mean. I would like to compare the effects of various parameter changes by examining the goodness-of-fit of the linear models compared with each other. Confidence intervals around the linear model parameters are not needed. Comparing R squareds does not seem appropriate given the heteroscedasticity. Would distance correlation or R squared after HCCM correction be appropriate ways to test if one method is superior to another ?


title: Cross-validation vs empirical Bayes for estimating hyperparameters

text: Given a hierarchical model $p(x|\phi,\theta)$, I want a two stage process to fit the model. First, fix a handful of hyperparameters $\theta$, and then do Bayesian inference on the rest of the parameters $\phi$. For fixing the hyperparameters I am considering two options.   1. Use **Empirical Bayes (EB)** and maximize the marginal likelihood $p(\mbox{all data}|\theta)$ (integrating out the rest of the model which contains high dimensional parameters).   2. Use **Cross Validation (CV)** techniques such as $k$-fold cross validation to choose $\theta$ that maximizes the likelihood $p(\mbox{test data}|\mbox{training data}, \theta)$. The advantage of EB is that I can use all data at once, while for CV I need to (potentially) compute the model likelihood multiple times and search for $\theta$. The performance of EB and CV are comparable in many cases, and often EB is faster to estimate. Question: Is there a theoretical foundation that links the two (say, EB and CV are the same in the limit of large data)? Or links EB to some generalizability criterion such as empirical risk? Can someone point to a good reference material?


title: Confidence interval(s) for randomization design

text: I am developing a method to do something with permutation testing. I have two permutation / randomization schemes and do not know which one is better than the other. Can I find the confidence interval for each permutation scheme? If so, is the one with the narrower confidence region better?


title: Sample size calculation for study aiming to demonstrate non-inferiority of a drug

text: I need to calculate the sample size required for an observational study in which incidence of disease is 19-29%. Population affected is 600,000 people. The study has two samples similar on baseline characteristics treated with two different drugs. Statistical analysis will be chi-square and Fisher's test. I have to demonstrate non-inferiority of one drugs with respect of the other one. I don't have any other info. Can anyone help me to calculate the sample sizes required in order to obtain power of 0.80 with an alpha 0.05? Thanks a lot in advance


title: What else can be deduced from the following class grade summary information?

text: Lets say for example a class of students and their grades are the data set. Lets say there are around 35 students. Here is what you know:               Your mark     The class average     The median mark     The standard deviation      Are there any other conclusions about this data that can be made given this information? Thanks


title: How do I test whether I can properly apply LDA?

text: I have some data which works nicely with JMP's canned linear discriminant analysis (LDA), but after reading about LDA I'm not sure if the analysis is valid. The Wiki article notes a fundamental assumption of LDA is that independent variables are normally distributed. Does this also mean that I cannot use LDA if my variables have any correlation whatsoever? Also, supposing I have or transform to independent variables-- then I must also test if they are normally distributed (say, Shapiro-Wilk test) and only use those that pass in the LDA? Thanks


title: Power calculation for disease modeling in iPS cells

text: I apologize in advance if the question is rather simple, but I have been having a difficult time figuring it out. I am an experimental scientist interested in modeling human disease with induced pluripotent stem cells. The disease is fairly common, but doesn't have a common genetic basis. It is mainly driven by rare genetic variants of high penetrance but in a small number of individuals. For (let's say) 5 patients out of 1000 screened affected cases, we are able to strongly associate mutations in one gene with the condition. These mutations are completely absent from a control population of 5000 individuals. We have made stem cells from 2 of those individuals to determine whether these mutations result in cellular dysfunction (let's call this phenotype A) when compared to controls. We have several stem cell lines per individual. In this situation, phenotype A is an 'intermediate phenotype'. The big conceptual question for me is this: Am I modeling the disease or am I only modeling the effect of the disease-associated mutations on cellular function? How can I determine the sample size to test the former claim with, say, 80% power?


title: "non-monotonic profile": what should I do when getting this when running confint(lmer)?

text: I am new to R and I have run the `lmer` in lme4. In the model summary, there is no any warning message, but when I ask for `confint(lmer)`, R gives me this message at the bottom of the output:               Warning messages:     1: In profile.merMod(object, signames = oldNames, ...) :       non-monotonic profile     2: In profile.merMod(object, signames = oldNames, ...) :       non-monotonic profile     3: In profile.merMod(object, signames = oldNames, ...) :       non-monotonic profile      What does it mean? and what should I do?


title: Dirichlet process/Chinese restaurant process for clustering in R

text: I recently read a fascinating article describing methods for clustering data without assuming a fixed number of clusters. The article even includes some sample code, in a mix of Ruby, Python, and R. However, the meat of the analysis is performed using scikit-learn's Dirichlet Process Gaussian Mixture Model to actually find clusters in some sample data taken from McDonald's menu. Obviously, this a a great excuse to learn some more python, but I'm lazy and would like to find a ready-made R package that can take a dataframe and return clusters, in a manner similar to the kmeans function. A quick search on CRAN reveals the packages dpmixsim and profdpm. Any suggestions for the best place to start?


title: Likelihood score function 101

text: I have some trouble with score functions in likelihood calculation. I'm not good at statistics or probability, so I'm still confused on formalism and mathematical-probabilistic language. **Some background:** I'm on particle filter, so inferring a pdf that is analytically intractable by a Monte Carlo sampling: put some random particles, evaluate them with a likelihood function, and the normalised-weighted set of particles can be an approximation of the pdf. I have to give a weight to each particle and compute the likelihood with a likelihood function or score function. If my hidden state described by this difficult pdf that I have to infer is denoted as $X$, and I have some observation $z$, the likelihood can be formalised by: $$ p( X | z ) = \text{likelihood} $$ **What I have now:** a nice distance function that give me the distance between the particle and the target. It works. It is the Bhattacharyya distance (I'm working with images, and colour histogram is a relevant feature). So that: $$ d = \text{Bhattacharyya(particle, target)} $$ This $d$ is a distance, and for a likelihood function I need a probability, so I'm putting it into a Gaussian with zero means and variance = 0.4: $$ \text{likelihood} = \frac{e^{-\frac{(d-\mu)^2}{2 \sigma^2}}}{\sqrt[2]{ \pi \sigma}} $$ **What I want:** Now I'd like to improve my experiments, adding similarly another score function so that the total likelihood is calculated with two score functions (let's say a shape-distance or something else) $s_1$ and $s_2$. I remember from school (but don't remember exactly why) that if you want to combine 2 probabilities, you have only to multiply, not sum. So I know that: $$ \text{likelihood} \propto s_1 s_2 $$ where $s_1$ is the Gaussian with the Bhattacharyya with the colour histogram, and $s_2$ is the shape distance. The two score functions are Gaussian, and they are not independent (they depend on the same hidden state). **Question 1** : Is it true I have to multiply them? Why? So that my final formula will be: $$ \text{likelihood} = \frac{e^{-\frac{(d_1-\mu)^2}{2 \sigma^2}}}{\sqrt[2]{ \pi \sigma}} \frac{e^{-\frac{(d_2-\mu)^2}{2 \sigma^2}}}{\sqrt[2]{ \pi \sigma}} $$ **Question 2:** What if I want to empirically add some more importance to one of the two scores because I think it is more descriptive of my data? Like: $$ \text{likelihood} = \alpha s_1 (1-\alpha)s_2 $$


title: How can I plot the components of a mixture model obtained from pymix?

text: I'm trying to wrap my head around mixture modelling, and I've come across a small matlab script that seems relevant. In order to familiarize myself with pymix, I've decided to try rewriting the matlab script in python. I appear to have successfully fitted/decomposed the data with the expectation maximization algorithm, but I'm not sure I understand the nature of the method's output. My question is as follows:   1. How can I replicate the second figure from the matlab code (below) with the output of `pymix.MixtureModel.EM`?   2. Can someone explain the output of the above method? I should mention that I'm horribly under-educated in statistics and mathematics in general (I'm working on it!), so please assume I know nothing =) In any case, here's the matlab script:               % Generate some data drawn from two Gaussians     data = [0.4+randn(100,1).*0.15; 1+ randn(200,1).*0.25]';     data(data < 0.05) = 0.05;     [n,x] = hist(data);     bar(x,n);          % Make the mixture model pdf     mixtureGauss = ...         @(x,m1,s1,m2,s2,theta) (theta*normpdf(x,m1,s1) + (1-theta)*normpdf(x,m2,s2));          % Set up parameters for the MLE function     options = statset('mlecustom');     options.MaxIter     = 20000;     options.MaxFunEvals = 20000;          % Get max likilihood parameters for our mixture model (start with some     %   reasonable guesses about the parameters)     p = mle(data, 'pdf', mixtureGauss, 'start', [0.5 0.1 0.5 0.1 0.5], ...         'lowerbound', [-Inf 0 -Inf 0 0], 'upperbound', [Inf Inf Inf Inf 1], ...         'options', options);          % Plot and print information     hold on;     x = linspace(min(data),max(data),100);     plot(x, mixtureGauss(x,p(1),p(2),p(3),p(4),p(5))*max(n), 'r', 'LineWidth', 2);     fprintf('Gauss 1: %0.2f (+/- %0.2f)\n', p(1), p(2));     fprintf('Gauss 2: %0.2f (+/- %0.2f)\n', p(3), p(4));     fprintf('Mix: %0.2f proportion first gaussian\n', p(5));      And here's what I've done in python so far. Note that I'm running this code in iPython with the `\--pylab=inline` option, thereby importing pyplot into my main workspace:               import numpy as np     import mixture          # Generate some data drawn from two Gaussians     data = np.concatenate((0.4 + np.random.randn(100) * 0.15, 1 + np.random.randn(200) * 0.25))          data[np.nonzero(data < .05)] = .05     print type(data), len(data)          plt = Figure()     hist(data, bins=50)     show()          #  Create DataSet object     mixdat = mixture.DataSet()     mixdat.fromArray(data)          # reasonable-guess mixture, akin to random starting point in K-Means     n1 = mixture.NormalDistribution(-2, 0.4)     n2 = mixture.NormalDistribution(2, 0.6)          # Mixture model and EM clustering     mix = mixture.MixtureModel(2, [.5, .5], [n1, n2])     postmat, _ = mix.EM(mixdat, 40, 0.1)          fig = Figure()     hist(data, bins=50)     x = np.linspace(np.min(data), np.max(data), 100)     # Now what?     show()      Any other comments, criticisms, or explanations are more than welcome. Thanks for putting up with my ignorance ;-)


title: How to compare species communities between two years?

text: I would like to compare the data of a community of species between two years. I used to compare communities based on different levels of a factor (for example, a treatment). The case of a treatment is easy because each site has a corresponding level of the treatment. However, if the factor is a year, we face a paired design with each site sampled twice. In this case I cannot just assign a year to a site and, therefore, I do not know which method I should apply to analyse community differences between two years.


title: References suggested for multivariate analysis of several similar time series

text: I have a time series dataset that reports the hourly page views and social media shares of online news stories. What I hope to obtain is the relationship between the two variables. I would imagine that the more shares a story gets, the more page views it will attract, and vice versa. That means the two variables would demonstrate a similar fast/slow growth curve. In terms of a single story, I can use VAR (vector autoregressive) model to estimate the coefficients. However, how could I make use of hundreds of stories to build a more generalized model? What I've thought of so far is to fit each story with the same model and then observe the distribution of the coefficients. Would it be an appropriate method? That would be great if you could kindly cite some papers and/or R packages that would inspire potential solutions. Thank you!


title: sparse coefficient time series models

text: I am trying to build sparse coefficient time series models but couldn't find a good resource/book to learn from. Can anyone share with me a good resource/book? It will be great if it comes with R codes.Can anyone use this kind model? Thanks.


title: About SVM cost and gamma parameters tuning

text: I am using `R` and `e1071` package to tune a C-classification SVM. My question is: regardless of the kernel type (linear, polynomial, radial basis or sigmoidal), is there any good criterion to choose the range in which cost and $\gamma$ parameters should range over and/or to choose what the granularity should be (that is, as an example, `gamma = 10 ^ (1:2)` or `gamma = 1:2` or `gamma = 100 ^ (1:2)`)? I add a second question: can `tune.svm()` return the best kernel type, too? Thanks,


title: How to incorporate constraints in random forest output

text: Suppose I am doing random forest classification of labels $A$,$B$,$C$,$D$. There is some theoretical ordering to this output such that when $A$ is more likely than $B$, $B$ is also more likely than $C$, etc. Also, if $P(D) > P(C)$, we also have that $P(C) > P(B) > P(A)$. There are other such conditions that need to be met. The issue is that a real random forest may give something silly that completely violates the above constraints, _even if_ it is able to predict the most likely outcome successfully. For my use case the ordering is important since decisions are made not only on the most likely outcome. It also seems intuitive that I should be able to improve generalization if I can somehow enforce this prior knowledge into the model. How do I account for this in a decision forest? Despite this structure to the output I do not think it is possible to construct a real-valued response variable since they are _still class labels_ with no natural real value, even if there is some type of ordering to them.


title: Jensen-Shannon divergence for finite samples

text: I have two _finite_ samples $s_1$ and $s_2$ and two distributions $p_1(s_1)$ and $p_2(s_2)$ that are associated to these samples. I'm essentially interested to measure the distance or similarity between these two distributions. I'm currently using the Jensen-Shannon (JS) distance, which involves in calculating the entropies of the two distributions separately and jointly. I need to bin the data for calculating the basic entropy functions and it looks like the resulted JS distance is also a function of the data binning strategy that I use. It sounds quite arbitrary and I don't know how much I can trust the results. Is there any way around this? How can I measure the distance between $p_1(s_1)$ and $p_2(s_2)$ in a more rigorous way that needs no data binning or any other sort of arbitrariness?


title: Reporting chi-square tests for weighted data

text: I was wondering if you could share your experience in reporting chi-square tests for complex survey data in journal publications. Normally Chi-square tests are reported as $\chi^2_1(2,\text{ N} = 90)= 0.89, \text{ p} = .35$. (for example, although I guess there can be some variations). However the concept of N becomes rather tricky with complex survey design, and I'm not sure whether I should report the probability weighted sample size or just the sample size.I'm also interested in whether people would report different Ns if the sample was just subject to probability weighting (e.g. for example if strata or cluster information was unfortunately missing) or probability weighting with clustering and stratification. And if the sample size was reported in the format above should this be reported as n rather than N? Just to give an example: Lumley's R package returns the following:               data(api)     dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)     svychisq(~sch.wide+stype,dclus1, statistic="Chisq")                    Pearson's X^2: Rao & Scott adjustment          data:  svychisq(~sch.wide + stype, dclus1, statistic = "Chisq")      X-squared = 11.94, df = 2, p-value = 0.005553      Would you report as N               sum(xtabs(~sch.wide+stype, data=apiclus1)) #or     sum(svytable(~sch.wide+stype, dclus1))      See this page for details on computation and Rao-Scott corrections to the Pearson chisquared test. Very interested in your opinion. Many thanks


title: Multivariate regression estimation when the variables' variances are known a priori / sourced seperately

text: I'm looking to use a multivariate regression for prediction, but making use of (possibly) superior estimates of variance for both the independent and extraneous variables. My approach is to standardize the dependent and extraneous variables (by dividing their respective standard deviations derived from the full data history). Once I have the standardized regression coefficients, I use these, together with my separately sourced variance estimates, to get back to an equation that will be used for generating the predictions. Do I simply multiply each coefficient by my separately sourced variance estimate? But what about the dependent variable variance estimate? As background info, my separate sourcing of variances estimates is because I believe I have better (more up-to-date) estimates than available from the full data history (let's say, because the variances aren't stable through time, I can estimate a more timely measure of variance using higher frequency data over a shorter, recent period).   * Is this approach a sensible standard practice?   * How to go from the standardized betas, to un-standardized ones using the high frequency variance estimates? I read here (http://stats.stackexchange.com/questions/29781/when-should-you- center-your-data-when-should-you-standardize) about WHEN to standardize, but not clear to me if this covers the case of standardizing both dependent and extraneous variables. Thank you!


title: Using central limit theorem for approximation

text: > Let $X$ be a random varaible from a distribution with pdf $$ f(x) = \theta > x^{\theta-1}, \quad 0< x < 1. $$ > > a) Name the distribution of $U=-\ln(X)$ by first finding its density > > b) Let $X_1, X_2, \ldots,X_n$ be independent and identically distributed > random variables with pdf given by earlier with $\theta$= 3. Using the > result from a) and by the central limit theorem (CLT) > > i) find an approximation to $P(X_1 \cdot X_2 \cdot \ldots\cdot X_{30} \leq > 1.85 \cdot 10^{-5})$ i.e. for the probability of the product of the r.v's. **My attempt** I found part a) by doing the transformation, and got an exponential with parameter $\theta$ where my pdf is: $$ f(x)=\theta e^{-u\theta} . $$ Now where I am struggling is with part i) where I am supposed to find an approximation using the CLT and together with the mean equal to $1/3$ and variance equal to $1/9$. I am able to show that with 30 observation we just take the product of every mean and variance from each individual observation giving us: mean = $(\frac{1}{3})^{30}$   variance = $(\frac{1}{9})^{30}$ and then we can just substitute these into the z score by using central limit theorem $$ P(X_1 \cdot X_2 \cdot \ldots \cdot X_{30} \leq 1.85 \cdot 10^{-5}) = P\left(\frac{X -\mu }{\sigma}\leq \frac {1.85 \cdot 10^{-5} - (\frac{1}{3})^{30}}{(\frac{1}{9})^{15}} \right) $$ Hence giving me an answer of $ P(Z \leq 3,808,985,943)$ which definitely cannot be correct. Would appreciate it if somebody could point out my mistake. **Reattempt:** Using the hint i got i managed to deduce the following $P(X_1 \cdot X_2 \cdot \ldots \cdot X_{30} \leq 1.85 \cdot 10^{-5})$ $P(log X_1 \cdot log X_2 \cdot \ldots \cdot log X_{30} \leq log 1.85 \cdot 10^{-5})$ P($\sum_{k=1}^{30} log X_i \leq log 1.85 \cdot 10^{-5})$ And since X random variable can be normally distributed X~ N( $\frac{1}{3}$ , $\frac{1}{9}$) Substitute into the Z score $ P(\frac{ log X -\mu }{\frac{\sigma}{\sqrt(n)}})\leq \frac { \frac{log1.85 \cdot 10^{-5}}{30} - (\frac{1}{3})}{(\frac{\frac{1}{3}}{\sqrt(n)})} $ $P(Z \leq 11.45)$ Am I on the right track? Is it correct to use $\mu$ = $\frac{1}{3}$ and $\sigma$ = $\frac{1}{3}$ in the z score?


title: Normalization of likelihood

text: If I'm not wrong, likelihood functions are sensitive to the size of the sample, i.e. the larger the sample, the lower the likelihood value. Given a sample $x$ of a random variable $X \sim f(\theta)$, and a parameter estimate $\hat\theta$, suppose I want to test the hypothesis that the likelihoods of different subsamples of $x$, let's call them $x_a$ and $x_b$ are equal. The problem is $x_a$ and $x_b$ have a different number of elements, say $n_a$ and $n_b$ respectively, and $n_a \neq n_b$, so I assume the likelihoods must be normalized in some way. Is it enough to divide by the sample size? For instance, can I use a test statistic such as $T = \frac{\ell(\hat\theta|x_a)}{n_a} - \frac{\ell(\hat\theta|x_b)}{n_b}$ and then test the hypothesis that $T = 0$. A related question would be how do I find the probability distribution of the above statistic, but that is another story.


title: Legitimate IV for panel?

text: In a panel regression using fixed effects, is it correct to use an instrument that does not vary by year? To use an old but simple example, Angrist and Kruger 1991 use quarter of birth as an instrument for education. If one had access to similar data that was in panel form, is it still correct (and I don't mean specifically quarter of birth, I know there are problems with the validity of this IV) to use quarter of birth (same for each respondent regardless of panel wave) as an instrument in the panel IV regression? Why or why not?


title: Can I use boostrapping, why or why not?

text: I am currently working on biomass estimates using satellite imagery. I'll quickly define the background of my question, and then explain the statistical question I am working on. **Background** _Problem_ I am trying to estimate biomass over an area in France. My response is the steamwood volume density (in $m^3/ha$), which is more or less proportional to biomass (depending on the wood densities...). The independent variables that I have are vegetation indices derived from measured reflectances over this area (the satellite used in the study is MODIS for those who know it). These indices are for example NDVI, EVI, etc. I have maps of the indices, and the resolution of the maps are 250m. There are strong correlations between these indices and the volume in a same forest type (biome and climate). So I am trying to regress the volume density against these indicators (actually their time series) on inventory plots where I know the volume. _Forest inventories_ The volume on these plots is estimated with the following sampling method:   1. Inventory nodes are placed on a regular grid covering the area.   2. A plot is attached to each node, and the inventory process (tree types, volumes, canopy height, etc.) occurs on this plot. Of course I am interested only in the inventory plot and the values of my vegetation indices is the value of the pixel containing the plot.   3. The inventory process on a plot is the following:  ![http://i.stack.imgur.com/DeHdC.png](http://i.stack.imgur.com/qcHtu.png)     * Measure of the trees that have a diameter > 37.5cm in the 15m radius circle     * Measure of the trees that have a diameter > 22.5cm in the 9m radius circle     * Measure of the trees that have a diameter > 7.5cm in the 6m radius circle The volume density is then calculated using expansion factors. For each plot I have access to the data for all the measured trees. Moreover, for each single tree, I have an uncertainty on the volume due to the use of allometric equations (let us say 10%). **Where statistics are important...** For my regressions to be more accurate, I need for each estimate of volume the variance/CI of this measure. This depends, IMO, on the number of trees sampled and the volume density found. So I have two problems:   1. How to account for the fact that my vegetation indices are measured over a pixel of 250m ? > I can assume that the volume density is constant over one pixel, and that I > sample this pixel with one inventory plot.   2. How to estimate the variability of my volume density?  > I think I could use bootstrapping on the tree population. But my total > number of trees measured can be pretty small (from 7 to 20...). Moreover, > how can I take into account the fact that I am measuring the trees on > different circles depending on their sizes ? And how should the variability > change if I am looking on an entire pixel ? > I was also thinking that I could use a Monte Carlo Simulation to simulate a > forest, and then randomly sample this forest with plots to see what is going > on... I do not have a strong statistical background, so I am a little bit lost!


title: How is a prediction interval generated for many input observations?

text: I am aware that the "prediction interval", as defined in most textbooks on linear models, is focused on the uncertainty in the model being fit and is used to estimate an output prediction range for exact inputs. However, as is normally the case, what about situations where the input isn't exact? How do you calculate a "prediction interval" that accounts for both the uncertainty in the model and the uncertainty in the input data? For example, Assume the model is: `y = a0 + (a1 * x1) + (a2 * x2)` Where `y`, `x1`, and `x2` are time series vectors. I only have one observation of the `y` time series with its associated `x1` and `x2` time series. That data is used to fit the model. However, I also have 1000 additional observations for `x1` and `x2`. I can easily calculate individual "prediction intervals", using the model and each `x1` `x2` pair, however I want to estimate a "prediction interval" that allows for all `x1` and `x2` observations. To be more specific, I generated the example below. The questionable code is below the `#======` comment line. Basically, I used `lm(..)` to fit the `y` vector to its associated `x1` and `x2` vectors. Next, I used `predict.lm(...interval="prediction", level=alpha)` to generate its typical `fit`, `lwr`, and `upr` vectors for each additional `x1` `x2` pair (see Graph 1, `fit` is the solid line, `upr` and `lwr` are the dashed lines for two observations). I then collected these vectors for all `x1` `x2` pairs and used `alpha` to extract the upper and lower tails of all of the `lwr` and `upr` vectors (see Graph 2). Is this scheme valid? Does the `alpha` that was used in `predict.lm(... level=alpha)` apply directly to counting-up/sorting the results to generate a "prediction interval" that allows for all `x1` `x2` pairs? Can a 5% range from a competing model (for example an ARIMA model) be compared to this 5% "prediction interval"? I'm fairly sure that the following scheme isn't right, but so far, I haven't figured out what I need to fix.               set.seed(1)          numpoi <- 10 #Number of data points in a time series vector          #First independent variable "x1", first observation     x1mea <- 0.03     x1sta <- 0.05     x1 <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))          #Second independent variable "x2", first observation     x2mea <- -0.01     x2sta <- 0.1     x2 <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))          #Dependent variable "y", first observation     a0 <- 3     a1 <- 2     a2 <- 1     noimea <- 0.0     noista <- 0.1     y <- a0 + (a1 * x1) + (a2 * x2) + rnorm(numpoi, mean=noimea, sd=noista)          #Build a data frame of the "first observation" data     datfra <- data.frame(y=y, x1=x1, x2=x2)          #Fit the model for the "first observation"     mod <- lm(y ~ x1 + x2, data=datfra)     summary(mod)          #Set up desired alpha value     alpha <- 0.95     onetai <- (1 - alpha)/2 #Convert the two tail "alpha" to a one tail value for use later          #Generate some new data "a" for a second observation of "x1" and "x2".     x1a <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))     x2a <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))     datfraa <- data.frame(y=rep(NA, numpoi), x1=x1a, x2=x2a)     modprea <- predict(mod, newdata=datfraa, interval="prediction", level=alpha)          #Generate some new data "b" for a third observation of "x1" and "x2".     x1b <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))     x2b <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))     datfrab <- data.frame(y=rep(NA, numpoi), x1=x1b, x2=x2b)     modpreb <- predict(mod, newdata=datfrab, interval="prediction", level=alpha)          #Plot the results for both new data "a" and "b"     plot(modprea[, 1], type="l", ylim=c(min(modprea, modpreb), max(modprea, modpreb)), main="Graph 1 - Second and Third Observations for x1 and x2", lwd=2, col="red")     lines(modprea[, 2], lwd=2, lty=2, col="red")     lines(modprea[, 3], lwd=2, lty=2, col="red")     lines(modpreb[, 1], lwd=2,  col="green")     lines(modpreb[, 2], lwd=2, lty=2, col="green")     lines(modpreb[, 3], lwd=2, lty=2, col="green")          #===========================================================================     #The code below is where my question lies.    Is this the appropriate method     #to account for all observations?          #Run the above calculation scheme on "all" observations.     numtri <- 1000 #All observations     modprecfit <- matrix(0, nrow=numpoi, ncol = numtri) #Matrix to hold all "fit" vectors     modpreclwr <- matrix(0, nrow=numpoi, ncol = numtri) #Matrix to hold all "lwr" vectors     modprecupr <- matrix(0, nrow=numpoi, ncol = numtri) #Matrix to hold all "upr" vectors          #Fill up the modprecfit, modpreclwr, and modprecupr matricies.     for (i in 1:numtri) {            #Generate the new data and put it in a dataframe       x1c <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))       x2c <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))       datfrac <- data.frame(y=rep(NA, numpoi), x1=x1c, x2=x2c)            #Predict the new "y" for new input data       modprec <- predict(mod, newdata=datfrac, interval="prediction", level=alpha)            #Store the "ft", "lwr", and "upr" vectors so they can be processed later       modprecfit[, i] <- modprec[, 1]       modpreclwr[, i] <- modprec[, 2]       modprecupr[, i] <- modprec[, 3]          }          #Extract the average "fit", the lower quantile "lwr", and the upper quantile "upr"     modprecfitfin <- apply(modprecfit, 1, quantile, 0.5)     modpreclwrfin <- apply(modpreclwr, 1, quantile, onetai)     modprecuprfin <- apply(modprecupr, 1, quantile, (1 - onetai))          plot(modprecfitfin, type="l", ylim=c(min(modpreclwrfin), max(modprecuprfin)), main="Graph 2 - All Observations for x1 and x2", lwd=2, col="red")     lines(modpreclwrfin, lwd=2, lty=2, col="red")     lines(modprecuprfin, lwd=2, lty=2, col="red")      ![enter image description here](http://i.stack.imgur.com/dHaNa.png) ![enter image description here](http://i.stack.imgur.com/oMcGL.png)


title: Output of one-tailed Wilcoxon sign rank test in R

text: I'm trying to assess the difference between two algorithms. They are stochastic, so I've run them multiple times against the same input files and noted their results. I want to determine whether my algorithm offers a statistically significant improvement over a previous approach. From my (potentially flawed) understanding this sort of situation is where I should (or rather could) use the Wilcoxon Sign Rank test. The data is naturally paired, so I can join together outputs of both algorithms by the particular input problem they used. Since I want to test that my algorithm is not just different, but actually better (in this case returning lower values) than the other algorithm this appears to be a one- tailed assumption. I'm using the R `wilcox.test` function to perform the test and I'm slightly confused about how I should interpret the results. I have read the help page for `wilcox.test` and it doesn't seem to offer much information about the results, more focus is made to the function's arguments. I've made a minimal working example with a small subset of my data:               x <- structure(list(instance = structure(c(1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L), .Label = c("competition01", "competition02",      "competition03", "competition04", "competition05", "competition06",      "competition07", "competition08", "competition09", "competition10",      "competition11", "competition12", "competition13", "competition14",      "competition15", "competition16", "competition17", "competition18",      "competition19", "competition20"), class = "factor"), nhoods = structure(c(1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("nhoods1", "nhoods2",      "nhoods3", "nhoods4", "nhoods5"), class = "factor"), run.no = structure(c(1L,      1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L), .Label = c("1", "2", "3",      "4", "5"), class = "factor"), partition = structure(c(1L, 2L,      1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("Hard-Soft", "Full"     ), class = "factor"), VNS = c(984L, 1445L, 1033L, 1445L, 1035L,      1318L, 1058L, 1445L, 913L, 1445L), `VNS-Skip` = c(1083L, 1425L,      1099L, 1230L, 1077L, 1363L, 1102L, 1442L, 1093L, 1252L)), .Names = c("instance",      "nhoods", "run.no", "partition", "VNS", "VNS-Skip"), row.names = c(NA,      10L), class = "data.frame")          wilcox.test(x[,5], x[,6], paired = TRUE, conf.int = TRUE, alternative = "greater")      My algorithm's results are in the 6th column and the original ones are in the 5th column. Since I want to assess whether my algorithm is better I've used the alternative = "greater" option which should mean that the test is checking for 1st arg > 2nd arg. This results in the output:                 Wilcoxon signed rank test            data:  x[, 5] and x[, 6]        V = 22, p-value = 0.7217       alternative hypothesis: true location shift is greater than 0        95 percent confidence interval:       -72 Inf        sample estimates:       (pseudo)median              -21       In this case the p-value is not less than 0.05 so there is not enough evidence to discard the null hypothesis. What does having infinity as the upper bound of a confidence interval mean? Is this because I'm using the one-tailed version of the test? All the tests I've run on my data always have the upper confidence interval as Inf. If the p-value were less than 0.05 would that mean I would be justified in saying "with a 95% confidence x[,5]'s mean will be within -72 of x[,6]'s?" What does the V value mean with regard to my data? From what I can see it is the difference between `median(x[,5])` and `median(x[,6]` but how would describe that in prose? Does anyone actually use the V value is the discussion of their data analyses?


title: Classifiers with post-training constraints on the prediction space

text: I'm familiar with using tools like SVMs and decision trees for discrete classification problems. But one detail that I have not encountered in that domain is: what do you do if your classifier must conserve some quantity when it is predicting outcomes for new input data? For example, suppose that you took all teams in the NFL and broke them into quintiles (20% buckets) based on their rushing yards per game (or say you took teams from the EPL and broke them down into quintiles based on their time-of- possession or some other statistic). Further let's say you are doing this on a weekly basis, so every week you collect your stats and then rank the teams 1 to N, chop them into quintiles, and each team's "class label" is its quintile. Now let's say you build a bunch of features for each team that you think will predict the team's quintile for next week. It could be the team's historical quintile, coaching changes, injuries to key players, whatever. But the point is you can make a feature vector for each team and you want to train a classifier that will predict which class label (1 to 5) to assign based on the features. What are some common or empirically effective ways to constrain the classifier so that the output makes sense? That is to say, it should only be possible to classify 20% of the teams in any particular quintile, and the classifier should not be able to logically contradict itself by putting more teams in a quintile than 20% of the total. (That is, the theory behind the classifier ought to demonstrate how we can _know in advance_ that the classifier won't make contradictions). I'm not looking for post-hoc hacky ways to do rounding or anything, and I am also specifically trying to avoid the situation where you build a regression model of the underlying continuous quantity and then just predict the quantity and manually assign quintiles yourself. I want the classifier to discretely choose quintiles for each feature vector while respecting external constraints on the class labels. More generically you could wonder, how do you train an SVM-like classifier such that it has internal rules like "Only M of the N objects can receive label L, and I need to take that into account all the time for _every other_ label that I assign too..." **Some thoughts:** (1) Use something like simulated annealing to shuffle the boundary points for the quintile after prediction. That is, use the trained classifier to classify every team's quintile. Use the distance from the margin as the "energy function" to be minimized, and then start randomly shuffling members of over- populated quintiles to nearby under-populated quintiles. You either accept or reject such a shuffling based on the extent to which it worsens the overall distance from the margin. I don't like this solution because it probably would be computationally inefficient after you already solve the SVM optimization problem. Also, it feels like you're actually just using the margin distance as a continuous quantity to determine the classification, at which point you might as well try to derive a more model-driven or physically plausible continuous quantity and just do regression or something and chop the predicted output into quintiles at the end, which is something I wanted to avoid.


title: Imputation of missing response variables

text: I am doing multiple imputation on a database of observations on hospital patients. There is one observation of many covariates per patient. There are 2 binary outcome variables:   1. Alive/Dead after 30 days   2. Died in hospital, or survived/discharged Two seperate analysis models (logistic regression), with identical covariates, are to be run, each with one of these outcomes as the response. There are 7 covariates in the analysis models. There is missing data of between 3 and 11% in the covariates and the outcomes. In addition there are 7 further covariates that are to be used in the imputation model to predict the missingness in the covariates and outcomes. My questions concerns the imputation of the two outcome variables. They are to be used as predictors for missingness in the covariates, as per standard practice, but they are highly collinear with each other. Is this a concern for the imputation model ? Is it valid/recommended to impute them both in the same model (to generate several complete datasets all containing both outcomes) or should seperate imputations be performed for each of the outcomes (to generate two distinct sets of complete datasets, each distinct set having one of the outcomes) ? Any other suggestions for how to proceed would be welcome.


title: Visual display of multiple comparisons test

text: Suppose, the data below shows the mean response time on a task for respondents among four different groups:               A     B     C    D       1.2   2.3   4.5  6.7      In order to assess which one of the means are different from one another I do a multiple comparisons test (after an omnibus ANOVA test is cleared) and the multiple comparisons test tells me that the mean for group D is significantly different from the ones for groups A and B and no other pair of differences is significantly different. What is the best way to present this information visually?


title: Building an input matrix in R for fisher's exact test; arrangement of variables

text: I am trying to implement a test as described in a paper in R. The aim is to calculate enrichment for a particular type of mutation that occurs in a given base in a given sequence context. I quote > Statistical evaluation of the over-representation of APOBEC signature > mutations in each sample was performed using a one-sided Fisher's exact test > comparing the ratio of the number of cytosine-to-thymine or cytosine-to- > guanine substitutions and guanine-to-adenine or guanine-to-cytosine > substitutions that occurred in and out of the APOBEC target motif (TCW or > WGA) to an analogous ratio for all cytosines and guanines that reside inside > and outside of the TCW or WGA motif within a sample fraction of the genome. Reference - http://www.nature.com/ng/journal/v45/n9/full/ng.2702.html This means they calculate enrichment by comparing CtoT or CtoG mutations in the base context TCW to those outside that context to the number of Cytosines found in and out of the TCW context in the genome. So we've got the following values.               TCW.mutations = the number of C to T or G mutations in the context TCW     non.TCW.CtoTG.mutations = the number of C to T or G mutations outside the context TCW     TCW.cytosines = the number of cytosines in the base context TCW in the human exome     nonTCW.cytosines = the number of cytosines outside the base context TCW.      I need to do a one-tailed fisher's exact test comparing the following using the fisher.test function in R               TCW.mutations/non.TCW.CtoTG.mutations to TCW.cytosines/nonTCW.cytosines      How is my input matrix supposed to look?


title: Multidimensional scaling using Python

text: I have 6,000 points for which I have all pairwise distances in a distance matrix. I want to get an idea whether these data were generated by a mixture of Gaussian distributions so I'm trying to get a visualization. I am trying to apply multidimensional scaling in 2 dimensions using sklearn in Python. I have four questions:   1. In all the examples I found, the number of points were very limited (around 20). Is it possible to apply MDS with 6,000 points?   2. If not, is there any way I can get a better idea whether my data were generated by a mixture of Gaussians? (I want to cluster this data using GMM.)   3. sklearn uses the function manifold.MDS. One of the parameters of this function is 'random state'. What does it represent?   4. In this example, the matrix used (X_true) has dimension 20x2. Shouldn't MDS take a distance matrix (i.e., a square matrix)?


title: Reproducing the Table 3.2 from the Elements of Statistical Learning

text: I was able to reproduce table 3.1 from ESL. However, when I tried to reproduce table 3.2, my estimated coefficients were way off (shown below):                        [,1]      [1,]  0.4292      [2,]  0.5765      [3,]  0.6140      [4,] -0.0190      [5,]  0.1448      [6,]  0.7372      [7,] -0.2063      [8,] -0.0295      [9,]  0.0095      The results from table 3.2 in ESL are as follow: ![enter image description here](http://i.stack.imgur.com/epSz4.png) I have attached my code here:               res <- read.table("prostate.data", sep = "")     fix(res)     XTraining = subset(res, train)     XTesting = subset(res, train == FALSE)               nrow = dim( XTraining )[1]     p = dim( XTraining )[2] - 1 # the last column is the response           D = XTraining[,1:p-1] # get the predictor data          # This gives the raw data for Table 3.1 from the book:     #      print(cor(D),digts=3)          library(xtable)     xtable( cor(D), caption="Duplication of the values from Table 3.1 from the book", digits=3 )          #     # Duplicate Table 3.2 from the book      #          # Append a column of ones:     #      Dp = cbind( matrix(1,nrow,1), as.matrix( D ) )          lpsa = XTraining[,p]          library(MASS)          betaHat = ginv( t(Dp) %*% Dp ) %*% t(Dp) %*% as.matrix(lpsa)          # this is basically the first column in Table 3.2:     #     print('first column: beta estimates')     print(betaHat,digits=2)      Does anyone have similar issues ? I would like to know what went wrong with my output.


title: Likelihood ratio test - lmer R - Non-nested models

text: I am currently reviewing some work and have come across the following, which seems wrong to me. Two mixed models are fitted (in R) using lmer. The models are non-nested and are compared by likelihood-ratio tests. In short, here is a reproducible example of what I have:               set.seed(105)     Resp = rnorm(100)     A = factor(rep(1:5,each=20))     B = factor(rep(1:2,times=50))     C = rep(1:4, times=25)     m1 = lmer(Resp ~ A + (1|C), REML = TRUE)     m2 = lmer(Resp ~ B + (1|C), REML = TRUE)     anova(m1,m2)      As far as I can see, `lmer` is used to compute the log-likelihood and the `anova` statement tests the difference between the models using a chi-square with the usual degrees of freedom. This does not seem correct to me. If it is correct, does anyone know of any reference justifying this? I am aware of methods relying on simulations (Paper by Lewis et al., 2011) and the approach developed by Vuong (1989) but I do not think that this is what is produced here. I do not think that the use of the `anova` statement is correct.


title: Determining sample size with a proportion and binomial distribution

text: I am trying to learn some statistics using the book, Biometry by Sokal and Rohlf (3e). This is an exercise in the 5th chapter which covers probability, the binomial distribution, and Poisson distribution. ![enter image description here](http://i.stack.imgur.com/T0Tth.jpg) I realize there is a formula to produce an answer to this question: $$ n = \frac 4 {( \sqrt{p} - \sqrt{q} )^2} $$ However, this equation is not in this text. I'd like to know how to calculate sample size knowing only the probability, the desired level of confidence, and the binomial distribution. Are there any resources covering this topic that I can be pointed to? I've tried Google, but what I've seen so far requires information I don't have access to in this problem.


title: Probability of intersection from multiple sampling of the same population

text: Here is an example case:   * I have a population of 10,000 items. Each item has an unique id.    * I randomly pick 100 items and record down the ids    * I put the 100 items back into the population    * I randomly pick 100 items again, record down the ids and replace.    * In total, I repeat this random sampling 5 times What is the probability that $X$ number of items appear in all 5 random samplings? I am not very well versed in statistics. Would this be correct for $X = 10$?   * For each sampling, the number of possible combinations of 100 items from 10,000 is ${\rm binom}(10000, 100)$   * Out of all possible combinations of 100 items, ${\rm binom}(9990, 90) * {\rm binom}(100, 10)$ combinations contain 10 specific items   * The probability of having 10 specific items is $({\rm binom}(9990, 90) * {\rm binom}(100, 10)) / {\rm binom}(10000, 100)$   * The calculated probability to the power of 5 would represent 5 indepenent samplings.  So essentially we are just calculating 5 independent hypergeometric probabilities and then multiplying them together? I feel like I am missing a step somewhere.


title: Continuous time Markov chain forward equation

text: This is a homework question and I need suggestion how to approach it. We have given the transitions   1. $\ i\rightarrow i+1$ with rate $\lambda(i)$ where $\ i \ge 1$    2. $\ i\rightarrow i-1$ with rate $\mu(i)(i-1)$ where $\ i \ge 2$ I am starting the forward equation like this: $$\ p_j(t) = [1-(\lambda_jh+\mu_jh +o(h))p(t)]+ \lambda(j-1)p(j-1)(t)h + ...$$ and I cannot really continue from here. I need to get to the point where I can show that $$\ G(z, t) = \sum P(N(t)= j|N(0) = a) z^j $$ satisfies $$\frac{\partial G}{\partial t} = z \left( z-1\right) \left(\lambda\frac{\partial G}{\partial z}- \mu\frac{\partial^2 G}{\partial z^2}\right)$$


title: How to generate predictions with rjags?

text: I've used rjags to run MCMC on a model, specified in the JAGS language. Is there a good way to extract that model and perform predictions with it (using the posterior distributions of my parameters)? I can re-specify the model in R and plug in the modes of my parameter posteriors; I'm just wondering if there's a less redundant way of doing this. I believe http://sourceforge.net/p/mcmc-jags/discussion/610037/thread/0ecab41c is asking the same question.


title: What is the meaning of a confidence interval taken from bootstrapped resamples?

text: I've been looking at numerous questions on this site regarding bootstrapping and confidence intervals, but I'm still confused. Part of the reason for my confusion is probably that I'm not advanced enough in my statistics knowledge to understand a lot of the answers. I'm about half-way through an introductory statistics course and my math level is only about mid-Algebra II, so anything past that level just confuses me. If one of the knowledgeable people on this site could explain this issue at my level it would be extremely helpful. We were learning in class how to take resamples using the bootstrap method and use those to build up a confidence interval for some statistic we'd like to measure. So for example, say we take a sample from a large population and find that 40% say they'll vote for Candidate A. We assume that this sample is a pretty accurate reflection of the original population, in which case we can take resamples from it to discover something about the population. So we take resamples and find (using a 95% confidence level) that the resulting confidence interval ranges from 35% to 45%. My question is, what does this confidence interval actually _mean_? I keep reading that there's a difference between (Frequentist) Confidence Intervals and (Bayesian) Credible Intervals. If I understood correctly, a credible interval would say that there's a 95% chance that in _our situation_ the true parameter is within the given interval (35%-45%), while a confidence interval would say that there's a 95% that in _this type of situation_ (but not necessarily in our situation specifically) the method we're using would accurately report that the true parameter is within the given interval. Assuming this definition is correct, my question is: What's the "true parameter" that we're talking about when using confidence intervals built up using the bootstrap method? Are we referring to (a) the true parameter of the _original population_ , or (b) the true parameter of the _sample_? If (a), then we'd be saying that 95% of the time the bootstrap method will accurately report true statements about the original population. But how could we possibly know that? Doesn't the whole bootstrap method rest on the _assumption_ that the original sample is an accurate reflection of the population it was taken from? If (b) then I don't understand the meaning of the confidence interval at all. Don't we already know the true parameter of the sample? It's a straightforward measurement! I discussed this with my teacher and she was quite helpful. But I'm still confused.


title: Expected value equal zero in Chi-square ($\chi^2$) test

text: I'm working on my final project, and it's a QUEST algorithm, and it uses the chi-square test. I have a problem where the expected value = 0.   1. What should I do?    2. What does it imply to the chi-square test or to the QUEST algorithm?


title: Estimating population size and a proportion

text: I am interested in the following sampling problem, which I will try to describe by a motivating example. Suppose we want to estimate how many people in a certain area, has blue eyes, how many have them has brown eyes etc.(think of a color scale) We do have an estimate on total number eye colors that the population contains (let's say 100 colors). We then sequentially observe samples - in a motivating context let's say a cruise arrives from that region each month. We disrupt those passengers and examine the eye color of each of them. The cruise does not always contain same number of passengers, and we don't know the decision rule on putting those people on that cruise. A reasonable assumption is that each person has an equal chance (p) of being in that cruise and reasonably we tend to see more the people with common eye colors. The problem is   1. Can we estimate total number of people living on the area?   2. Can we estimate the proportion or chance "p" ?


title: Building a forecast model based on past year data in R

text: I am attempting to build a model to forecast attendance in a given week in the current year based on this year's attendance values up until the present, and data from two previous years. My data looks like this:                  Week 11-12 Cumulative ADA    12-13 Cumulative ADA    13-14 Cumulative ADA        1    0.9941                  0.9941                  0.9914        2    0.9907                  0.991                   0.989        3    0.9888                  0.9888                  0.9879        4    0.9877                  0.987                   0.9869        5    0.9869                  0.9865                  0.9867        6    0.9862                  0.985                   0.9859        7    0.9856                  0.9842                  0.9857        8    0.9856                  0.984                   NA        9    0.9852                  0.9839                  NA        10   0.9848                  0.9834                  NA      Any guidance on how to predict the three NAs based on the past two years data and this year's values would be much appreciated. Thanks!


title: How to do linear mixed effects modeling with this data?

text: **Study Design:**   Below is a clinical trial in a longitudinal dataset. All subjects (n=34) attended `V1` (baseline) and then they were assigned to either a `Placebo` (n=15) or a `Drug` (n=19), both then were investigated at visit `V2`. Based on an investigation's result at `V2` called `grpVar1` they were grouped into two groups: `grpVar1=P` (n=16) and `grpVar1=N` (n=18).   Those `N` were routed for visit `V3`, whereas those `P` were routed for visit `V4` and later visit `V5`.   `Par1` was measured in all of the visits for all subjects at two anatomical sites: `site1` and `site2`. So each subject has two measurements at the two sites in all of the visits. Additionally, `Days` variable is the number of days of subsequent visits based on the baseline. Even though the visit is the same, different subjects attended the visit at different days and that was recorded. **It looks like this:**               > head(df,12)        id id_site Gender Visit Days grpVar1 Therapy        Par1     1   1   site1      F    V1    0       P Placebo  0.61030187     2   1   site2      F    V1    0       P Placebo  0.43833820     3   1   site1      F    V2  244       P Placebo  7.45492517     4   1   site2      F    V2  244       P Placebo  0.63105623     5   1   site1      F    V3  322       P Placebo  0.22600916     6   1   site2      F    V3  322       P Placebo  0.23180883     7   2   site1      M    V1    0       P    Drug  0.18130448     8   2   site2      M    V1    0       P    Drug  0.14165622     9   2   site1      M    V2  223       P    Drug  1.02354748     10  2   site2      M    V2  223       P    Drug  2.05987466     11  2   site1      M    V3  314       P    Drug  0.10194901     12  2   site2      M    V3  314       P    Drug  0.22315864      **Class of variables**               > sapply(df, class)            id   id_site    Gender     Visit      Days   grpVar1   Therapy      Par1       "factor"  "factor"  "factor"  "factor" "numeric"  "factor"  "factor" "numeric"       **Conception:**   So as you can see, repeated measures design with a twist; there is some kind of routing of subjects based on `V2` results, I call it "bifurcation of subjects" into two different routes of visits after `V2`.   Clearly there should be some individual variations that need to be accounted for by allowing both intercept and slope to change. Another source of variation I believe is the anatomical site that need to be accounted for alike. So `Par1` is measured in a hierarchy of site and subjects. Both `id_site` and `id` can be looked at as contextual variables; therefore, they can be treated as the random effects in a mixed model.   Although there were 5 visits as shown in `Visit` variable, but can we make use of the related `Days` variable to account for variation in the time as days after the baseline. I wonder how one can let the model account for this also? **My trial in R**               library(nlme)     lm_form_1 <- formula (Par1 ~ Visit + Gender + Therapy + id_site + grpVar1)     mod1 <- lme(lm_form_1, random = ~ 1|id , data = df)      **Question:**   How to account for variation of `id`, `id_site`, and `Days` as random effects. How to set the var-covar structure of `pdMat` class based on this study design. Is the layout of the dataset fit for modeling with `lme()` function or should be modified, what would be suitable? I would be very grateful to suggestions by the statistics community. Please limit your answer to the `nlme` R package and not `lme4`.


title: Regression using SVMs

text: I have set of time series data tuples: {$(1,22), (2, 25), (3, 18), (4, 26), ...$ so on} I want to estimate $Y_{t+1}$ using support vector machine regression technique. I have found numerous publications that mostly refer to Smola's paper. I want to understand the technique to implement it in my programming. There are numerous examples available for simple regression techniques, but is there any simple example available that solves the problem of finding regression with support vector machine on a given data set using the equations and also shows the steps?


title: How to apply t-test on MATLAB?

text: I have accuracy vectors of two classifiers. How can I make a t-test for them on Matlab and is there a good plot type to show it?


title: How to compute the power of Wilcoxon test?

text: I am using Wilcoxon test to compare two paired sets of data for whether their means differ. Besides the p-value, I would also like to know the power of this test. How to compute it in R? Thank you in advance.


title: How to assess whether a correlation between play and number of friends holds in different contexts?

text: I am trying to ascertain if a relationship exists between **hours of play** and **number of friends**. I have two contexts for the study: "at school" and "at a place other than the school". My data set relates to the "at school" context and the correlation r=0.8766 (significant). Now the difficult part.... As I want to focus on the **same group of students** for the "at a place other than the school" context, I have asked a question that filters my respondents and puts them in the "at a place other than the school context". The question is: Did you make any friends in the last three weeks? The answer can be either "yes" or "no". If it is yes, then I have asked them to indicate whether they made the friends "at the school", "at a place other than the school" or "at both places". So I have four data sets as follows:   1. Those who did not make any friends   2. Those who made one or more friends at the school   3. Those who made one or more friends at a place other than the school   4. Those who made one or more friends at both places (i.e. school and a place other than the school) I have ignored data set 1 because there is obviously no relationship between hours of play and number of friends because no new friends are made. I have also ignored data set 2 because I already have this finding (in the correlation above, which was at the school) I have also ignored data set 4 because it does not distinguish between "at school" and "at a place other than the school". I have used data set 3 and worked out the correlation, which is 0.7687 (significant). My interpretation of this finding is the relationship between hours of play and number of friends exists at the school and also outside the school. Does this make sense?


title: Can paired samples be treated as un-paired for 1-way ANOVA

text: This may be too simple a question for this forum - it is possibly a very basic statistics question as commonly found in biomedical research. But since scientists often lack good understanding of statistics, I would appreciate a lot if I could hear the thoughts from expert statisticians on this forum and be educated a little bit. I want to test if the treatment X affects the parameter P in a sample group of subjects. The design is paired: each subject is tested with no treatment (negative control) and with treatment X (order is randomized). To analyze such data, a paired t-test is obvious. However, I also have a positive control group - which is a separate independent group of subjects with sort of intrinsically different level of P, the one we would expect if treatment X works. So now I have three groups of measurements - two are paired, and one is un-paired. For three groups, I would need to use 1-way ANOVA. However, this will treat three groups as independent. How can I test if there is a difference in means between these three groups. Which test to use? Is it always appropriate to treat paired samples as un-paired? I understand that this leads to lower power, however - is this the only problem with this approach and is it still correct? Thank you!


title: Principal components analysis: how independent should the input measures be?

text: I have a waveform, for which most people will measure either a peak amplitude or a slope. I have included the area under the curve and several other measures involving the same original waveform. I now have 10 measures on such waveforms (and I have measured these waveforms in 500 patients). Is it legitimate to put all these measures side by side in one big principal component analysis (PCA) to distil a smaller set of principal components. My worry is about the fact that some of the measures are inherently interdependent (for example, peak amplitude will be inherently correlated with area under the curve), so I was afraid it might fail to meet some of the assumptions behind PCA, such as independence of measures. Is there another analysis method I should use instead of PCA? Independent Component Analysis? Another multi-dimensional scaling method? Cluster analysis? Thanks very much! A bit more detail: In my case, I have made 9 raw measures, _x1_ .. _x9_ , in each patient. I already know that _x2_ and _x4_ are highly correlated, but their difference has important prognostic implications. On one hand, I am afraid PCA would fail to note this difference (because _x2_ and _x4_ are so correlated), so I would like to include _x2-x4_ as a 10th measure. However, I am bothered by the fact that I am thus introducing trivial correlation, which will inflate the importance of _x2_ and _x4_ and any noise in them. Thanks again.


title: Detecting reflection in non-orthogonal rotation

text: I've known that, in orthogonal rotation, if the rotation matrix has determinant of -1 then reflection is present. Otherwise the determinant is +1 and we have pure rotation. May I extend this "sign-of-determinant" rule for non-orthogonal rotations? Such as orthogonal-into-oblique axes or oblique- into-orthogonal axes rotations? For example, this matrix                .9427   .2544   .1665   .1377     -.0451  -.0902  -.9940  -.0421      .3325   .3900   .1600   .8437      .4052   .8702   .2269   .1644      is an oblique-to-orthogonal rotation (I think, because sums of squares in rows, not in columns, are 1). Its determinant is -0.524. May I state that the rotation contains a reflection? Thanks in advance.


title: Multidimensional scaling of random matrix

text: Suppose I have a random symmetric matrix W of size $n\times n$, with i.i.d. coefficients uniformly distributed in [0,1], and I set $W_{ii} = 0$. Then I apply a Multidimensional Scaling of dimension $k$, which I define as minimizing the quantity: $H=n^{-2}\sum_{i,j} (W_{ij} - ||x_i-x_j||^2)^2$ with respect to the $n$ points $x_i \in R^k$. Say I denote $H^*$ the minimum. The question is : what is the behavior of $\mathbb{E}[H^*]$ as a function of $k$ and $n$ ? Is there a mathematical article out there dealing with this question ? or is it an open question ? Thanks.