title: Differences between lsmeans and difflsmeans

text: I've created an `lmer` model. The effects of one of my treatments, "Fertilizer" varies depending on whether I use `lsmeans` or `difflsmeans` from `lmerTest`. Which should I trust, `lsmeans` or `difflsmeans`? I have a split plot experiment. Fertilizer is the whole plot, Harvest is the split plot. Within each plot, there are 10 samples (quad). Data collected over 4 years, which I have made ordinal. The response is total biomass, which I have transformed.               > model <- lmer(logTotal ~ ordYear*Fertilizer*Harvest + (1|(Block:Fertilizer)) + (1|(Plot:Quad)), REML = FALSE, data = df)                   > difflsmeans(model, test.eff="Fertilizer")     Differences of LSMEANS:                                   Estimate Standard Error   DF t-value Lower CI Upper CI p-value        Fertilizer none-recommended       -0.3         0.0868 16.0   -3.27  -0.4677  -0.0996   0.005 **     Fertilizer none-half              -0.2         0.0868 16.0   -2.19  -0.3742  -0.0060   0.044 *      Fertilizer none-manure            -0.1         0.0868 16.0   -0.84  -0.2572   0.1110   0.412        Fertilizer recommended-half        0.1         0.0868 16.0    1.08  -0.0905   0.2777   0.297        Fertilizer recommended-manure      0.2         0.0868 16.0    2.42   0.0264   0.3946   0.028 *      Fertilizer half-manure             0.1         0.0868 16.0    1.35  -0.0671   0.3011   0.197        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1               > lsmeans(model, list(pairwise~Fertilizer))     NOTE: Results may be misleading due to involvement in interactions     $`lsmeans of Fertilizer`      Fertilizer    lsmean         SE    df lower.CL upper.CL      none        4.602075 0.07090356 21.33 4.454764 4.749387      recommended 4.885723 0.07090356 21.33 4.738411 5.033034      half        4.792150 0.07090356 21.33 4.644838 4.939462      manure      4.675189 0.07090356 21.33 4.527877 4.822500          Results are averaged over the levels of: ordYear, Harvest      Confidence level used: 0.95           $`pairwise differences of Fertilizer`      contrast                estimate        SE    df t.ratio p.value      none - recommended   -0.28364706 0.1002728 21.33  -2.829  0.0455      none - half          -0.19007464 0.1002728 21.33  -1.896  0.2595      none - manure        -0.07311311 0.1002728 21.33  -0.729  0.8843      recommended - half    0.09357242 0.1002728 21.33   0.933  0.7875      recommended - manure  0.21053396 0.1002728 21.33   2.100  0.1853      half - manure         0.11696153 0.1002728 21.33   1.166  0.6537          Results are averaged over the levels of: ordYear, Harvest      P value adjustment: tukey method for a family of 4 means       Thank you


title: How many counts for good frequency estimate (binomial)?

text: This is a pretty basic question, I think, but I'm finding it difficult to locate an answer or decent treatment of it on the web. For a binary dataset (binomial) with some unknown parameter p (probability of "heads" or success or something observed), how does the accuracy of the estimated p depend on the number of trials? I've read a number of ways to calculate confidence intervals around the estimated p, but these seem to all depend on assuming a large number of trials and a normal distribution. I'm interested in particular here with small numbers of trials. Thanks for any solutions or pointers.


title: How can I use the results of GARCH in order to improve a forecast?

text: I am kind of confused with what I should actually do with predicted volatility values that I obtained via a ARCH/GARCH model other than feeling happy that I know when volatility rises/falls. Is there a way that I can incorporate the predicted volatility values via a GARCH/ARCH model into a prediction model for my actual time series or is what I am saying erroneous? I use R as my primary tool.


title: How to learn the relationship between two variables

text: I have a data set where there are two variables X and Y, X being a set of values and Y being a set of corresponding values, like X = [1, 2, 3, 4] and Y = [100, 200, 300, 400]. Ideally I should have a directly proportional relationship. When X increases Y does too, and vice versa. What I would like to do is learn what might be the relationship, is it linear or non-linear or something else. How are such things explored? More generic questi![enter image description here on when give two such variables how do you start exploring it? How to check if there is correlation? Attaching the scatter plot image ![enter image description here](http://i.stack.imgur.com/MgMiC.png)


title: Which statistics to use?

text: I have a medical research in which I have two sets of data   1. Control group: People who have no ailment.   2. Study group: People who have some particular ailment. Within each category there are subgroups like age, gender, blood pressure ranges, cholesterol ranges, etc. Now I got a few observations (around 50) for both control group and study group. I have observation of various parameter like blood pressure, cholesterol level etc. I want to prove statistical significance of these parameter (Blood pressure/cholesterol) in study group. What is the best way to do it? Should i use $t$ distribution or normal distribution? Also when I calculate the $p$-value, should I compare each subgroup in study group (based on range) with corresponding subgroup in control group or calculate $p$-value as a whole for study vs control group?


title: Testing for Differences in Two Datasets

text: Assume that I have a dataset that contains the dates that measurements were made on a group of people. Thus, unique people identifiers, measurement results, and dates of measurements are available. Now assume that I have a another version of the same dataset from a different source, and that there are differences between the two datasets. Is there a statistical test that can be used to tell if the differences between the two datasets are statistically significant?


title: Stuck with interactions in logistic regression analysis

text: I would like to get some tips, because I am currently stuck with my regression part. I am using binary logistic regression. The aim is to compare two separate groups' responces. I have a MV (moderator variable, 1-one group, 0-other group). My IV (independent variables are 11 and all of them were measured on 7-point Likert scale). My DV (dependent variable or outcome is choice of 2 categories). First problem I face: I used the median split approach to obtain my DV for each respondent, as I do not have my DV, which I have to obtain to run regression. Second, I tried to run 2 seperate regressions for two groups, though my coefficients were poorly significant. Than I tried to run a one model with interactions. Calculating product term for each IV (IV1mean_centered* MVmean_centered) etc. for all IV I have. While running a regression I entered all the IV's vs. my MV and also calculated product terms with MV_centered (also tried without mean centering and with only MV). My issue is that my coefficients are even less significant than with seperate models!!! Does anyone have some suggestions how to solve this issue? Maybe I have made some mistake regarding interactions in logistic regression? Help!!!!


title: Validation error less than training error -- implications?

text: I am running a neural net to predict used car prices, sample size is 800. Using both 10-fold cross validation (10 times) and 1/3 holdback (10 times), the $R^2$ for training is about 0.60 and for validation is about 0.68 for all 20 runs. The smallest difference in the 20 runs is training $R^2$ = 0.64 and validation $R^2$ = 0.68, so the training $R^2$ is always less than the validation $R^2$. I am very used to seeing training $R^2$ bigger than validation $R^2$, which means overfitting. In the past when I have seen training $R^2$ less than validation $R^2$, it has been a transient phenomenon that disappeared when I re-ran the model. This is the first time that I have seen validation $R^2$ systematically larger than training $R^2$. I have no idea what this means. Any thoughts?


Can a percentage be lets say 0.172% even after multiplying it by 100? The reason is, in my calculation i did the following; (4.3601/2329.475495)x100 giving me the above answer. Is that correct or am i in the wrong?


title: Discussion about proxy- and instrument variables and endogeneity in the context of a multi equation model

text: Assume two equations   1. $Y_1 = X_1\beta_1 + X_2\beta_2 + U_1$   2. $Y_2 = X_1\alpha_1 + X_2\alpha_2 + U_2$ Further assume that $ \ U_1 = X_4 + E_1$ and $U_2 = X_4 + E_2$ with $ \ corr(Y_1,X_4)\ne 0, \ \ corr(Y_1,E_2)= 0 \ \ \ $ and $ \ \ \ corr(Y_2,X_4)\ne 0, \ \ corr(Y_2,E_1)= 0$ So if we would estimate eq. 1 and 2 by OLS we would get a omitted-variable- bias-problem since the regressor $X_4$ does not appear in either equation. **My question is:** Since $Y_1$ and $Y_2$ do not appear in those equations as explanatory variables, i.e. are not correlated with each other, then could one use $Y_1$ as a proxy for $X_4$ in eq. 2 and analogously $Y_2$ as a proxy for $X_4$ in eq. 1? My answer would be, yes if one controls for the potential endogenity of those variables since the two requirements for a proxy variable are given   1. the proxy variable is correlated with the unobservable variable   2. the proxy variable is not partially correlated with the dependend variable * * * Secondly assume the two equations   1. $Y_1 = Y_2\beta_1 + X_1\beta_2 + X_2\beta_3 + U_1$   2. $Y_2 = Y_1\alpha_1 + X_1\alpha_2 + X_2\alpha_3 + U_2$ Further assume that $ \ U_1 = X_4 + E_1$ and $U_2 = X_4 + E_2$ with $ \ corr(Y_1,X_4)\ne 0, \ \ corr(Y_1,E_2)= 0 \ \ \ $ and $ \ \ \ corr(Y_2,X_4)\ne 0, \ \ corr(Y_2,E_1)= 0$ So if we would estimate eq. 1 and 2 by OLS we would get a omitted-variable- bias-problem paired with a endogenity-problem since the regressor $X_4$ does not appear in either equation and the dependend variables $Y_1$ and $Y_2$ are jointly determined within the model. **My question is:** For instance: The endogenity of $Y_2$ in eq. 1 is due to   1. the partial correlation to $Y_1$   2. the partial correlation to the omitted regressor $X_4$ If I would use a control function approach which controls for this endogenity of $Y_2$ in eq. 1 dont I get a quantitiy which not only controls (1) but also for (2)? In my opinion this does not yield a solution for a omitted-variable- bias-problem but yield a solution for the endogenity? **So bottom line:** If $Y_1$ and $Y_2$ are actually not correlated I would actually be able to solve the omitted-variables-bias-problem. In any case I need to justify which of my problems is more serious in order to find a adequat solution.


title: Confusion related to linear dynamic systems

text: I was reading this book Pattern Recognition and Machine Learning by Bishop. I had a confusion related to a derivation of the linear dynamical system. In LDS we assume the latent variables to be continuous. If Z denotes the latent variables and X denotes the observed variables $p(z_n|z_{n-1}) = N(z_n|Az_{n-1},\tau)$ $p(x_n|z_n) = N(x_n,Cz_n,\Sigma)$ $p(z_1) = N(z_1|u_0,V_0)$ In LDS also alpha beta forward backward message passing is used to calculated the posterior latent distribution i.e $p(z_n|X)$ $\alpha(z_n)=p(x1...xn,z_n)$ $\hat\alpha(z_n) = \alpha(z_n)/P(x1....xn)$ My first question is in the book it is given as $\hat\alpha(z_n) = N(z_n|u_n,V_n)$ How come we got the above. I mean $\hat\alpha(z_n)$ = $N(z_n|u_n,V_n))$. I mean how did we get this? My next question is related to the derivation as you can follow along the screenshots of the pages of the book attached. I didn't get what where $K_n$ came from and what Kalman filter gain is $u_n = Au_{n-1} + K_n(x_n - CAu_{n-1})$ $V_n = I - K_nC)P_(n-1)$ $c_n = N(x_n|CAu_{n-1},CP_{n-1}C^T + \Sigma$ $K_n$ is the Kalman gain matrix $P_{n-1}C^T(CP_{n-1}C^T + \Sigma) ^ {-1}$ How did we derive the above equations, I mean how come $u_n = Au_{n-1} + K_n(x_n - CAu_{n-1})$ I am just confused how the above derivation is made. ![enter image description here](http://i.stack.imgur.com/1Uutt.png) ![enter image description here](http://i.stack.imgur.com/hBijX.png) ![enter image description here](http://i.stack.imgur.com/C0xWn.png)


title: How many parameters in this specific linear model with interaction?

text: I have a question where I am not sure about the answer: > A linear model has the following characteristics: > > *A dependent variable ($y$) > > *One continuous variable ($x_l$), including a quadratic term ($x_1^2$) > > *One categorical ($d$ with 3 levels) predictor variable and an interaction > term ($d \times x_1$) > > How many parameters, including the intercept, are associated with this > model? So I thought about: $y=\beta_0 + \beta_1 x_1 + \beta_2 x_1^2 + \beta_3 d$ + $\beta_4dx_1$ I thought I need 5 parameters for this model, but I was told the answer would be 7. So I guess I did the interaction wrong. So it seems to me I have to model each level of d with x1 with a single interaction. But how would I write this down and why do I need this?


title: Bayes' Theorem - Probability Pants problem

text: I'm having an issue with a question regarding Bayes' Theorem. Here is the question: > An online clothing store carries three brands of jeans. 40% of sales are > brand A, 20% are brand B and the remainder are brand C. 20% of brand A jeans > cost over 100, 40% of brand B jeans cost over $100 and 90% of brand C jeans > cost over 100. Given that a pair of jeans is purchased for over 100, what is > the probability that they are brand A? My work is as follows: Using Bayes Theorem I classified two events $A$: Picking brand A $B$: Jeans cost over 100 So I would need to find $P(A|B)$ $$P(A) = 2/5$$ $$P(B|A) = 1/5$$ $$P(\neg A) = 3/5$$ $$P(B|\neg A) = 11/25$$ When I do calculations with Bayes' theorem, I get; $$\frac{2/25}{2/25 + (11/25)\times(3/5)} = 10/43$$ But the answer is $2/13$. Now, when I first did the question, I forgot to multiply $P(B|\neg A)$ by $P(\neg A)$ in the denominator, and I got the right answer. Is there a reason I should leave $P(\neg A)$ out? Or did I approach the problem in the wrong way completely?


title: Weighted average of precision (or recall) of all classes?

text: How to compute the global precision given the precision calculated for each class? Is it just the average over classes precisions? When I use Weka the global precision is not computed as the average one, but as a "weighted average" and I don't know how this later is computed (weighted by what?). See the example bellow of the results that we get with Weka, can you please tell me how "Weighted Avg." is computed in this result (last line)?               Correctly Classified Instances         608               93.5385 %     Incorrectly Classified Instances        42                6.4615 %     Total Number of Instances              650               === Detailed Accuracy By Class ===                     TP Rate   FP Rate   Precision   Recall  F-Measure   ROC Area  Class                  0.889     0.005      0.727     0.889     0.8        0.995    APP05179028                  0.991     0.032      0.866     0.991     0.924      0.999    APP05179007                  0.633     0.005      0.864     0.633     0.731      0.989    APP05179012                  0         0          0         0         0          0.954    APP05179014                  0.972     0.013      0.936     0.972     0.954      0.998    APP05179010                  0.957     0.002      0.957     0.957     0.957      0.999    APP05179009                  0         0          0         0         0          ?        APP05179018                  0.6       0          1         0.6       0.75       0.988    APP05179027                  0         0          0         0         0          ?        APP05179023                  1         0.002      0.75      1         0.857      1        APP05179025                  0.8       0          1         0.8       0.889      0.991    APP05179016                  0         0          0         0         0          ?        APP05179021                  0.889     0          1         0.889     0.941      1        APP05179029                  0.918     0          1         0.918     0.957      1        APP05179020                  0         0          0         0         0          0.958    APP05179022                  1         0.002      0.963     1         0.981      0.999    APP05179015                  0.846     0          1         0.846     0.917      0.998    APP05179008                  0.992     0.002      0.992     0.992     0.992      1        APP05179011                  0.947     0.002      0.973     0.947     0.96       1        APP05179013                  1         0          1         1         1          1        APP05179026                  0.857     0          1         0.857     0.923      1        APP05179019                  0.714     0.011      0.417     0.714     0.526      0.952    APP05179017                  0.969     0          1         0.969     0.984      1        APP05179030                  1         0          1         1         1          1        APP05179024                  0.935     0.008      0.939     0.935     0.933      0.998    **Weighted Avg.**                       a   b   c   d   e   f   g   h   i   j   k   l   m   n   o   p   q   r   s   t   u   v   w   x         8   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 |   a = APP05179028        0 110   0   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 |   b = APP05179007        0   9  19   0   2   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 |   c = APP05179012        0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0   1   0   0   0   0   0 |   d = APP05179014        0   3   0   0 103   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 |   e = APP05179010        0   1   0   0   0  22   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 |   f = APP05179009        0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 |   g = APP05179018        3   0   0   0   0   0   0   6   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0 |   h = APP05179027        0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 |   i = APP05179023        0   0   0   0   0   0   0   0   0   3   0   0   0   0   0   0   0   0   0   0   0   0   0   0 |   j = APP05179025        0   1   0   0   0   0   0   0   0   0   4   0   0   0   0   0   0   0   0   0   0   0   0   0 |   k = APP05179016        0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 |   l = APP05179021        0   0   0   0   0   0   0   0   0   0   0   0  32   0   0   0   0   0   0   0   0   4   0   0 |   m = APP05179029        0   1   0   0   0   1   0   0   0   0   0   0   0  45   0   0   0   0   0   0   0   2   0   0 |   n = APP05179020        0   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 |   o = APP05179022        0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  26   0   0   0   0   0   0   0   0 |   p = APP05179015        0   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0  11   0   0   0   0   1   0   0 |   q = APP05179008        0   0   0   0   1   0   0   0   0   0   0   0   0   0   0   0   0 126   0   0   0   0   0   0 |   r = APP05179011        0   1   0   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0  36   0   0   0   0   0 |   s = APP05179013        0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   3   0   0   0   0 |   t = APP05179026        0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   6   0   0   0 |   u = APP05179019        0   0   0   0   2   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   5   0   0 |   v = APP05179017        0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0   0  31   0 |   w = APP05179030        0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  12 |   x = APP05179024


title: What is the requirement on the instrumental density in importance sampling?

text: 1. From Rubinstein's _Simulation Monte Carlo Method_ ![enter image description here](http://i.stack.imgur.com/hYF9R.png) I was wondering why $g$ is required to dominate $Hf$?   2. From Wikipedia > The basic idea of importance sampling is to change the probability $P$ so > that the estimation of $E[X;P]$ is easier. Choose a random variable $L\geq > 0$ such that $E[L;P]=1$ and that $P$-almost everywhere $L(\omega)\neq 0$. > The variate $L$ defines another probability $P^{(L)}=L\, P$ that satisfies > $$ \mathbf{E}[X;P] = \mathbf{E}\left[\frac{X}{L};P^{(L)}\right]. $$ The conditions on the instrumental density $L$ are $L\geq 0$ such that $E[L;P]=1$ and that $P$-almost everywhere $L(\omega)\neq 0$. I think I can understand them, because they are equivalent to say that $L$ is a density of $P^{(L)}$ wrt $P$, and the denominator in $\frac{X}{L}$ is zero only on a subset of probability measure zero so that it won't affect the expectation/integration. Given that the two sources seem to give different answers, I wonder what conditions should be put on the instrumental density? Thanks and regards!


title: Resampling or Basic Simulation and confidence intervals

text: I have a population of sales that might be won or lost. I know the rate that they are won from historical data. This case 30% of them historically win. To figure out how much money I will be making off of these sales I need to weight them. I could just add up all the values of the sales then take 30% of that number, but there is high variance in what the value of single sales are.               > describe(Dataset$V1)           var    n     mean      sd median  trimmed   mad  min     max    range  skew         1   1 2314 185090.7 1129744  32000 52077.27 37065 1000 3.5e+07 34999000 19.26           kurtosis       se         1    478.5 23485.42         > sum(Dataset$V1)*.3          [1] 128489970      So I figured that if you win 30% of them then let's try sampling the population at a win rate of 30%. I ran the following code               #Put the Percentage Value you want to use here     P <- .30     X <- round(nrow(Dataset)*P)     #Figured out how many wins will be based on different sample sizes.     boot <- data.frame(       boot = replicate(10000, sum(sample(Dataset$V1, size=X)))     )     write.csv(boot, file=ThePath, row.names=FALSE)      I have replacement as False because a single sale can't be won more than once. I do know that only the very first sale is picked at a 30% rate and the rest are picked at a slightly lower rate due to no replacement. So it randomly takes 30% of the data and then adds it up because I only care about the total expected value. I have it do this 10,000 times then I turned it into a histogram to see the distribution. Then have it set up to calculate what % of the 10,000 runs each bin has. Then look for about 8000 or 80% of the values in sequential bins. Can I then say that 80% of the time in a simulation of 10,000 the expected total value of the sales is expected to be within X and X? Or should I be using a normal bootstrap confidence interval calculation? Does the fact that the original data is skewed create issues and in some of the datasets when I run this code the Bootstrapped data comes out almost bimodal, does that lead to problems also?


title: What is the difference between pooled cross sectional data and panel data?

text: They seem so similar. Are they the same thing but just referred to as different names?


title: Confused about proper way to normalize two variables

text: I have two variables of interest:   * Residential Vacancies (res_vac)   * Commercial Vacancies (com_vac) I also have two variables with which I might normalize the above:   * Total Residences (res_tot)   * Total Businesses (bus_tot) Commercial Vacancies I get the general theory of how to normalize data, but the dichotomous split is giving my brain a headache. If I wanted to normalize res_vac (say, to compare it against a normalized com_vac variable), for example, which is correct:   * res_vac / res_tot   * res_vac / (res_tot + bus_tot) I'm leaning towards the first one. A great answer will also help me understand why it's one or the other.


title: Kullback–Leibler divergence between two Wishart distributions

text: The result is shown in: [1] W.D. Penny, KL-Divergences of Normal, Gamma, Dirichlet, and Wishart densities, Available at: www.fil.ion.ucl.ac.uk/~wpenny/publications/densities.ps But could anyone help me out to understand the lines on top of page 3: $L=\int \log(\Gamma) Wishart(\Gamma|a,B) d\Gamma$ $=\log(\tilde\Gamma(a,B))$ $=\sum\limits_{i=1}^{d}\Psi((a_{s}+1-i)/2)-\log|B_{s}|+d\log(2)$ ($\Psi(.)$ is digamma function) So what is $\tilde\Gamma(a,B),a_{s},B_{s}$?


title: Regarding the sampling procedure in Adaboost algorithm

text: The AdaBoost algorithm states that it is to train a classifier based on the training data according to a weight vector. Assume the size of training data is N, the weight vector is of dimension N as well. I have three questions regarding this sampling procedure, 1) Will the size of sampled data be the same as the original data set? 2) What does the weight vector look like? If it is a distribution, then the sum of them has to be 1. Is that possible to have a weight vector with entries of integer number? 3) Generally, which algorithm can be used to sample a data set based on a given weight vector or a distribution?


title: Do low silhouette widths mean the data has little underlying structure?

text: I am new to sequence analysis, and I was wondering how you react if the average silhouette widths (ASW) from cluster analyses of Optimal Matching- based dissimilarity matrices are low (around.25). Would it seem appropriate to conclude that there is little underlying structure that would allow the sequences to be clustered? Might you ignore the low ASW based on other measures of cluster quality (I have pasted some below)? Or is it likely that choices made during the sequence analysis or subsequent cluster analyses might be responsible for the low ASW numbers? Any suggestions would be appreciated. Thanks. In case more context is needed: I am examining 624 sequences of work hour mismatches (i.e., mismatches between the number of hours a person prefers to work in a week and the number of hours they actual work) among people in their 20s. All the sequences I am examining have a length of 10. My sequence object has five states (M=wants more hours, S=wants the same hours, F=wants fewer hours, O=out of the labor force, and U=unemployed). I have not done a systematic accounting of how ASW results vary with different combinations of approaches. Still, I have tried low and medium indel costs (.1 and .6 of the max substitution cost--I care more about the order of events than their timing) and different clustering procedures (ward, average, and pam). My overall impression is that the ASW numbers remain low. Perhaps low ASW results make sense. I would expect these states to come in a variety of different orders, and the states can be repeated. Removing duplicate observations only lowers the N from 624 to 536. Studying the data reveals that there is indeed a good bit of variety and sequences that I would consider very different e.g., people who wanted the same hours the entire time, developed a mismatch, resolved a mismatch, and oscillated back and forth between having and not having a mismatch. Perhaps lack of clearly differentiated clusters is not the same thing as a lack of interesting variation. Still, the weak cluster results seem to leave me without a nice way to summarize the sequences. Results from Ward's method with indel set at .1 of the substitution cost of 2 These statistics seem to suggest a 6 cluster solution might be good. The ASW, however, is low-- at least for solutions that have a reasonable number of clusters (2 or 3 is too few).                          PBC   HG HGSD  ASW ASWw     CH   R2   CHsq R2sq   HC     cluster2  0.56 0.78 0.75 0.38 0.38 110.76 0.15 241.65 0.28 0.14     cluster3  0.51 0.68 0.65 0.27 0.27 108.10 0.26 237.60 0.43 0.17     cluster4  0.54 0.74 0.71 0.25 0.25  88.66 0.30 203.72 0.50 0.14     cluster5  0.59 0.83 0.79 0.25 0.25  75.85 0.33 183.21 0.54 0.09     cluster6  0.59 0.85 0.82 0.24 0.25  66.94 0.35 164.51 0.57 0.08     cluster7  0.47 0.79 0.75 0.18 0.19  64.09 0.38 154.47 0.60 0.12     cluster8  0.47 0.81 0.77 0.20 0.21  59.47 0.40 152.36 0.63 0.11     cluster9  0.48 0.84 0.80 0.19 0.21  56.68 0.42 147.83 0.66 0.10     cluster10 0.47 0.86 0.82 0.19 0.21  53.24 0.44 140.18 0.67 0.08


title: Whether the samples from Latin hypercube method can be used for bootstrap?

text: This question is a simple version for my previous question: Choice of sampling method before using bootstrap validation method for computer models Can the samples from Latin hypercube (LHS) method be used for bootstrap? If so, the advantages of using LHS will lose since the new bootstrap samples would be non-LHS. It is possible to overestimate the variances of statistics interested using the combinations of LHS and bootstrap. Is the random sampling or quasi-random method appropriate to be used for bootstrap? Any thoughts would be greatly appreciated.


title: Conditional table, how to fill it

text: I found this question and since I'm learning probability I'm not sure how to go about it:               Of a company’s employees, 30% are      women and 6% are married women.      Suppose an employee is selected at      random. If the employee selected is a      woman, what is the probability that she is      married?      I have done it like this:                   W       W'          M  0.018          M' 0.282          p(W)=0.3  p(W')0.7          p(W')=1-p(W)=1-0.3=0.7     p(M and W)=p(W)p(M|W)=0.3x0.06=0.018     p(W and M')=p(W)-p(M and W')= 0.3-0.018=0.282      Well that's where I'm at, how do I go about it next? Can I get this numbers? Thanks.


title: How do I analyse data with 2 independent variables and 2 dependent variables?

text: Im confused about what to do. I was thinking of running two seperate multiple regressions with a DV in each.. its after this, that I'm stuck. How do I see what effect the two DV's have combined? Or am I going about the whole thing wrong. Help! My IV's are gender and group. My DV's will be scores on two seperate psychometric tests (likert scales) I hope to have at least 100 people in each group (3 groups) so sample size will be roughly 300.


title: Percent correctly predicted of logit model

text: Is there a standard way to report the percent correctly predicted when predicting a binary outcome? Using glm in r, the results are predicted probabilities. However, in order to make a comparison to another model, I want to report a single percent correctly predicted value from my binary model. Do I simply choose a cutpoint, and if so, how? Here is a simple example of the code.               model.results <- glm(binary.outcome ~ predictor1 + predictor2, family=quasibinomial)      Thanks,


title: Interpreting interactions between two treatments

text: I have two treatments A & B. Here are my groups, where X represents the appropriate control for that particular treatment: Group 1: XX Group 2: AX Group 3: XB Group 4: AB The hypothesis is that the treatment B will have an effect, but that that effect will no longer be apparent when combined with treatment A. So, if run my experiment and run an ANOVA on the data, and the results of the analysis show that only Group 3 was significantly different than the others, is it correct to say that "treatment B had an effect and that effect was lost when combined with treatment A"? Or, do would I also need to show a significant different between XB and AB?


title: Assessing significance of correlation

text: I have two variables, and I can calculate e.g. the Pearson correlation between them, but I would like to know something analogous to what a t-test would give me (i.e. some notion of how significant the correlation is). Does such a thing exist?


title: Estimating probability distribution function of a data stream

text: I have a very large number of observations. Observations arrive sequentially. Each observation is an $n$-dimensional vector (with $n \ge 100$), is independent from the others and is drawn from the same unknown distribution. Is there an optimal policy to estimate the unknown distribution, given some space bounds on the number of observations that can be stored? I would leave the estimation criteria open-ended, (in terms of expected or minimax error, asymptotic consistency, aysmptotic efficiency etc.).


title: Why is KNN not "model-based"?

text: ESL chapter 2.4 seems to classify linear regression as "model-based", because it assumes $f(x) \approx x\cdot\beta$, whereas no similar approximation is stated for k-nearest neighbors. But aren't both methods making assumptions about $f(x)$? Later on in 2.4 it even says: >   * Least squares assumes $f(x)$ is well approximated by a globally linear > function. >   * k-nearest neighbors assumes $f(x)$ is well approximated by a locally > constant function. > The KNN assumption seems like it could also be formalized (though not sure if doing so would _lead_ to the KNN algorithm in the way assuming $f$ is linear leads to linear regression). So if KNN actually isn't model-based, why? Or am I misreading ESL?


title: Combining auto.arima() and ets() from the forecast package

text: I've been using the ets() and auto.arima() functions from the forecast package to forecast a large number of univariate time series. I've been using the following function to choose between the 2 methods, but I was wondering if CrossValidated had any better (or less naive) ideas for automatic forecasting.               auto.ts <- function(x,ic="aic") {         XP=ets(x, ic=ic)          AR=auto.arima(x, ic=ic)              if (get(ic,AR)<get(ic,XP)) {             model<-AR         }         else {             model<-XP         }             model     }      /edit: What about this function?               auto.ts <- function(x,ic="aic",holdout=0) {         S<-start(x)[1]+(start(x)[2]-1)/frequency(x) #Convert YM vector to decimal year         E<-end(x)[1]+(end(x)[2]-1)/frequency(x)         holdout<-holdout/frequency(x) #Convert holdout in months to decimal year         fitperiod<-window(x,S,E-holdout) #Determine fit window              if (holdout==0) {             testperiod<-fitperiod         }         else {             testperiod<-window(x,E-holdout+1/frequency(x),E) #Determine test window         }              XP=ets(fitperiod, ic=ic)         AR=auto.arima(fitperiod, ic=ic)              if (holdout==0) {             AR_acc<-accuracy(AR)             XP_acc<-accuracy(XP)         }         else {             AR_acc<-accuracy(forecast(AR,holdout*frequency(x)),testperiod)             XP_acc<-accuracy(forecast(XP,holdout*frequency(x)),testperiod)         }              if (AR_acc[3]<XP_acc[3]) { #Use MAE             model<-AR         }         else {             model<-XP         }         model     }      The "holdout" is the number of periods you wish to use as an out of sample test. The function then calculates a fit window and a test window based on this parameter. Then it runs the auto.arima and ets functions on the fit window, and chooses the one with the lowest MAE in the test window. If the holdout is equal to 0, it tests the in-sample fit. Is there a way to automatically update the chosen model with the complete dataset, once it has been selected?


title: What is the reason for getting so many ”Warnings: NAs produced” while using the function ”pop.sim” in the ”demography” package

text: When so many warnings, what does it in fact means? Is there a problem with the validity of the stochastic forecasting models? I am doing a stochastic forecasting of a small population with approx. 50.000 people. I am using one year age intervals 0-90+. Because so small dataset, I am borrowing the mortality rates from a population which is similar in life expectancy. Fertility rates with 0 are replaced with 1/10000. Net migrations are calculated with the function ”netmigration” in the ”demography” package. Then I use the function ”pop.sim” for simulating say 1.000 sample paths of population 40 years ahead. When doing so with set.seed(505) and N=1000, I don't get any warning. But nearly all other values in set.seed and or N=1000 give me the warning messages: NAs produced. For example when my code is               set.seed(300)     sim300 <- pop.sim(mort=mort.fcast, fert=fert.fcast, mig=mig.fcast,                        firstyearpop=mort.fo,N=300, mfratio=mfratio,                       bootstrap=FALSE)      The first ten messages looks like this:               Warning messages:          1: In rpois(rep(1, length(fert$age)), lambda) : NAs produced           2: In rbinom(1, B, mfratio/(1 + mfratio)) : NAs produced           3: In rpois(1, Ef0 * mort.sim$female[1, j, i]) : NAs produced       4: In rpois(1, Em0 * mort.sim$male[1, j, i]) : NAs produced           5: In rpois(rep(1, p), Ef * mort.sim$female[, j, i]) : NAs produced       6: In rpois(rep(1, p), Em * mort.sim$male[, j, i]) : NAs produced           7: In rpois(1, Ef0 * mort.sim$female[1, j, i]) : NAs produced       8: In rpois(1, Em0 * mort.sim$male[1, j, i]) : NAs produced         9: In rpois(rep(1, p), Ef * mort.sim$female[, j, i]) : NAs produced     10: In rpois(rep(1, p), Em * mort.sim$male[, j, i]) : NAs produced        How shall I interpret these messages? Does this indicate that something is wrong in the data, models, forecasting or simulation procedures? My first guess is that, in message nr. 1, a negative value is assigned to ”lambda”?, if so, is there a way to prevent that? The most important question is whether these NAs produced is an indicator of the validity of the models and or the forecast of the population? Is anyone out there who can say something about what is going on here?


title: Where do the descriptors for Cronbach's alpha values come from (e.g., poor, excellent)?

text: It seems fairly common to describe Cronbach's alpha values as follows:   * α ≥ 0.9 Excellent   * 0.7 ≤ α < 0.9 Good   * 0.6 ≤ α < 0.7 Acceptable   * 0.5 ≤ α < 0.6 Poor   * α < 0.5 Unacceptable Where do these values come from? I cannot find an original research article describing these. Edit: I am 90% sure its merely based on convention and there is no classic research article outlining these.


title: Questions on bootstrap methods on calculating p-value

text: I have a silly question on bootstrap methods for calculating p-value. Suppose I have a dataset, which is likely normal distributed. Now I want to determine whether a particular value X is signifiant high in the dataset by giving a pvalue. For me, there are two options   1. since the dataset is roughly normal distributed, I could calculate the sd and mean, get the area for z>(X-mean)/sd from Z table, and use it as pvalue.   2. by bootstrap. I could generate a new sample S (length N) from the original dataset by bootstrap, then get the p-value by sum(S>X)/ N.  My question is am I right to use bootstrap in 2nd option? Could you recommend books for better understanding bootstrap and permutation test? Thanks.


title: Formal errors from non-negative least-squares?

text: I have seen a couple of similar questions, but none with answers, so perhaps I can formulate this question in a way that will be answerable. I am computing a standard linear regression subject to a positivity constraint using non-negative least squares (`lsqnonneg` in Matlab, actually). Is it possible to compute formal errors from a non-negative least squares problem, and if so how to do it? If no analytical form of the errors exists, I could just bootstrap, I suppose. For a regular problem I would use something like:               mcov = sigma^2*inv(G'*G);      where sigma is the data variance and G is the system matrix. But, this is very close to singular for my problem, and thus mcov blows up.


title: Comparing two sets of data

text: I have a research on the effectiveness of the new curriculum for our school. One of my indicators is the students' academic standing (passing/failure percentages). Available data for the new curriculum is only for a year, while I have data for the last five years for the old curriculum. I was planning to use the averages from the old curriculum and then compare it to the recent one. How can I determine whether the changes computed are statistically significant?


title: lme4 - extracting the diagonal of hat matrix

text: I'm trying to generate an autoplot for mer objects in the same vein as the the autoplot.lm example. I can extract the original data frame, the residuals and the linear predictors directly from the returned object...               > random.model <- lmer(a ~ b + c + (1 | d), data = example, family = binomial     > diagnostics  <- cbind(random.model@frame, random.model@eta, random.model@resid)      ...and after reading a suggestion here I can calculate Cooks Distance using the influence.ME package...               > library(influence.ME)     > cooks       <- cooks.distance(influence.ME::influence(random.model))     > diagnostics <- cbind(diagnostics, cooks)      ...and elsewhere (sorry can't find link) found that I could derive the standardised residuals using the HLMdiag package...               > library(HLMdiag)     > stdresid    <- HLMresid(random.model, level = 1, standardize = TRUE)     > diagnostics <- cbind(diagnostics, stdresid)      But I've hit a problem as the autoplot example calls ggplot2's `fortify()` function to calculate these and two additional measurements, sigma, the estimate of the residual SD when the corresponding observation is dropped from the model and hat, the diagonal of the hat matrix. Reading around I thought the `hatTrace()` function would be one part of the solution, but found posts on R-help indicating that it was removed from the lme4 package. Can anyone advise on how to calculate sigma and hat vectors? Or if anyone has canned solutions for diagnostic plots for lme4 that would be useful too.


title: Outline of benefits/costs of using Bayesian rather than Frequentist OLS and time-series?

text: I'm reading up on Bayesian techniques for Linear models and time series. While the texts are great at teaching the theory I would like to get a better handle on the pro's/con's of Bayesian analysis vs their frequentist equivalents. I can't find an article or even articles that do this. I really want it as a fairly high level rather than getting caught up in the technical detail right away. Ideally highlighting where and in what situations Bayesian is better and vice versa. Most just hint that the frequestist approaches also rely on prior information but are less honest about this, but it didn't go into much detail.


title: Sample size with respect to prediction in classification and regression

text: With respect to hypothesis testing, estimating samples sizes is done through power, and it is intuitive that increasing the same size increases the precision of estimated effects. But what about prediction for both classification and regression? What aspects of the prediction problem are influenced by sample size other than estimating the generalization error or RMSE for regression. In sum, properties that contribute to power in the hypothesis-testing setting differ from those that those that enable successful prediction through penalized regression/data mining/algorithmic modeling. How does sample size influence the success of these techniques? One paper that describes this idea is this one. Can anyone provide references for their comments? Thanks.


title: Testing Spearman's rho between-groups

text: I am analysing ordinal, non-normally distributed dependent measures using Spearman's rho. Specifically, we have an independent groups design, and I have found a significant relationship between the variables in one group but not the other (running separate Spearman's analyses for each group). I wish to test the difference in associations between the two groups, but am unable to find information concerning what procedure/test to use. I understand the Fisher r to z transformation may be used in the case of a Pearson's Product Moment Correlation, allowing differences in associations to be tested between-groups - is there a means of running a similar test of between-groups association with Spearman's rho?


title: Predicting ordered logit in R

text: Pardon my naïveté if this is a dumb question, but I'm new to R. I'm trying to do an ordered logit regression. I'm running the model like so (just a dumb little model estimating number of firms in a market from income and population measures). My question is about predictions.               nfirm.opr<-polr(y~pop0+inc0, Hess = TRUE)     pr_out<-predict(nfirm.opr)      When I run predict (which I'm trying to use to get the predicted y), the outputs are either 0, 3, or 27, which in no way reflects what should seem to be the prediction based upon my manual predictions from the coefficient estimates and intercepts. Does anyone know how get "accurate" predictions for my ordered logit model? **EDIT** To clarify my concern, my response data has observations across all the levels               >head(table(y))     y     0  1  2  3  4  5      29 21 19 27 15 16       where as my predict variable seems to be bunching up               > head(table(pr_out))     pr_out     0     1   2   3   4   5      117   0   0 114   0   0


title: The way an MA(q) model works

text: I am trying to understand the way MA(q) models work. For this purpose I have created a simple data set with only three values. I then adapted a MA(1) model to it. The results are shown below:               x<-c(2,5,3)     m<-arima(x,order=c(0,0,1))          Series: x      ARIMA(0,0,1) with non-zero mean           Coefficients:               ma1  intercept           -1.0000     3.5000     s.e.   0.8165     0.3163          sigma^2 estimated as 0.5:  log likelihood=-3.91     AIC=13.82   AICc=-10.18   BIC=11.12      While the MA(1) model looks like this: $$X_t = c +a_t - \theta*a_{t-1}$$ and $a_t$ is White Noise. What I cant figure out is how to get the fitted values:               library(forecast)     fitted(m)     Time Series:     Start = 1      End = 3      Frequency = 1      [1] 3.060660 4.387627 3.000000      I tried different ways, but I cant find out how the fitted values (`3.060660`, `4.387627` and `3.000000`) are calculated. I would be very thankful for an answer!


title: Clustering with constraint size

text: I want a clustering method which can let me apply constraints on the maximum and the minimum cluster sizes which I can implement. Otherwise suggestions for a combination of existing techniques which I could use to get the same result are welcome.


title: Beginner to prediction/statistics: Where do I start?

text: I sincerely apologize if there is another thread already that will answer this question. I'm so incredibly out of my league here that I don't even know what keywords to search for :-). I'm a computer programmer by trade, and while I have a basic background in math, statistics was never really my cup of tea. I currently work at a school and just finished developing a basic set of tools to help automatically collect and analyze data on our student's behaviors (this is a school for children with autism and other disabilities). So, we have a couple of year's worth of data for things like: given Billy, how frequently did he have Aggressions, Self-Injurious Behaviors, Drop, etc. Probably 6 - 10 "inputs" (I think that's the correct term) per student. We'll be adding more as well in the future. **What I'm curious about is this:** Are there any beginning tutorials out there that might show me some interesting things to do with this data (besides just graphing it?) For example, it would be interesting to be able to predict when Billy is likely to have a long string of aggressions given that these `x` other factors have been increasing lately. Or, there is an increasing trend of this behavior which is way out of whack with its previous values, that should raise a big red flag. I've been doing some basic Googling and this seems to be in the realm of "Statistical Data Mining"--some brief tutorials were found on Andrew Moore's site, but these just aren't detailed enough for me to really learn anything. I realize that this is akin to someone walking into Stack Overflow and saying "Hey, tell me how to write the next Facebook." So, if these are the sorts of things that I can only do with years and years of statistical experience, just let me know and I'll be on my way. However, I also know that while someone couldn't walk into SO and write the next Facebook in a few weeks, we could probably point them in the right direction to create a basic site for their dad's business, even if it would be a pretty basic site. Likewise, I'm not looking to create a genius AI capable of predicting student behavior down to the millisecond; rather, I'm just curious if there's any low-hanging fruit that a guy like me could pick up in a few weeks or months of diligent reading that might make for some interesting uses of this new data we've unlocked. **I'm open to online tutorials, books, textbooks, videos, open source programs and libraries, etc.**


title: How to exploit periodicity to reduce noise of a signal?

text: 100 periods have been collected from a 3 dimensional periodic signal. The wavelength slightly varies. The noise of the wavelength follows Gaussian distribution with zero mean. A good estimate of the wavelength is known, that is not an issue here. The noise of the amplitude may not be Gaussian and may be contaminated with outliers. **How can I compute a single period that approximates 'best' all of the collected 100 periods?** I have no idea how time-series models work. Are they prepared for varying wavelengths? Can they handle non-smooth true signals? If a time-series model is fitted, can I compute a 'best estimate' for a single period? How? A related question is this. Speed is not an issue in my case. Processing is done off-line, after all periods have been collected. **Origin of the problem:** I am measuring acceleration during human steps at 200 Hz. After that I am trying to double integrate the data to get the vertical displacement of the center of gravity. Of course the noise introduces a HUGE error when you integrate twice. I would like to exploit periodicity to reduce this noise. Here is a crude graph of the actual data (y: acceleration in g, x: time in second) of 6 steps corresponding to 3 periods (1 left and 1 right step is a period): ![Human steps](http://i.stack.imgur.com/44q8d.png) My interest is now purely theoretical, as http://jap.physiology.org/content/39/1/174.abstract gives a pretty good recipe what to do. It does not address periodicity. **Note:** I have asked this question on stackoverflow but it seems to be off- topic there.


title: Calculating p-value for a two-way ANOVA

text: I would like to compute the p-value of my two-way ANOVA. The score that I am using to detect significant samples is the eta score which is computed as SS(between)/SS(total). In many sites I saw that F=Var(between)/Var(within). I want to know if I can consider the eta score as the F-statistic value, and then compute the p-value? What if I want to compute the p-value for other scores such as eta-partial, or omega?   Is it meaningful if I calculate the p-value for them? Or is it only meaningful for the ratio Var(between)/Var(within)? What difference does it make if I use SS(between)/SS(within) instead of Var? The formula that I am using to calculate p-value is:               pvalue=-log10(betai(0.5 * df2, 0.5* df1, df2 / (df2 + df1 * eta)));        based on the book _Numerical Recipes in C_ , where df1 and df2 are degrees of freedom. Thanks for your help.


title: Finding similarity between a reference and few working documents

text: I have to find the similarity between a reference document and a set of documents in a repository . Here is my method :               1. I find the term document matrix for all the documents including the reference document.     2. The svd is calculated for this matrix.     3. I take the v array (the third result).     4. I transpose this matrix so that the each row represents a document.      5. The first row represents the reference document.      6. I find the cosine similarity between this row and the rest of the rows.       My doubts:   1. Since I have around 7 documents in my database, I get only a 8x8 `v` array (document matrix). So will I get a correct result if I find the cosine similarity with these 8 values alone?    2. Is such a method adopted generally?


title: What is a sparse estimator?

text: Can anyone point me to a reference, either book or paper, where I can find the precise definition of sparse estimator? Thanks


title: Central limit theorem with unknown variance

text: In my experiment I compute the average latency of operations per second. I would like to define N, i.e: how many times do I need to run my experiment to compute a close-to-real average latency? I figured that I could apply the CLT here, because If I repeat the same experiment 1000 times and plot a histogram, I get a normal distribution curve. Is the central theorem useful in my case? In the definitions I found, to be able to compute an estimated mean with a certain error, one needs to know the variance beforehand, and I don't know it.


title: Minimizing number of questions of questionnaire from past binary responses

text: We have data from 600,000 users that describes whether they observe each of 80+ binary features. That is, our data are a 600,000 x 80 binary matrix with user-profile. We know from inspection that some features are positively and negatively correlated. Some positive/negative features exclude others. Most users have less than 10 positive features in their profiles. We want to retrieve the profile of new users by asking them the minimum set of questions from the 80+ potential ones given this previous data. The idea is to give a small set of questions (5-10) to new users. Those should provide the maximum amount of information in order to "cut" the number of plausible subsequent questions. After a user has answered the first set of questions, we would like to ask a next set that, again, "cut" the number of subsequent questions faster. It seems reasonable to take into account positive and negative responses. Could you please provide me with some hints to how to implement this model? We would like to have:   * A way to measure the distribution of the expected number of questions given to each user.   * Some way to tune the "number of initial questions" provided.   * The model should be preprocessed in order to be able to react fast to user input.   * If possible, visualize the relationship between questions.   * If possible, be able to control the expected number of questions (I guess by discarding low correlations).   * If possible, update the model incrementally using new respondents (not critical) We plan to prototype in Python and then implement in .Net, but any other hint/code will be welcome.


title: comparing KNN with SVM on ROC chart

text: 1) I wish to compare the performance of KNN and SVM on ROC chart. a)For KNN, I obtained a linear line on ROC Chart instead of curved one.From x=0 until x=0.275, y= 0.65 and from x=0.35 until x= 0.95, y = 1.0. b)For SVM , the ROC curve behave normally . x=0, y=0 and at x=0.95, y=0.98. Which one is the better model of the two?


title: What is the method used in statistical classification when dealing with multi-dimensional discrete target labels?

text: If the training set is a set of $n$-tuples, with discrete labels, one can standardly use multinomial logistic regression (softmax), but what if the target labels are pairs of discrete values, or more generally, an $m$-tuple of discrete values? Is there any method I can use, or am I allowed to perform $m$ separate softmax regression classifications, for each target dimension?


title: How to calculate inter-rater reliability for just one sample?

text: I'm trying to compute a reliability rating for potentially malicious emails. Essentially a system where multiple people will view a suspect email, and assign it a "malice score". The idea is for emails where the raters reliably rate it malicious, to flag it and do stuff to it (discard, mark for review etc.) I've (tried to) read up about the various options...Kippendorf's Alpha, ICC etc, but these all need the same raters to rate multiple samples (which will not occur in my scenario - there will just be lots of samples being rated by lots of randomly selected users). Can anyone suggest something to try? (And/or let me know if I've completely misunderstood!) On consideration, I think I need to elaborate more: The goal is to quantify the degree of consensus among the random sample of raters for **each email**. With that information, we can automate an action for each email: e.g. If there is consensus the the email is bad/good, discard/allow. If there is significant disagreement, quarantine. Any interest in whether or not the whole system of randomly selected raters + randomly selected emails is overall reliable is more of a side-topic. Does that change anything?


title: Confidence interval and statistical significance in comparison of AUC

text: Recently, I have compared two correlated AUC with the method of Delong. Someone said that since the CI’s overlap, we cannot state the two models were different. I know that the method of Delong calculated correlated AUC, but I don't know how to answer his question in a formal way. Could someone give some help?


title: How to compare both precision and recall of two classifiers over one dataset?

text: Then two classifiers are evaluated on the same dataset, McNemar's test could be used to compare their overall classification accuracy, telling if the difference between them are significant. But what if we like to compare them separately on type I and type II errors, quantifying the difference on both precision and recall? It seems like McNemar's test could not be used at least for precision because of the randomness of the denominator. What are the suitable statistical tests for this case?


title: Transformations(CDF technique)

text: Consider the following short example of transformations. Let the joint density of X and Y be given by the unit square, i.e. $$f_{X,Y}(x,y) = \begin{cases} 1\ \quad 0<x<1\ \text{ and }\ 0<y<1 \\ 0 \quad \text{elsewhere} \end{cases}$$ Then the Cumulative Distribution Function of $Z=X+Y$ is given by: $$ F_Z = \begin{cases}\begin{array}{ll} 0\ & \text{ for }\ z<0 \\ \int_0^{z} \int_0^{z-x} dy\,dx\ & \text{ for }\ 0\leq z <1 \\1-\int_{z-1}^1 \int_{z-x}^1 dy\,dx\ & \text{ for }\ 1\leq z<2 \\1\ & \text{ for }\ 2\leq{z} \end{array}\end{cases} $$ I understand why we have to partition our CDF, what I am having trouble figuring out is why for the interval $[1,2)$ that specific form. What is the intuition here? Thanks.


title: Establishing the lack of randomness in lottery number selection

text: Related to this question, if I have 1500 or so jackpot results from a 6/49 lottery (numbers drawn, number of winners and prize per jackpot winner), how can I demonstrate that some numbers are less likely to be chosen by players than others? I don't have direct access to the distribution of numbers actually picked, of course. Take, for example, the hypothesis that players are more likely than would be expected by chance to pick numbers corresponding to dates, ie from (1-31): I find a positive correlation between the number of numbers > 31 and the prize-per-winner, and a negative correlation between the number of numbers < 32 and the prize-per-winner. Also, the sum of the numbers is positively correlated with the prize-per- winner, and the number of numbers > 31 is positively correlated with a rollover event, when no-one wins. What is the best way of approaching this with the data available?


title: Comparing linear vs. non-linear regression model for "best" fit

text: How, if at all, can one compare the "fit" of a simple linear vs. non-linear regression model to observed data? I apologize if I didn't search long/hard enough for the answer, but I cannot find anything concrete. -Patrick


title: Testing association between two categorical variables, with repeated experiments

text: Say you are counting sharks and whales in the Atlantic and Indian occeans and want to check if eg sharks prefer the Atlantic and whales the Indian occean, as described in the example here. But instead of just total counts (eg sharks observed in the atlantic), you keep counts for individual days. For example:               +--------+-----------+-----------+-----------+---------+---------+---------+       |        | Atlantic1 | Atlantic2 | Atlantic3 | Indian1 | Indian2 | Indian3 |       +--------+-----------+-----------+-----------+---------+---------+---------+       | whales |         9 |        3  |       11  |       4 |       5 |       6 |       | sharks |         2 |        5  |        1  |      10 |       5 |       8 |       +--------+-----------+-----------+-----------+---------+---------+---------+        You could take the means to get a 2*2 contingency table, and proceed as in link to the fisher exact test above, but wouldn't you then remove valuable information about the variance? Is there a way to test for association using all the information in the table above? Thanks


title: Coloring by principal component loading

text: In a diagram, I would like to indicate an item's principal component loadings by color. The idea is to map the item's loadings on three principal components (e.g., the first three) onto RGB values. The items in question can be represented by 20-dimensional unit vectors (in the specific case motivating this question, the items are protein sequences and I am analyzing their amino-acid composition). The principal components are derived from eigendecomposition of the correlation matrix of these vectors for a large number of items. One can envision simple mappings like: eigenvectors = $n \times n$ matrix $\mathbf{T}$, target $n \times 1$ vector = $\mathbf{v}$, projection $\mathbf{p} = \mathbf{v}'\mathbf{T}$, first three loadings of projection $\mathbf{n} = \mathbf{p_{i \in 1\dots 3}}$, and a vector of RGB values between 0 and 1 obtained from $(\mathbf{n/||n||}+[1\quad1\quad1]')/2$ (which first normalizes $\mathbf{n}$, yielding entries guaranteed to lie in [-1,1], and then maps [-1,1] entries to [0,1] entries). This is pretty simple and seems to work, but may have non-obvious blindspots; I am particularly worried about whether the simple normalizations lead to poor coverage of color space and thus poor contrast. A major goal, as with all color-coding, is to provide visual contrast. Is there a standard solution to this problem? If not, suggestions to improve this one? Updated: the intended application is to allow automatic and informative coloring of items outside of a scatterplot or other quantitative format. An example of such a format is shown below; each row is a protein sequence, and each colored cylinder is a region of the protein sharing certain characteristics, with the same color used to indicate repeated instances of a type of region. The limitation is that the similarities or differences in colors, in this example, have no meaning and are chosen for graphical design reasons. ![Fig. 3 from Schreiber et al. 2006](http://i.stack.imgur.com/seBo1.jpg)


title: Probability limit calculation

text: My class notes list the following steps for calculating a plim under classic errors in variables: $$ {\rm plim}\ \beta_1 = \frac{{\rm cov}(\beta_0 + \beta_1 x_1 + \epsilon - \beta_1 e, x_1)}{{\rm var}(x_1)}=\frac{\beta_1 {\rm var}(x_1) - \beta_1 {\rm cov}(x_1, e)}{{\rm var}(x_1)}... $$ How was the simplification in the numerator done?


title: Running standardized/unstandardized regression equation on SPSS

text: I am working on a regression assignment using SPSS. I have raw, centered, and z scores for the data. I already regressed the DV and IV and see the unstandardized and standardized coefficients under the coefficients output box. My standardized regression equation is: Y’=1727.528X+-6290.967 and unstandardized regression equation is: Y’=.633x How do I write out standardized and unstandardized regression equation for each score? Do I need to check the standardized and unstandardized boxes under the Linear Regression> Save..._menu for each IV?


title: Statistical test for testing equal proportions

text: I have a questionnaire item (N=6572) for which: 5064 (77.1%) subjects **responded** to the question, and 1508 (22.9%) subjects **did not respond** to the question. What statistical test should i use to test if the **proportion of subjects who responded = proportion of subjects who did not respond?** I am thinking one sample z test for proportions. Please advise.


title: How to fit a negative binomial distribution in R while incorporating censoring

text: I need to fit $Y_{ij} \sim NegBin(m_{ij},k)$, i.e. a negative binomial distribution to count data. However, the data I have observed are censored - I know the value of $y_{ij}$, but it could be more than that value. The log- likelihood is \begin{equation} ll = \sum_{i=1}^n w_i (c_i \log(P(Y_{ij}=y_{ij}|X_{ij})) + (1- c_i) \log(1- \sum_{k=1}^32 P(Y_{ij} = k|X_{ij}))) \end{equation} where $X_{ij}$ represent the design matrix (with the covariates of interest), $w_i$ is the weight for each observation, $y_{ij}$ is the response variable and $P(Y_{ij}=y_{ij}|X_{ij})$ is the negative binomial distribution where the $m_{ij}=exp(X_{ij} \beta)$ and $\alpha$ is the over-dispersion parameter. Does anyone know of an R package to tackle this problem?


title: How should I normalize my accelerometer sensor data?

text: I'm working with a large set of accelerometer data collected with multiple sensors worn by many subjects. Unfortunately, nobody here seems to know the technical specifications of the devices and I don't think they have ever been recalibrated. I don't have a whole lot of info about the devices. I'm working on my master's thesis, the accelerometers were borrowed from another university and altogether the situation was a bit intransparent. So, preprocessing on-board the device? No clue. What I do know is that they are triaxial accelerometers with a 20Hz sampling rate; digital and presumably MEMS. I'm interested in nonverbal behavior and gesturing, which according to my sources should mostly produce activity in the 0.3-3.5Hz range. **Normalizing** the data seems quite necessary, but I'm not sure what to use. A very large portion of the data is close to the rest values (raw values of ~1000, from gravity), but there are some extremes like up to 8000 in some logs, or even 29000 in others. **See the image below**. I think this makes it a bad idea to divide by the max or stdev to normalize. **What's the usual approach in a case like this?** Divide by the median? A percentile value? Something else? As a side issue, I'm also not sure if I should clip the extreme values.. Thanks for any advice! **Edit** : Here is a plot of about 16 minutes of data (20000 samples), to give you an idea of how the data is typically distributed. ![enter image description here](http://i.stack.imgur.com/nXQBU.png)


title: Prove this theorem related with specification tests

text: I'm not a mathematician and I need to prove what this theorem says. I think is easy and I know how it works, but definitively I'm not too much rigorous to make a demonstration. Can anybody help me? I enounce it: If the test statistic has a continuous distribution, then under H0 : θ = θ0, the p-value has a Uniform[0, 1] distribution. Therefore, if we reject H0 when the p-value is less than α, the probability of a type I error is α.   * Therefore, when H0 is true, the p-value is like a random draw from a Uniform[0, 1].   * On the other hand, if H0 is not true, the distribution of the p-value will tend to concentrate closer to 0.   * A large p-value can occur for two reasons:     * H0 is true, or     * H0 is false, but the test has low power   * Do not confuse the p-value with P(H0 | data). The p-value is not the probability that the null hypothesis is true. Thank you so much!


title: Computing a marginal posterior of a hierarchical model

text: This is my first time posting an actual homework question. Usually I have more local resources such as office hours and student peers but this time I am a little short on those. I also need to rest but I can't so I'm spinning my wheels a lot more than I ordinarily would. I have a hierarchical model of $m$ different binomials with $m$ different sample sizes, whose $m$ different $p$ parameters are governed by a population beta distribution. The beta distribution, of course, also has a prior. I constructed what I think is the correct joint posterior distribution for this model: $$P(p_1, p_2,\dots,p_{m}, \alpha,\beta|\mathbf X)\propto \pi(\alpha, \beta)\prod_{i=1}^{m} P(p_i|\alpha, \beta)P(x_i|p_i)$$ $$\propto\pi(\alpha, \beta) \prod_{i=1}^{m}p_i^{x_i+\alpha-1}(1-p_i)^{n_i-x_i+\beta-1} $$ The problem I am facing is that if I want to compute $$P(\alpha,\beta|\mathbf X)$$ It seems as though, since I have m parameters to remove I am looking at m integrals (integrals because the beta is continuous). This can't possibly be what I am expected to do unless one of the following is true:   1. I got the joint posterior incorrect   2. I should be doing this by sampling   3. There is an easier trick to mathematically computing a marginal posterior on this   4. I have the wrong idea of what marginalizing m parameters from a population distribution amounts to (eg. we just integrate once since they are all from the same distribution; I can't see what to make of the individual p parameters in that case)


title: Specifying a covariance structure: pros and cons

text: What are the benefits of specifying a covariance structure in a GLM (rather than treating all off-diagonal entries in the covariance matrix as zero)? Aside from reflecting what one knows of the data, does it   1. improve goodness of fit?   2. improve predictive accuracy on held-out data?   3. allow us to estimate the extent of covariance? What are the costs of imposing a covariance structure? Does it   1. add computational complications for estimation algorithms?   2. increase the number of estimated parameters, also increasing AIC, BIC, DIC? Is it possible to determine the right covariance structure empirically, or is this something that depends on your knowledge of the data-generative process? Any costs / benefits I didn't mention?


title: What's the best way to fit this very strange-looking data?

text: ![enter image description here](http://i.stack.imgur.com/4S96q.jpg) Hi all, Could anybody please help me? I am trying to learn if there is any good approach in modern statistics that can fit this data well, without over-fitting. Ok, here is another plot that I used "pch=".""... thanks for pointing that out to me! ![enter image description here](http://i.stack.imgur.com/PyPvK.jpg) Here is a 3D one: again, I want to be able to classify the {positive, negative}-ness of the z-axis... ![enter image description here](http://i.stack.imgur.com/g4JlY.jpg)


title: Cohen's d in one-way ANCOVA

text: I would like to calculate Cohen's d for pairwise comparisons in a one-way ANCOVA (4 groups). I have found this on the net: d-FAMILY MEASURE A common d-family measure (effect size index) is Cohen’s d. Like with the Oneway ANOVA, the effect size consists of the (numerator) mean difference between the contrasts (pairs) divided by the (denominator) averaged or pooled variability. We can estimate the average variability by taking the square root of MS’error from the analysis of covariance, which would standardize the mean difference in the metric of the adjusted scores. As such, we can estimate Cohen’s d with the following formula: d = difference of the 2 adjusted means/square root of MS'error. Is that correct?


title: References for use of symplectic geometry in statistics?

text: I have heard that many problems in Mathematical Statistics can be stated and solved in terms of Symplectic Geometry. Unfortunately this was a pretty vague statement and I am interested in something more concrete. Also there seem to be some books written form this perspective, but I couldn't find any.


title: SVM training method (or alg.)

text: I'm using SVM classification (Matlab) within my research works, and I want to know:   1. The advantages and disadvantages of each training algorithm, i.e., SMO, LS and QP     2. In general case, what is the suitable algorithm?     3. If the choice of one of them can impact the classification performance?     4. If the are any relationship between the choice of training method (or algorithm) and the choice of the kernel function.


title: How to write program that takes a vector and a CDF and performs the Kolmogorov-Smirnov test?

text: I've been reading the Wikipedia articles about various statistical tests, and trying to code up programs that perform them. However, the article on the Kolmogorov-Smirnov test is rather unclear to me. Can somebody explain, in simple language, how I would write a program that takes a list of real numbers and a CDF, and computes a p-value for how closely they match?


title: How do you find a peak rate in time series?

text: I have a series of timestamps, and I would like to know what's the maximum hourly rate of the event that these timestamps represent. You can generalize this to any units. I'm looking for a window (interval) of a certain size that contains most of these values. The obvious approach is to divide the set by this interval and then pick the one with most values. This is however imprecise, especially if the peak is near the limits of interval. Is there some better solution?


title: How to specify pulses/level-shifts in data when creating ARIMA in R?

text: I am new to ARIMA modelling. I understand most of the basic concepts, and I've read a lot of topics about ARIMA on this site. At present I am pretty comfortable with analysing ACF and PACF graphs, checking for seasonality and creating model from my data. Next step I want to deal with outliers and periodic pulses in my data. So, how do I do that using R? Should I just create some dummy variable and code periodic pulses with, for example 1, other outliers with 2 and everything else with 0, and then pass this variable to "xreg" option in arima command ? I am also afraid of mixing things up by not understanding terminology correctly. There is a lot of talk about pulses, level-shifts, seasonal pulses, local time trands, etc. Is there any guide or paper discussing all this stuff and how to deal with it ?


title: Unable to fit repeated measures in R

text: I've been trying to fit a repeated measures model in R using `lmer` but I get an error message. I've used SAS in the past and I was able to fit a somewhat similar model, but I can't seem to be able to do it in R. This is the SAS code:               proc glimmix data=dat;     class size day position subject;     model y = day size size*day position size*position/ ddfm=kr;     random day /subject=subject(size) type=ar(1) residual;     lsmeans size|day position ;      The following is the R code:               m1 = lmer(y ~ 1  + Size + Position + Size*Day + Position + (1 + Day|ID),           data = dat, REML = F)      and these are the error messages I get:               Warning messages:     1: In checkZrank(reTrms$Zt, n = n, control, nonSmall = 1e+06) :       number of observations <= rank(Z); variance-covariance matrix will be unidentifiable     2: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :       failure to converge in 10000 evaluations      The data consists of 242 observations taken from 22 individuals over a 11 day period, plus some other variables were observed as well (i.e. position). `Size`, `day`, `subject` and `position` are factors and `y` is the observed response. The data structure (I can provide the whole 242 observations if that would help):               > dput(head(dat,4))     structure(list(y = c(7.12820702168194, 7.70152948380659, 6.42263644198494,      6.9249926063685), Day = structure(1:4, .Label = c("1", "2", "3",      "4", "5", "6", "7", "8", "9", "10", "11"), class = "factor"),          Size = c(1, 1, 1, 1), Subject = structure(c(1L, 1L, 1L, 1L         ), .Label = c("101X", "104X", "111X", "112X", "118X", "119X",          "126X", "139X", "26X", "41X", "44X", "45X", "46X", "48X",          "55X", "63X", "67X", "68X", "79X", "84X", "87X", "92X"), class = "factor"),          Position = c(1, 1, 1, 2)), .Names = c("y", "Day", "Size",      "Subject", "Position"), row.names = c(NA, 4L), class = "data.frame")      I'd appreciate any help with this issue. Regards,


title: What is the posterior probability of the data given the model used for model averaging with Bayesian Logistic Regression?

text: I am trying to learn about Bayesian Model Averaging using Bayesian Logistic Regression (Genkin, A., Lewis, D. D., & Madigan, D. (2007). Large-scale Bayesian logistic regression for text categorization. Technometrics, 49(3), 291-304.), as found in the Weka data mining tool. I have two questions—   1. What are the ways in which models are generated using Bayesian Logistic Regression? (are subsets of predictive variables used at a time? are hyperparameters of the priors altered for each model?)   2. What is the posterior probability of the model (weights for model averaging) used? Is it the likelihood of the data and the parameters of the logistic regression? or, is it computed with information theoretic methods like AIC/BIC (I am unfamiliar with this)? Thank you.


title: Determing which sensor type is more informative

text: I have two sets of observations that correspond to the same set of quantities. The first set of observations were taken with one sensor type and the other set of observations with a different sensor type. Based on the observations and their associated quantities, how can I determine which sensor type is more informative?


title: What is the correct way of generating a p-value for a correlation with resampling?

text: I have a vector of gene expression values, across 20 patients. Each patient also has a glucose measure (a continuous numeric vector). I want to find how significant the real correlation value is compared to randomly distributed correlation values.   1. Real Pearson calculated by:               real<-cor(gene1, glucose)        2. I generate a distribution of 10K random correlations using sample with replacement              random_glucose <- replicate(10000,               cor(sample(gene1,replace=T),sample(glucose,replace=T)))      **How do I get a p-value of how significant the real is compared to the random values?** I have tried the following:   1. Using `pnorm` with the code:              pnorm_p <- 2*pnorm(q=abs(real),mean=mean(random_glucose),sd=sd(random_glucose),lower.tail=F)        2. I am not convinced the random values are normally distributed, so am also using `ecdf` with the code:              ecdf_calc<-function(real, random_glucose){       funny<-ecdf(random_glucose)       value<-funny(real)       if(value<=0.5){         res<-2*value         return(res)       }else{         res<-((1 - value)*2)         return(res)         }     }      What are the right ways to do it with `ecdf` and `pnorm`? I see different ways of doing this via Google and I don't understand the difference. E.g., methods using Z-scores as here; http://www.cyclismo.org/tutorial/R/pValues.html#calculating-a-single-p-value- from-a-normal-distribution And sometimes I see `1-pnorm(...)` and others I see `pnorm(..)` alone without the `1-`.


title: What is MA(q) model input in real world?

text: I understand the AR(p) model: its input is the time series being modelled. I'm completely stuck when reading about the MA(q) model: its input is _innovation_ or _random shock_ as it's often formulated. The problem is I can't imagine how to get an _innovation_ component having no model of the (perfect) time series already (i.e. I think $\varepsilon=X_{\rm observed}-X_{\rm perfect}$, and that's probably wrong). Moreover, if we can get this innovation component in-sample, how can we get it when doing a long- term forecast (model error term as a separate additive time series component)?


title: what is the return value of predict in the fGarch package

text: I have a question about a quit sophisticated model for a time series. Suppose $ \{X_t:0\le t\le T\}$ is a time series. The plot of autocorrelation function and partialcorrelation function suggest and ARMA model. However, I also want to model the volatility, hence I use a ARMA(p,q)-GARCH(1,1) model, say. This means $$ X_t=\mu_t+\sigma_t Z_t$$ $$\mu_t=\mu +\sum_{i=1}^p\phi_i(X_{t-i}-\mu)+\sum_{j=1}^q\theta_j(X_{t_j}-\mu_{t-j}) $$ $$ \sigma_t^2=\omega+\alpha_1(X_{t-1}-\mu_{t-1})^2+\beta_1\sigma_{t-1}^2$$ Then $\mu_t$ models the conditional expectation and $\sigma_tZ_t$ the conditional variance, where $\{Z_t\}$ is a strict white noise. In `R` we can use               model <- garchFit(formula=~arma(p,q)+garch(1,1),cond.dist="std",trace=F)      command where we have to specify $p,q$ and we decided to assume the white noise has a student $t$ distribution. Using the built in functin `predict` we can get a forecast for the next day:               predict(model,n=1)      However what is the return value of `predict`? Does it forecast $\sigma_{t+1}^2$, $\mu_{t+1}$ or $X_{t+1}$ and how can I get the other two forecasts within `R`?


title: Markov Transition Probabilities

text: Appreciate this is a very easy question - but I've only seen markov transition matrices with specified probabilities so am having a bit of a block with this. In a simple 2-state regime, if state 1 is normally distributed with mean 2 and variance 1 and state 2 norm distributed with mean 4 and variance 36, and an initial probability of being in state 1 of 0.75 - how do I work out the transition matrix? Thanks


title: How is the impulse-response function of a given system related to the autocorrelation function?

text: If I have the autocorrelation function of an observed system output, how does this relate to the impulse-response function of that system, if I don't have information about the input?


title: fitting multiple curvers at same time to set of points

text: Can we fit multiple second order curves to a set of points in $R^2$? Let's say we have set of points in the above mentioned space. I know that there could be at most 10 curves that can be fit in that space. I want to have the best possible fit to those points using arbitrary number of curves.


title: Z test could be a good choice ? and z-test for multiple samples?

text: I've two samples with mean value $\mu $ and standard deviation $\sigma$ (normal distribution). Can I use a Z-test to compare them ? Also if I don't know the sample size ? If I've five samples to compare ..what do I do ? A z-test for each comparison :\ or an ANOVA test (as for t-test) ?


title: Random forest - proof of convergence

text: I'm having some trouble understanding Leo Breiman's proof that the generalization error of a random forest converges as the number of trees increases (here's a link to the paper). At Appendix I he states that: > **Proof of theorem 1.2** : It suffices to show that there is a set of > probability zero C on the sequence space $ \Theta_1, \Theta_2, ... $ such > that outside of C, for all x, $$ \frac{1}{N}\sum_{n=1}^{N}I(h(\Theta_n,x)=j) > \rightarrow P_\Theta(h(\Theta,x)=j) $$ I'm failing to understand the following:   1. The sequence space of $ \Theta $'s is a sequence space as defined in wikipedia? are the elements of the sequence space infinite sequences of real or complex numbers?   2. What does a "set of probability zero C on the sequence space ..." mean? is there another way to describe C? Any help or explanation, partial or otherwise, is much welcomed. Thanks !


title: How do I estimate a prediction interval from a weighted sample?

text: I have a vector of data with their standard errors:               Estimate <- c(0.254719513441046, 0.130492717014416, 0.0386710035855823, 0.14118562325405, 0.160649388742147, 0.60363287936294, 0.173485345603584, 0.425817607348994, 0.128802795868366, 0.104136474748465)     SE <- c(0.126815201703205, 0.240179692822184, 0.248612907189712, 0.379800224602374, 0.0799874163236805, 0.170568135051654, 0.108163615496468, 0.0585237357996271, 0.16702614577514, 0.124308993809982)      I need to form a prediction interval for new elements that will be drawn from the same parent population. At first I thought I would do a linear regression to estimate the mean and the error on the mean:               summary(lm(Estimate ~ 1, weights = 1/SE^2))     Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  0.28014    0.04996   5.607 0.000331 ***      But when I try naively to call `predict()` here is what I get:               predict(lm(Estimate ~ 1, weights = 1/SE^2), newdata = data.frame(x = 0), interval = "prediction")             fit       lwr      upr     1 0.2801405 -2.860802 3.421083     Warning message:     In predict.lm(lm(Estimate ~ 1, weights = 1/SE^2), newdata = data.frame(x = 0),  :       Assuming constant prediction variance even though model fit is weighted      Did I do this right? Is the 95% confidence interval really [-2.86 3.42]? Do I need to worry about the warning?


title: I have an easy Text mining frequency question

text: I have a very simple question. I am doing a simple frequency analysis where I am looking at two documents and comparing the frequency of the top used words in each document. I want to do a comparison of the two documents in R, but the problem is that one document has 10x more words than the other document. What would be the best way to normalize so that I have two side by side frequency bar graphs that allow me to compare the most frequent terms. I hope that makes sense. No need to provide code. I can figure out the code I just want to know how to present the bar graphs. I was thinking (frequency of word x)/(total number of words in document y).   Thanks!


title: Time dummies in panel data -- absorbing effects?

text: I am conducting a data analysis. I have a panel with individual firms with firm-specific and macroeconomic variables. I would like to run an OLS regression adjusted for firm clustering effects and accounting for time effect (I have 449 firms and 36 time periods, unbalanced panel) The question arises: When I include time-dummies will the effect of the "overall" macroeconomic variable absorbed by the dummies, in other word: Will including time-dummies lead to multicollinearity between the time-dummies and the "macroeconomic" variables? I am wondering because different studies use this approach, although it appears more logical to me, to exclude the macroeconomic variables when including time-dummies. Did I fail to consider anything important here?


title: Which statistics to use in order to understand a dataset?

text: So I have a dataset that I will use to train a bunch of classifiers. I need to do that for my thesis. However I'm not sure which statistics are good to use to better understand the dataset and the features. I will mention here the statistics that I compute over the dataset, but if you know better ones or if I forgot something, or if there is something cool to use, so **please** tell me! It would be really nice to have in the thesis ;). I need to note that all features' values are between 0-1.0. So here is what I'm doing: 1- Covariance plot for the features. 2- Scatterplot matrix for the features. 3- Pearson correlation between the features. I have an additional question please. I though of calculating the information gain over the feature. But does it make sense to compute and report the information gain for _all_ the features only at the root level? Or is it better to report only the information gain on the split nodes? I think the latter makes more sense, right?


title: Fit measures for GMM Arellano-Bond estimator in R

text: A colleague and I have been working with difference GMM, i.e. the Arellano- Bond estimator, in R. Our option has been to use the pgmm command from the plm package. However, now I am struggling to test the fit of my models, since the package plm itself does not bring such functionality. Would there be a good way to test for the model fit of a difference-GMM model? Of course, an approach would be to test the square of the correlation between Y and Yhat, as in Bloom, Bond et al. (2001). However, the pgmm command does not give the predictions, only the fitted-values. Could anyone give any tips on that?


title: Adding coefficients to obtain interaction effects - can I add standard errors?

text: I posted this question earlier and am rewriting it in hopes of getting some guidance. I am using a weighted regression (after propensity score matching) to obtain estimates of the effects of a treatment (treat) on an outcome for different income quintiles. I include interactions of the quintiles and the treatment in the equation, as shown below:               regress outcome treat quintile2 quin3 quin4 quin5 treatXquin2 treatXquin3 treatXquin4   treatXquin5      To obtain the effect of treatment on quintile2 I add the coefficient on treat with the coefficient on treatXquin2. However, what do I do with the standard errors? Can I simply add the standard errors on the two coefficients together? Can I add the t-stats together? Any advice would be much appreciated. I'm not a stats expert but have a basic understanding of econometrics.


title: Monte Carlo experiment to estimate coverage probability

text: I'm working on a problem as follows for a course that I'm auditing: > Suppose a 95% symmetric t-interval is applied to estimate a mean, but the > sample data are non-normal. Then the probability that the confidence > interval covers the mean is not necessarily equal to 0.95. Use a Monte Carlo > experiment to estimate the coverage probability of the t-interval for random > samples of $\chi^2(2)$ data with sample size $n = 20$. Here is the current state of my R code:               alpha = 0.05;     n = 20;     m = 1000;          UCL = numeric(m);     LCL = numeric(m);          for(i in 1:m)     {         x = rchisq(n, 2);         LCL[i] = mean(x) - qt(alpha / 2, lower.tail = FALSE) * sd(x);         UCL[i] = mean(x) + qt(alpha / 2, lower.tail = FALSE) * sd(x);     }          # This line below is wrong...     mean(LCL > 0 & UCL < 0);      The problem is that the result is $0$. Am I approaching this question incorrectly? What exactly does coverage probability mean...? Thank you :)


title: How can i make a fraud detection dataset (I have the data ready but unordered)

text: I'm a little confused with the creation of the dataset for a fraud detection predictive model. Here i put a link with a sample of the dataset that I made. (the real dataset have ~950.000 clients). https://www.dropbox.com/s/r2r0f8s7nqvi59d/sample.csv I'm developing a Logistic Regression and the dependent variable is Fraud. So i have a problem with the dataset because I don't know how to arrange this, because I have 67 months with consumption for each client. And I'm thinking that the right way to arrnge the dataset is something like this. https://www.dropbox.com/s/ejj7v405xj0s9m7/Arrange.PNG In the real dataset in 1 row I put the ID (Unique Number) and the 67 months with consumption, so in the picture that I attached (That is the way that i think could be the best) there are 67 rows for each client, that means tha in row 1 i have January.08, in the row 2 i have February.08 and so on. So for each client I will have 67 rows and the total dataset will be 67(months)*950.000(clients) So I want an advice because I'm learning from this and I want to improve my work. I'm very thankful for any help Regards. Cristian Valenzuela


title: Moving average confusion

text: I have come to understand that a 3x5 MA (moving average) is equivalent to a 7 term MA. However, on applying the same to a time series such as below, I did not get the exact final time series values (one value, in this case) but approximately close. Made me wonder, if I did it correctly. The series I considered was: 0.067, 0.133, 0.2, 0.2, 0.2, 0.133, 0.067 Could someone highlight how exactly are the two calculated. My approach: 7-MA : add them up and divide by 7 = 0.1423 3x5 MA: a. Get the 5-MA as 0.16, 0.1732,0.16 b. Get the 3-MA of above as: 0.1644 (!= 0.1423) Am I doing it correctly?


title: Examples of studies using p < 0.001, p < 0.0001 or even lower p-values?

text: I come from the social sciences, where p < 0.05 is pretty much the norm, with p < 0.1 and p < 0.01 also showing up, but I was wondering: what fields of study, if any, use lower p-values as a common standard?


title: Operations on probability distributions of continuous random variables

text: How do probability distributions of continuous random variables transform under functions? I.e. I have a random variable, X, drawn from a normal distribution with mean 0 and variance 1. What is the probability distribution associated with sin(X)? ![Histograms mimicking probability density functions of X and sin\(X\)](http://i.stack.imgur.com/KY5Wc.png) More generally, what are the rules for transforming continuous random variables? If we know the PDF and CDF of two random variables, X,Y what is the PDF/CDF of Z=X*Y ? How about Z=X^Y ? How about Z = sin(X+Y)+3? Are there Computer Algebra Systems which can compute this symbolically? Is this possible generally? If not, for what class of probability distributions and functions is it possible? Note: excuse the plot. These are obviously noisy histograms and obviously not to scale (areas under the curves do not match and so could not both sum to one). Hopefully the plot does get the point across though. Given the blue distribution describing X, I want the red distribution describing sin(X).


title: How to implement optimal scaling of categorical variables for fuzzy clustering?

text: I am referring this paper "Fuzzy e-Means Clustering of Mixed Databases Including Numerical and Nominal Variables" from IEEE for implementation and learning. Kindly correct my understanding so far:   1. The category score is updated using eq: 18.    2. Have to normalize the category score obtained as above. To normalize, subtract mean, and divide by standard deviation.   3. To normalize, mean is calculated from left hand side of eq: 19,variance is calculated from right hand side of eq: 20 [from variance calculate standard deviation] while calculating mean from L.H.S of eq: 19: $1_{n}^{T}U_{c}^{\theta}q_{j}$ does the variable $U_{c}^{\theta}$ asks us to use entire membership matrix or membership of any one cluster [one column of membership matrix]   4. Compute the normalized category score vector q_j. I am stuck here in implementation so that these constraints in eq:19,20 are satisfied. Any help is greatly appreciated. Thanks in advance. Yogalakshmi.J


title: Probability of Multiple Observations of Multiple Objects

text: **pr** = physical RFID readers that are deployed on the production line in the factory.   **p(pr)** = probability of the physical RFID reader.   **o** = product part moving on the production line in the factory.   **O** = {o1, o2, o3 ... on} is the set of product parts that would move on the production line in the factory.   **re** is a read event and is defined as re = (o, pr, t); where o is the object detected/read by the physical reader pr at time t.   **p(re)** = probability of the read event, In our scenario p(re) is equal to the probability of the physical RFID reader that triggered the read event i.e. p(re) = p(pr) **Product Part Movement and Sequences:** Product parts move in the factory in a sequence. Figure 1 displayed below shows this movement of product parts on the production line in the factory. ![enter image description here](http://i.stack.imgur.com/NQA7S.png) Physical RFID readers (pr) are electronic devices that are accurate with a certain probability. They may report the detection of a product part o when there may be no product part on the production line. Or they may not detect a product part o when there was actually a product part on the production line. The probability of RFID readers p(pr) to detect a product part o having an RFID tag depends on the number of tags it has to detect, the distance between the reader and the tag and other such things. **Partial Sequence:** A partial sequence (ps) is a sequence of product parts of the form (o1 < o2). The partial sequence (o1 < o2) denotes that product part o1 is ahead of product part o2 on the production line. **Sequence Detection Rules:** **Rule 1:** If pr1 detects a product part o1 at time t1 and then detects product part o2 at a later time t2, then o1 is ahead of o2 on the production line, i.e. (o1 < o2). In this scenario, the two read events would be re1 = (o1, pr1, t1) and re2 = (o2, pr1, t2). Note: In this scenario, the physical reader is same for both read events. ![enter image description here](http://i.stack.imgur.com/MDC5j.png) **Rule 2:** If pr1 is deployed before pr2 on the production line i.e. (pr1 < pr2), then a product part o1 detected by pr2 at time t1 is ahead of product part o2 detected by pr1 at the same time t1, i.e. (o1 < o2). In this scenario the two read events would be re1 = (o2, pr1, t1) and re2 = (o1, pr2, t1). Note: In this scenario, the time of detecting the product parts is same for both read events. ![enter image description here](http://i.stack.imgur.com/10TzN.png) **Rule 3:** If pr1 is deployed before pr2 on the production line (pr1 < pr2), then a product part o1 detected by pr2 at time t1 is ahead of product part o2 detected by pr1 at a later time t2, i.e. (o1 < o2). This is obvious from Rule 1 and 2. **p(ps)** = probability of the partial sequence. **Question 1:** What would be the probability of the partial sequence of two product parts (o1 < o2) ? **Extended Sequence:** Let (o1 < o2) be a partial sequence ps1.   Let (o2 < o3) be a partial sequence ps2. Since product parts move on the production line in sequence and one after the other, we can apply transitivity and extend the two partial sequences into an extended sequence (es) of 3 product parts (o1 < o2 < o3) ![enter image description here](http://i.stack.imgur.com/lfbEc.png) **p(es)** = probability of the extended sequence. **Question 2:** What would be the probability of the extended sequence of three product parts (o1 < o2 < o3) ? **Question 3:** What would be the probability of the extended sequence of n product parts (o1 < o2 < o3 .... < on) ? **Note:** This is not a homework assignment. I am working on a project that involves real-time product monitoring using RFID readers. The problem and the concepts are my personal work.