title: Adjusting a subset of P values to correct for multiple testing

text: I have a list of 300 $p$ values selected from all $p$ values obtained in multiple testing such that $p<.05$. I do not know the values of $p>.05$ or their number. The $p$ values are not corrected for multiple testing, and I would like to do so for the 300 raw $p$ values that I have. Is it OK to apply the Benjamini-Hochberg correction method, for instance with `p.adjust()` in R, on my list of $p$ values even though they have been pre- selected for being less than $.05$? If not, is there another method that can be used?


title: Understanding proof of a lemma used in Hoeffding inequality

text: I am studying Larry Wasserman's lecture notes on Statistics which uses Casella and Berger as its primary text. I am working through his lecture notes set 2 and got stuck in the derivation of lemma used in Hoeffding's inequality (pp.2-3). I am reproducing the proof in the notes below and after the proof I will point out where I am stuck. * * * **Lemma** Suppose that $\mathbb{E}(X) = 0$ and that $ a \le X \le b$. Then $\mathbb{E}(e^{tX}) \le e^{t^2 (b-a)^2/8}$. **Proof** Since $a \le X \le b$, we can write $X$ as a convex combination of $a$ and $b$, namely $X = \alpha b + (1 - \alpha) a$ where $\alpha = \frac{X-a}{b-a}$. By convexity of the function $y \to e^{ty}$ we have $e^{tX} \le \alpha e^{tb} + (1 - \alpha) e^{ta} = \frac{X-a}{b-a} e^{tb} + \frac{b-X}{b-a} e^{ta}$ Take expectations of both sides and use the fact $\mathbb{E}(X) = 0$ to get $\mathbb{E}(e^{tX}) \le \frac{-a}{b-a} e^{tb} + \frac{b}{b-a} e^{ta} = e^{g(u)}$ where $u = t(b-a)$, $g(u) = -\gamma u + \log(1-\gamma + \gamma e^{u})$ and $\gamma = -a /(b-a)$. Note that $g(0) = g^{'}(0) = 0$. Also $g^{''}(u) \le 1/4$ for all $u > 0 $. By Taylor's theorem, there is a $\varepsilon \in (0, u)$ such that $g(u) = g(0) + u g^{'}(0) + \frac{u^2}{2} g^{''}(\varepsilon) = \frac{u^2}{2} g^{''}(\varepsilon) \le \frac{u^2}{8} = \frac{t^2(b-a)^2}{8}$ Hence $\mathbb{E}(e^{tX}) \le e^{g(u)} \le e^{\frac{t^2(b-a)^2}{8}}$. * * * I could follow the proof till $\mathbb{E}(e^{tX}) \le \frac{-a}{b-a} e^{tb} + \frac{b}{b-a} e^{ta} = e^{g(u)}$ but I am unable to figure out how to derive $u, g(u), \gamma$.


title: Recognizing relationship/pattern

text: I have some data that I know are related but I can't seem to find the exact relationship. Was hoping someone could help find it. Left column data is related to the right column.               1160786507      E02K     1432646245      Udre     1211119187      H06S     1214345830      HARF     1211117653      H00U     1429221446      U00F     1430466865      UC11


title: TF-IDF cutoff percentage for tweets

text: I'm currently trying to analyze Tweets and classify them as either positive, negative, or neutral using the NLTK library in Python. I can see that there's potential in the approach that I'm taking, however, I'm having trouble with my feature selection process. Indeed, input from Twitter isn't exactly conventional, so there aren't many examples of using TF-IDF to remove uninteresting words in Tweets. My implementation of the algorithm is working _alright_ , but I don't now how many words I should be cutting out of my feature list based on the TF-IDF scores: 25%? 15%? 30%? Additionally, and this is the main problem, since Tweets are different than normal English, many uninteresting words, such as 'the' and 'is', don't occur that often, and are therefore assigned good scores by the TF-IDF algorithm. Now, I know I can use a stoplist, but if I could avoid it, I'd be happy (I have a multitude of reasons _not_ to use a stoplist). So, to summarize, is there any way to circumvent this issue of uninteresting words getting good scores with TF-IDF? Edit: I am using a Naive Bayes classifier.


title: Technical variation versus real signal

text: I'm testing for difference in a continuous outcome under three different conditions. Under condition A I take a measurement of the outcome. I do this twice for the same sample. Example values could be 2.2, 2.1. These are "technical" replicates that come from the same biological source I do the same for four "biological" replicates for condition A:               A1, measure 1: 2.2     A1, measure 2: 2.1     A2, measure 1: 2.0     A2, measure 2: 2.1     A3, measure 1: 1.9     A3, measure 2: 1.8     A4, measure 1: 1.5     A4, measure 2: 1.6      I also have conditions B, C, and D, with two "technical" replicates in each of four "biological" replicates. How would I test test for mean differences (ANOVA) that best accounts for both the technical and biological variation? I wouldn't want to fit a model counting each measurement as a separate observation, because each pair comes from the same biological sample. I'm assuming there must be a better way than just averaging over the pairs. Bonus: how do you do this in R? Assuming I have data that looks like this:               > data        condition sample measurement outcome     1          A      1           1     2.2     2          A      1           2     2.1     3          A      2           1     2.0     4          A      2           2     2.1     5          A      3           1     1.9     6          A      3           2     1.8     7          A      4           1     1.5     8          A      4           2     1.6     9          B      1           1     1.7     10         B      1           2     1.6     11         B      2           1     1.5     12         B      2           2     1.6     13         B      3           1     1.4     14         B      3           2     1.3     15         B      4           1     1.0     16         B      4           2     1.1     17         C      1           1     2.4     18         C      1           2     2.3     19         C      2           1     2.2     20         C      2           2     2.3     21         C      3           1     2.1     22         C      3           2     2.0     23         C      4           1     1.7     24         C      4           2     1.8      I probably wouldn't want to do something like this:               summary(lm(outcome~condition, data=data))      Thanks in advance.


title: Convert probit regression coefficient into correlation

text: Most coefficients (e.g., odd ratios, mean differences, etc.) can be converted into a (Pearson-like) correlation coefficient. Q: Is there a conversion from probit coefficient --> correlation coefficient?


title: How to fit a Weibull distribution to input data containing zeroes?

text: I'm trying to reproduce an existing prediction algorithm, handed down by a retired researcher. The first step is to fit some observed data to a Weibull distribution, to obtain a shape and scale which will be used for predicting future values. I'm using R to do this. Here's an example of my code:               x<-c(23,19,37,38,40,36,172,48,113,90,54,104,90,54,157,51,77,78,144,34,29,45,16,15,37,218,170,44,121)     f<-fitdistr(x, 'weibull')      This works fine unless there are any zeroes in the input array, which causes it to fail completely. The same thing happens in SAS. As I understand it, this is because one of the steps in calculating the Weibull distribution is taking the natural log, which is undefined for 0. Is there a reasonable way to work around this? The best I've found so far is to add 1 to all of my input values, fit the curve, and then subtract one from my predicted values ("shift" the curve up and then back down by 1). This fits the previously predicted data fairly well, but seems like it must be a wrong way of doing so. edit: The values in the input array are observed, real-world data (the number of occurrences of something) for a range of years. So in some years the number of occurrences was zero. Whether it's the best way or not (I agree that it may not be), the original algorithm author claims to have used the Weibull distribution, and I have to try to replicate their process.


title: Roy model question

text: I am referring to G.S. Maddala: Limited Dependent and Qualitative Variables in Econometrics, pages 257-258. I add the relevant screenshots here: ![p1](http://i.stack.imgur.com/fuVav.png) ![p2](http://i.stack.imgur.com/M8F6e.png) My question is, why is $\sigma_{2u}>\sigma_{1u}$? I am not seeing this? I tried the following: \begin{align} \sigma_{2u}>\sigma_{1u}\\ (\sigma^2_{2}-\sigma_{12})/\sigma>(\sigma_{12}-\sigma^2_{1})/\sigma\\ \sigma^2_{2}-\sigma_{12}>\sigma_{12}-\sigma^2_{1}\\ \sigma^2_{2}+\sigma^2_{1}>2\sigma_{12} \end{align} But why should $\sigma^2_{2}+\sigma^2_{1}>2\sigma_{12}$ be fulfilled? This is not necessarily the case?


title: Why do a density plot and a rug plot seem to disagree?

text: The second peak of the `density plot` is large in this example. Why does the rug representation of the data--which seems to show few high values--not appear to match the much higher density estimated there? How can one make rug plots less misleading? ![enter image description here](http://i.stack.imgur.com/Zj9ie.png) Here's its `R` code:               plot(density(Boston$tax))         rug(Boston$tax, col=2, lwd=3.5)


title: Hand hygiene observations

text: How much validity is there in doing external hand hygiene observations in a hospital where another facility observes for a total of 6 hrs. twice a year and reports that data out, when the observed facility already has their own internal hand hygiene observers to obtain there own data monthly?


title: Fourier phase randomization

text: Can anybody please explain or point to an online resource that explains the Fourier phase randomization technique? I encountered it in the context of comparing the cross-correlation of two signals before and after Fourier phase randomization. What does it mean when the cross-correlation disappears after Fourier phase randomization? Thank you.


title: Internal vs external cross-validation and model selection

text: My understanding is that with cross validation and model selection we try to address two things: **P1**. Estimate the expected loss on the population when training with our sample **P2**. Measure and report our uncertainty of this estimation (variance, confidence intervals, bias, etc.) Standard practice seems to be to do repeated cross validation, since this reduces the variance of our estimator. However, when it comes to reporting and analysis, my understanding is that **internal validation is better than external validation** because: It is better to report:   * The statistics of our estimator, e.g. its confidence interval, variance, mean, etc. on the full sample (in this case the CV sample).  than reporting:   * The loss of our estimator on a hold-out subset of the original sample, since: (i) This would be a **single measurement** ( ** _even if we pick our estimator with CV_** ) (ii) Our estimator for this single measurement would have been trained on a set (e.g. the CV set) that is smaller than our initial sample since we have to make room for the hold-out set. This results in a **more biased** (pessimistic) estimation in **P1** . Is this correct? If not why? ### Background: It is easy to find textbooks that recommend dividing your sample into two sets:   * The **CV** set, which is subsequently and repeatedly divided into **train** and **validation** sets.   * The **hold-out** (test) set, only used at the end to report the estimator performance My question is an attempt to understand the merits and advantages of this textbook approach, considering that our goal is to really address the problems **P1** and **P2** at the beginning of this post. It looks to me that reporting on the hold-out test set is **bad practice** since the analysis of the CV sample is more informative. ### Nested K-fold vs repeated K-fold: One can in principle combine **hold-out** with regular **K-fold** to obtain **nested K-fold**. This would allow us to measure the variability of our estimator, but it looks to me that for the same number of total models trained (total # of folds) repeated K-fold would yield estimators that are **less biased** and **more accurate** than nested K-fold. To see this:   * Repeated K-fold uses a larger fraction of our total sample than nested K-fold for the same K (i.e. it leads to lower bias)   * 100 iterations would only give 10 measurements of our estimator in nested K-fold (K=10), but 100 measurements in K-fold (more measurements leads to lower variance in **P2** ) What's wrong with this reasoning?


title: sign of loglikelihood during parameter estimation

text: I need to optimize parameters of a Gaussian shaped function to best fit my data points using maximum likelihood. I first make an initial estimate for the parameters (mu, sigma and peak) and apply a gradient descent (fmincon in Matlab) algorithm to optimize these parameters. The cost function is the sum of the negative log values of probabilities of observing each data point given the model parameters and the model noise parameter.               sum(-log( normpdf( Y - model( X, params) ,0,error_sigma) ));      where, Y: observations normpdf : normal probability densitiy function. model: a function that takes the x values and parameters to be optimized params : parameters of the model error_sigma : noise of the data I optimize the model parameters and the noise parameter simultaneously to obtain highest -loglikelihood values. I am puzzled by obtaining negative likelihood values. I would have estimated that -log of probabilities should always evaluate in positive values. However, the matlab function normpdf, especially for small spread values, returns probability values which are bigger than 1. This is because the integral of the normpdf function does nearly never sums up to 1. Therefore I am thinking that I should be using the normal cumulative function for optimization, instead of normpdf in order to get valid loglikelihood values. However, the parameters never converges as good as the normpdf case, when I use normcdf. I hope I am clear.


title: Creating Correlations

text: I have three variables, A, B, and C. It seems obvious that if X = A-B and Y = C-B that there should be a correlation between X and Y. When I've done this in matlab with random numbers, this seems to be the case with a mean r of about 0.5. The fact that there is a correlation makes intuitive since because X and Y are sharing the variance in B; but when I've tried to figure out why it's 0.5, I've had less luck. Any hints as to why this is the case would be greatly appreciated.


title: How to calculate pooled variance of two groups given known group variances, means, and sample sizes?

text: Say there are $m+n$ elements split into two groups ($m$ and $n$). The variance of first group is $\sigma_m^2$ and variance of second group is $\sigma^2_n$. The elements themselves are assumed to be unknown but I know the means $\mu_m$, $\mu_n$ and $\mu_{(m+n)}$. The variance doesn't have to be unbiased so denominator is $(m+n)$ and not $(m+n-1)$. Is there a way to calculate the combined variance $\sigma^2_{(m+n)}$?


title: How to summarize measurement variability for repeated measures?

text: I have a data set where different measures are taken repeatedly (3) times for multiple samples, like so:                  A B C      X  ...      Y  ...      So, for people X, Y, etc I have measures A, B, C and so on. Each of A, B, C are measured 3 times for each sample. That is, for each person, I have 3 measures of A (A1, A2, A3), 3 measures of B (B1, B2, B3), for each of the columns. What I'd like to know is, how can I summarize the variability or the reliability of each column? So, over the entire data set, how reliable is the measurement for A? Thanks for the help!


title: What analyses do I need to run?

text: Long post. Question at the bottom. My data look likes this: ![enter image description here](http://i.stack.imgur.com/BSKqG.png) It consists of the test set accuracies of a series of 50 models formed on training sets. Using random sampling of the raw observations, I created 100 pairs of test & train sets in a 1:2 ratio. Each test/train pair contained all the raw observations, but each pair was partitioned differently. Running the 100*50 = 5000 iterations took about 40 hours of wall-time, so in practical terms, I cannot significantly increase the number of randomized test/train sets from 100 to, say, 1000 (an order of magnitude). I calculated the mean of the 100 accuracies of each model. They look like this: ![enter image description here](http://i.stack.imgur.com/mOzUK.png) The series of 50 models is parametrized by a single parameter. By eye, it appears to have three regions: 1:4, 5:10 and 11:50. Here are the linear models for those three regions, as created by the R function `lm(y~x)`: ![enter image description here](http://i.stack.imgur.com/no3QI.png) But that is just 'by eye'. I'd like to know what analyses would allow me to say something more rigorously. I'm particularly interested in the counter-intuitive region 1:4. For example, the mean for model #2 is 0.59 greater than model #1. Thus it seems that #2 makes better predictions than #1. For each of the 100 test/train pairs, I calculated the difference in accuracy between model #2 and model #1. They look like this: ![enter image description here](http://i.stack.imgur.com/a3uKd.png) With 64 of the values greater than 0, it again seems that #2 makes better predictions than #1. Their histogram looks like this: ![enter image description here](http://i.stack.imgur.com/wDhav.png) I've no training in stats. I ran a KS test (appropriately or not) to test if the values #2-#1 follow a Normal distribution with mean 0, `ks.test((M2-M1)/sd(M2-M1),'pnorm')`. This gave a D value of 0.184 and P value of 0.0023. I'm guessing the D and P values mean that the #2-#1 does not follow a Normal distribution, maybe because the mean is not zero or the shape is not Normal. I ran a similar second KS test `ks.test(((M2-M1)-mean(M2-M1))/sd(M2-M1),'pnorm')`, which gave D=0.054, P-value=0.94, which I'm guessing imply that the shape of #2-#1 is close to Normal. What analyses will allow me to make rigorous statements about 1) several regions of model progression, and 2) model #2 being a better predictor than model #1?


title: Power analysis of repeated measures ANOVA using simulation in R?

text: I would appreciate another opinion on this simulation approach to estimating power for a within-subjects (or repeated-measures) ANOVA; in this case a 2x2 factorial design. In other words, is there anything WRONG with this method that I've missed? The code below is R. * * *               library(ez)          nsub = 30     nconds = 4          nsims = 1000          #create an empty matrix that will be populated with p-values      p = matrix(NA, nrow=nsims, ncol=3)          #subject vector     sub = sort(rep(1:nsub, nconds))          #2x2 factorial design     cond = data.frame(x1=c('a','a','b','b'), x2=c('c','d','c','d'))          # fixed effects      intercept = rep(-1, nconds)     true.effect.x1 = c(0, 0, .5, .5)     true.effect.x2 = c(0, .5, .5, .5)     X = rep((intercept + true.effect.x1 + true.effect.x2),nsub)          #simulation loop     for (x in 1:nsims){            #random effects            #relatively large subject-specific variance in intercepts       sub.intercept = rep(rnorm(nsub, mean=0, sd=2), times=1, each=nconds)             #relatively small by-subject adjustments to the effects       sub.effect = rep(rnorm(nsub, mean=0, sd=0.05), times=1, each=nconds)             #unexplained error       error = rnorm(nsub*nconds, mean=0, sd=1)            #simulated dependent variable       observed.y = X + (sub.intercept + sub.effect + error)            #place everything in a data frame       df = cbind(sub, cond, observed.y)       names(df) = c('sub','x1','x2','y')       df$sub = as.factor(df$sub)            #extract the p-values for each effect from a repeated measure ANOVA (ezANOVA from 'ez' package)       p[x,1] = ezANOVA(data=df, dv=.(y), wid=.(sub), within=.(x1, x2))$ANOVA[1,5]       p[x,2] = ezANOVA(data=df, dv=.(y), wid=.(sub), within=.(x1, x2))$ANOVA[2,5]       p[x,3] = ezANOVA(data=df, dv=.(y), wid=.(sub), within=.(x1, x2))$ANOVA[3,5]     }          ###### p-values < .05 ? ######     sig.x1 = ifelse(p[,1] <= .05, 1, 0)     sig.x2 = ifelse(p[,2] <= .05, 1, 0)     sig.int = ifelse(p[,3] <= .05, 1, 0)          ###### Histograms ######     par(mfrow=c(3,1))     hist(p[,1], 20, xaxp=c(0, 1, 20), col=2, main = paste('power X1:', mean(sig.x1 * 100), '%  with ', nsub, 'subjects'))     hist(p[,2], 20, xaxp=c(0, 1, 20), col=2, main = paste('power X2:', mean(sig.x2 * 100), '%  with ', nsub, 'subjects'))     hist(p[,3], 20, xaxp=c(0, 1, 20), col=2, main = paste('power interaction:', mean(sig.int * 100), '%  with ', nsub, 'subjects'))


title: Datasets constructed for a purpose similar to that of Anscombe's quartet

text: I've just come across Anscombe's quartet (four datasets that have almost indistinguishable descriptive statistics but look very different when plotted) and I am curious if there are other more or less well-known datasets that have been created to demonstrate the importance of certain aspects of statistical analyses.


title: Does K-means incorporate the K-nearest neighbour algorithm?

text: I was watching this tutorial on K-means clustering and from what I understand K-means is:   1. Randomly generate the centroids for k clusters   2. Create a classification model dividing into k regions (Do we use kNN here?).   3. Generate new centroids for each region   4. Repeat 2-3. In the video I was watching regarding K-means, it appears that the lecturer may have used kNN to create the dividing regions though I can't see anything online regarding if kNN is used in K-means.


title: Visualizing association results after adjusting for confounders

text: I'm trying to find a way to visualize the results of an association analysis where I corrected for confounding variables. I have a set of cytokine data (amount of protein in the blood) from a set of patients infected and uninfected with HCV. The difference between the two is minimal when we test/visualize the raw data. However, when I adjusted for things like Race, gender, age, and disease status (by building a simple linear model) the p-values improved drastically. However, I'm having trouble finding a way to show this visually. Normally I would just show a box-plot or bee-swam plot or something similar. But the raw data is unimpressive. Does anyone have any ideas on ways to show this beyond reporting a p-value and the confounders/effect sizes?


title: Numerical integration problem

text: I thought I would write this to Stack Overflow, but LaTeX Math is not working there. I am trying to program the following equation in R: $\hat{\sigma}^2=1/n\sum_{i=1}^{n}[\int \hat{K}(u)\sum_{i=1}^{n}W_{1i}(u)\{dN_{1i}(u)-\hat{\lambda}(u)Y_{1i}(u)du\}]^2+1/n\sum_{i=1}^{n}[\int (1-\hat{K}(u))\sum_{i=1}^{n}W_{2i}(u)\{dN_{2i}(u)-\hat{\lambda}(u)Y_{2i}(u)du\}]^2$ Where, $\hat{\lambda}(u)=\frac{dN_{1i}+dN_{2i}}{Y_{1i}+Y_{2i}}$ $\hat{K}(u)=\frac{W_{2i}(u)Y_{2i}(u)}{W_{1i}(u)Y_{1i}(u)+W_{2i}(u)Y_{2i}(u)}$ The whole thing is based on a lifetime data set and $dN_{ji}(u)$ and $Y_{ji}(u)$ for $j=1,2$ are indicators of death and being in the risk set for subject $i$ respectively like in usual survival counting processess. $W_{ji}(u)$ is just a weight function. I have discrete lifetimes and in the formula of $\hat{\sigma}^2$ I am planning to conduct a numerical integration. But at each death time point (say, $u_k$) it seems to me that $\int \hat{K}(u)\sum_{i=1}^{n}W_{1i}(u)\{dN_{1i}(u)-\hat{\lambda}(u)Y_{1i}(u)du\}$ $=\sum_k [\hat{K}(u_k)\sum_{i=1}^{n}W_{1i}(u_k)\{dN_{1i}(u_k)-\hat{\lambda}(u_k)Y_{1i}(u_k)\}]$ is a constant. So, why is that $1/n\sum_{i=1}^{n}$ over there before this quantity? Am I making any mistake in the numerical integration or in understanding it? Kind regards.


title: Unbiased estimator for variance or Maximum Likelihood Estimator?

text: I know that $\frac{1}{n-1}\sum_{i=1}^{n}(X-\bar{X_{n}})^{2}$ is an unbiased estimator for the variance. I thought of such an estimator to be useful when it is not know from which distribution the data at hand are coming from. Now given I have a data set and I do know which is the underlying distributions of the observations. Shouldn't I prefer a Maximum Likelihood Estimator for the variance even if it is biased since it obviously makes use of the underlying distribution of the data over the general variance estimator cited above, which is not for any specific distribution? Thanks


title: References for stability of Factor Analysis extraction methods

text: Does anyone have any good references comparing suitability / stability of different extraction methods for exploratory factor analysis on binomial variables? * * * I have had a search on-site and internet search but i have only been able to find one 20 year old paper on WLS extraction, which is still good. I do hope for a bit more. I am hoping someone has conducted some simulations, or even better some theory, regarding this due to rather large differences in the factor groupings i get depending on the extraction method used. Testament to the rather ill-conditioned correlations matrices that seem to go hand in hand with tetrachorics. Gung, stability is probably not the best word to use. It is just when using the same correlation matrix as input i get very different answer from ML, WLS and OLS. I do have zero cells, and when i use a correlation matrix with continuity correction some of the differences diminish. But i am surprised about the large differences that occur with different methods on the same input matrix. I do however need to justify the method that i am going to use.


title: Comparing two vectors from negative binomial distribution in R

text: I'm using R and have two vectors of discrete values. They are not strictly speaking categorical because the values themselves are number of dots counted on the image of a cell (whole vector is all the cells on the image). There are two vectors: reference and a vector with dot counts after some perturbation What I believe is that such data should follow negative binomial distribution and some sort of goodness of fit should give a p-value and some statistic describing whether the two distributions differ significantly. What people advised me is that chi square test would do the trick but in my understanding chi square considers all values only as a category and ignores the fact that these are numbers and if lets say number of cells with 5 dots decreased a bit while number of cells with 4 dots increased that's not the same as if same situation would happen with 0-dots and 6-dots categories. However what I didn't find is a test which could deal with negative binomial distributions. I hope I described the problem clearly. So if somebody know any test which would deal with such kind of data or if anybody thinks that my assumptions are wrong you are welcome to share your ideas. **Example 1**               library(ggplot2)     c.dots = c(0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 3, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0,                0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 2, 2, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 2,                0, 0, 0, 1, 0, 0, 1, 0, 0, 3, 0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 1, 1,                0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 1, 0,                0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 2, 0,                0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0,                0, 0, 1, 1, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 1, 0, 0, 1, 1, 0, 0, 3, 0, 0,                0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0)          w.dots = c(0, 0, 0, 1, 3, 1, 1, 1, 1, 0, 0, 2, 0, 0, 2, 1, 0, 1, 3, 0, 1, 0, 0, 0, 2,                0, 2, 2, 0, 3, 1, 2, 1, 0, 2, 1, 0, 2, 0, 1, 2, 1, 0, 0, 1, 0, 1, 1, 0, 0,                0, 1, 1, 0, 2, 0, 0, 1, 3, 0, 0, 1, 0, 2, 1, 0, 1, 1, 1, 1, 1, 1, 2, 1, 1,                2, 4, 1, 0, 0, 2, 2, 0, 1, 0, 1, 3, 0, 2, 1, 1, 2, 0, 0, 0, 0, 0, 0, 1, 0,                1, 1, 0, 1, 0, 0, 2, 0, 1, 0, 2, 1, 0, 1, 2, 0, 4, 2, 0, 1, 0, 2, 0, 1, 2,                1, 1, 2, 1, 1, 3, 1, 0, 1, 0, 1, 2, 0, 1, 2, 0, 1, 1, 2, 2, 0, 3, 0, 1, 1,                0, 0, 2, 0, 1, 1, 0, 1, 2, 0, 0, 1, 0, 1, 2, 0, 0, 4, 3, 0, 1, 0, 0, 1, 0,                4, 0, 1, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 1, 1, 0, 3, 1, 1, 0, 4, 1, 1, 3)          chisq.test(rbind(table(w.dots), table(c.dots)))          nbrand = rnbinom(length(c.dots), mu = 1, size = 1)     ggplot() +        geom_density(aes(x=x), data=data.frame(x=c.dots), fill="red", alpha=0.5) +        geom_density(aes(x=x), data=data.frame(x=w.dots), fill="blue", alpha=0.5) +       geom_density(aes(x=x), data=data.frame(x=nbrand), colour="green", alpha=0, linetype=3)      **Example 2**               library(ggplot2)     c.dots = c(1, 0, 0, 1, 0, 0, 3, 0, 1, 0, 3, 0, 2, 0, 0, 2, 2, 0, 0, 1, 1, 0, 0, 1, 0,                0, 0, 0, 1, 0, 1, 0, 2, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,                1, 1, 1, 2, 0, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 3, 4,                         0, 1, 1, 0, 1, 0, 2, 1, 2, 2, 3, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1,                1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 2, 0, 0, 3, 2, 2, 1, 0, 2, 0, 2, 2, 0, 0, 2,                1, 0, 2, 0, 0, 2, 2, 1, 0, 0, 0, 0, 0, 1, 0, 3, 0, 1, 0, 0, 1, 0, 0, 0, 0,                2, 1, 1, 0, 1, 0, 1, 1, 0, 1, 3, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 2, 1, 0,                1, 2, 0, 0, 3, 3, 0, 1, 2, 0, 0, 1, 1, 0, 1, 1, 3, 1, 3, 0, 2, 0, 0, 0, 0)          w.dots = c(1, 3, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 3, 0, 0, 0, 1, 2, 0,                1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 5, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 2,                0, 1, 0, 3, 0, 0, 1, 2, 3, 1, 0, 0, 0, 2, 1, 1, 2, 0, 2, 0, 3, 0, 2, 0, 0,                0, 0, 2, 0, 1, 0, 2, 0, 0, 1, 1, 2, 3, 0, 2, 2, 1, 0, 1, 0, 0, 1, 0, 1, 0,                0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 2, 0,                0, 1, 2, 1, 1, 1, 2, 1, 2, 3, 2, 0, 0, 0, 0, 0, 1, 2, 0, 0, 1, 0, 0, 1, 1,                0, 0, 1, 2, 1, 0, 1, 2, 1, 1, 1, 0, 1, 0, 0, 0, 2, 1, 1, 1, 0, 0, 0, 0, 0,                1, 2, 1, 2, 0, 1, 2, 1, 0, 1, 3, 2, 1, 0, 0, 0, 0, 0, 2, 1, 1, 2, 2, 1, 2)          chisq.test(rbind(table(w.dots), table(c.dots)))          nbrand = rnbinom(length(c.dots), mu = 1, size = 1)     ggplot() +        geom_density(aes(x=x), data=data.frame(x=c.dots), fill="red", alpha=0.5) +        geom_density(aes(x=x), data=data.frame(x=w.dots), fill="blue", alpha=0.5) +       geom_density(aes(x=x), data=data.frame(x=nbrand), colour="green", alpha=0, linetype=3)


title: What about throwing out 2/3rds of a data set?

text: I have a data set of presidential election results for the past 33 elections. I threw out everything pre 1972, however, because of the structural changes that occurred through the decades, and because I am assuming a degree of auto- correlation in the data, so I want to get rid of previous trends that are not material contemporarily. With the entire data set, I could not get any models with a decent R-squared, the highest was .068, and P value was 15 percent. However, after dropping all of the data up to the year 1972, I was able to fit a different model that had a P value of .028 and also had an R-squared of .45. The current data set consists of 10 points. If I add in 2 more data points, p value is still small, and the R-squared drops to about 30. Is it okay that I dropped this data?


title: Starting point for survival analysis (overall survival and recurrence-free or progression-free survival)

text: In a RCT we would like to compare oncological outcome of two surgical 2-stage strategies, i.e. patients require to undergo both stages to become free of cancer. The problem is that strategy A generally takes less time to come to the second stage becoming free of cancer. On the other hand, in strategy B some patients might not even get to the second stage of their therapy because of the longer time between stage 1 and 2, i.e. in this longer interval they run the risk of disease progression or not getting fit for stage 2. **My question is, from which time-point should our overall survival be calculated?** Starting at randomisation or starting at completion of stage 2. We believe at randomisation. **Should the starting point be different for recurrence-free OR progression- free survival?** **Is it at all possible to compare recurrence-free survival between the two groups as some patients in strategy B might not proceed to their second stage (and become free of disease at all, which is a requirement for being able to assess the event 'recurrence')?**


title: Two unbiased estimators for the same quantity

text: In several situations, I have two unbiased estimators, and I know one of them is better (lower variance) than the other. However, I would like to get as much information as possible, and I would like to do better than throwing out the weaker estimator. $$\newcommand{\Outcome}{\text{Outcome}}\newcommand{\Skill}{\text{Skill}}\newcommand{\Luck}{\text{Luck}}\Outcome = \Skill + \Luck$$ $\Outcome$ is observed. $\Skill$ is what I would like to determine. $\Luck$ is known to have the average value $0$. From other observables, I can estimate $\Luck$ by $L$ so that $\mathbb E(L) = 0$ and $\mathrm{Var}(\Luck-L) \lt \mathrm{Var}(\Luck)$. $\Outcome$ is an unbiased estimator for $\Skill$. A better estimate from variance reduction is $\Outcome - L$, which is also unbiased. For example, in one situation $\Outcome$ is the average of repeated trials, and I might produce a $95\%$ confidence interval of $[-5.0,13.0]$ without using variance reduction. Using variance reduction, I might get a confidence interval of $[-2.0,4.0]$. The typical practice is for people to use $\Outcome-L$ instead of $\Outcome$. However, this is unsatisfactory to me because in my experience, there is more information in the pair $(\Outcome, \Outcome-L)$ than in just $\Outcome-L$. Specifically, in some situations I know that if $\Outcome$ is low, then $\Outcome-L$ tends to be an underestimate for $\Skill$, and if $\Outcome$ is high, then $\Outcome-L$ tends to be an overestimate for $\Skill$. **What's a good way to take advantage of the extra information from knowing both estimators?**


title: Simulate data for path model

text: How does one simulate data for a path model? Below is an example path model with parameters b1, b2, b3, b4, b5, e1, e2 and e3. I would like to investigate how sample size affects my estimates and what happens to the estimates if a path is added or removed. ![enter image description here](http://i.stack.imgur.com/Sj5t7.png)


title: Kolmogorov's paper defining Bayesian sufficiency

text: I'm looking for a translation to either English, French or German of Kolmogorov's Russian paper _Kolmogorov, A. (1942). Sur l’estimation statistique des paramètres de la loi de Gauss. Bull. Acad. Sci. URSS Ser. Math. 6, 3–32._ This seminal, oft-cited paper (e.g. in Lehmann and Romano's Testing Statistical Hypotheses (2008) p. 20 and in Blackwell and Ramamoorthi's A Bayes But Not Classically Sufficient Statistic (1982) p. 1) originated the concept of Bayesian Sufficiency. Translations to other Romance languages or to Dutch are also welcome.


title: Cholesky update for removing row

text: What is most efficient way to updating cholesky factorization of for removing a column from the matrix?               T = Chol(X' *X)      If I remove a column from X, how to effienciently update cholesky. It is trivial to remove the last column of X as new cholesky will be removing last row and column form T. What about removing other column?


title: How do I interpret the figure output from package dlnm in R?

text: I am performing distributed non-linear lag models in R. I got the figure result of `dlnm` as shown in the vignette (pdf) on page 13: ![enter image description here](http://i.stack.imgur.com/Qjf7B.png) The X-axis is lag, which I can understand. However, I cannot understand what the label for the Y-axis means. For the figure on page 13, the Y-axis is called "RR" (relative risk) which is positive; however, I got a negative RR when I input my data. My $X$ is temperature, my $Y$ is pollutant concentration. Can I still use the model? Could anybody explain it to me a little bit? The R code that I am using is:               OrigData <- read.csv("Hourly data for use SHULE-2.csv", na.strings = ".",                           header=T, sep=",")     CO2H6 = OrigData$CO2H6      # $     OutT  = OrigData$OutT       # $     Year  = OrigData$Year       # $     Month = OrigData$Month          argvar <- list(type="ns", df=11)     arglag <- list(type="ns", df=5)          suppressWarnings(cb1 <-crossbasis(OutT, lag=12, argvar=argvar, arglag=arglag))     model1 <- glm(CO2H6~ cb1 +as.factor(Month))     summary(model1)          cp <- crosspred(cb1,model1,-22.83:34.79,by=1, cumul=TRUE)     par(cex.axis=0.7,cex.lab=0.7)     d3 <- plot(cp, xlab="Delay Effects of OutT in IN2H-H6", zlab="CO2 con.", r=90,                 d=0.3, col="red" , xlab="Out T")          plot(cp, "slices", var=10, ci="n", ylim= c(-1000,1000), ylab="CO2 con.", lwd=1.0,           xlim=c(0,10))     for(i in 4:6){       lines(cp, "slices", var=c(-20,-10,0,10,20,30)[i], col=i+1, lwd=1.5)     }     legend(7,1000, col=5:7, lwd=0.9, pch=1:1, legend=c("OutT =10℃", "OutT =20℃",                                                         "OutT =30℃"), cex=0.7)      ![enter image description here](http://i.stack.imgur.com/ZBJ8K.jpg) Since RR cannot be negative, could there be a way for me to transfer my x(Temperature) and y (Pollutant con.) to make it fitting the requirement of the model? Thanks for all your help!!


title: Statistics on three variables/factors

text: I'm trying to say something about sales of candy bars. I got data on sales for a population of children. There are three factors: age, income of parents, and name of a child. What I can compute from this data is, say, that 34% of all sold candy bars were bought by children of age 12, 25% by children of age 10, etc. At the same time, I can compute that 52% were bought by children coming from families with monthly income X, 35% by children coming from families with monthly income Y, etc. And finally, I know that say 40% of the bars were bought by Sams, 25% by Abigails, etc. What I would like to deduce from these three characteristics is, given a child of certain age and name, and coming from a family with monthly income Z, what is the "probability" that he/she will buy a candy bar?


title: When should I prefer kappa over accuracy for binary classification problems?

text: I'm using the `caret` package in `R` that provides two kind of metrics for classification problems - accuracy and kappa. When should I prefer one over another?


title: Glmnet and dummy variable trap

text: When feeding a categorical variable into `glmnet` do I code `n` or `n-1` dummy variables? For instance if using days of the week as an independent variable would I use 6 dummies or 7? If the answer is 6, how do I interpret coefficients, etc for dropped category? EDIT: Here's some example code:               library(glmnet)     library(caret)          df1 <- data.frame(id = 1:210, var1 = rep(c('Mon','Tues','Wed','Thurs','Fri','Sat','Sun'),30))     df1$targetVar <- runif(210)         df1$mktVol <- round(runif(210)*1000000,0)     df1$mktVol <- ifelse(df1$var1 %in% c('Sat','Sun'), 0, df1$mktVol)          df1     vtu <- c('mktVol','var1')     dv1 <- dummyVars( ~.,data = df1[,vtu])     df2 <- data.frame(predict(dv1,df1))          glmnet1 <- cv.glmnet(df2$targetVar, data.matrix(df2[,-c('targetVar')]), nfolds = 5)          glmnet1 <- cv.glmnet( data.matrix(df2[,-1]), df2[,"mktVol"] ,           family="gaussian", alpha=.95, nfolds=5, standardize = FALSE,           type.measure="mse")          Coefficients1 <- coef(glmnet1, s = glmnet1$lambda.min)     Active.Index <- which(Coefficients != 0)     Active.Coefficients <- Coefficients[Active.Index]     names(X1)[varsToUse[Active.Index]]          ##############################          df1 <- data.frame(id = 1:210, var1 = rep(c('Mon','Tues','Wed','Thurs','Fri','Sat','Sun'),30))     df1$targetVar <- runif(210)         df1$mktVol <- round(runif(210)*1000000,0)     df1$mktVol <- ifelse(df1$var1 %in% c('Sat','Sun'), 0, df1$mktVol)          df1     vtu <- c('mktVol','var1')     #dv1 <- dummyVars( ~.,data = df1[,vtu])     #df2 <- data.frame(predict(dv1,df1))     dv1 <- model.matrix(~.,data = df1[,vtu])          #glmnet1 <- cv.glmnet(df2$targetVar, data.matrix(df2[,-c('targetVar')]), nfolds = 5)          glmnet1 <- cv.glmnet( data.matrix(df2[,-1]), df2[,"mktVol"] ,           family="gaussian", alpha=.95, nfolds=5, standardize = FALSE,           type.measure="mse")          Coefficients2 <- coef(glmnet1, s = glmnet1$lambda.min)          ##############################          df1 <- data.frame(id = 1:210, var1 = rep(c('Mon','Tues','Wed','Thurs','Fri','Sat','Sun'),30))     df1$targetVar <- runif(210)         df1$mktVol <- round(runif(210)*1000000,0)     df1$mktVol <- ifelse(df1$var1 %in% c('Sat','Sun'), 0, df1$mktVol)          df1     vtu <- c('mktVol','var1')     #dv1 <- dummyVars( ~.,data = df1[,vtu])     #df2 <- data.frame(predict(dv1,df1))     dv1 <- model.matrix(~ 0+ .,data = df1[,vtu])          #glmnet1 <- cv.glmnet(df2$targetVar, data.matrix(df2[,-c('targetVar')]), nfolds = 5)          glmnet1 <- cv.glmnet( data.matrix(df2[,-1]), df2[,"mktVol"] ,           family="gaussian", alpha=.95, nfolds=5, standardize = FALSE,           type.measure="mse")          Coefficients3 <- coef(glmnet1, s = glmnet1$lambda.min)     Coefficients1      Coefficients2      Coefficients3


title: Do the commands logit and glm (Binomial family) in Stata use different fitting/estimation algorithm?

text: In my mind, the commands `logit` and `glm` (Binomial family) in `Stata` both use maximum likelihood estimation, and the ML estimates and SE values are the same. But one instructor stated that they use different fitting/estimation algorithm. So I would like to hear the difference.


title: Excel T-Test “Set X is better than Set Y by ___%”

text: Coming from this post, but their question wasn't fully answered. I understand how the t-test works – if I have two sets of data and want to see if one is better (or worse) than the other, I would say `ttest(setx, sety, 1, 2)` for a type 2, and that would tell me whether or not one set is better/worse than the other. What I don't understand is how to say "Set `X` is better than set `Y` by `__`% with statistical significance." My best guess would be using `ttest(setx, sety*.2, 1, 2)`, and if that still gives me a small enough p-value, then that .2 would turn into "Set `X` is better than set `Y` by 80%." Am I thinking that through correctly?


title: Which is more accurate glm or glmnet?

text: R glm and glmnet use different algorithms. I notice non trivial differences between the estimated coefficients when I use both. I am interested in when one is more accurate than another, and the time to solve/accuracy trade off. Specifically I am referring to the case where one sets lambda=0 in glmnet s.t. it is estimating the same thing as glm.


title: interaction terms in separate models?

text: I'm doing a multilevel mixed model with three cross-level interactions in LMER. Do all three interactions necessarily need to be included in the same model or can I evaluate them in three separate models?


title: Compare fit of ordinal to linear model

text: I´m doing multivariate modeling on a dataset where the dependent variable is ordinal-values running from 0-11. So far I´ve done general linear models and it runs alright. However, I´d also like to run it as an logistic ordinal regression.   * What is the recommended maximum number of categories in the dependent ordinal variable?    * Would I need to cut the variable into groups and reduce the number of possible categories?   * How can I compare the fit of an ordinal-outcome model to a linear regression model? Especially, if the ordinal variable has been cut into fewer groups?


title: Time series forecast by Principal Component Analysis

text: Suppose that I have a series of $M$ time-observations of $N$ "quantities" $z_1(t_1),...,z_1(t_M)$, ..., $z_N(t_1),...,z_N(t_M)$. I want to estimate the values of $z_1(t_{M+1}),...,z_N(t_{M+1})$. This a problem of interest, for example, in stock asset prediction. I want to use Principal Component Analysis (PCA) performed by a Singular Value Decomposition (SVD). My questions:   1. What is the physical meaning of the first singular vector (namely, that one corresponding to the largest singular value)? What is the physical meaning of the remaining singular vectors? I understand that the singular vectors provide uncorrelated linear combinations of the above random variables.   2. What is the physical meaning of the singular values? As long as I know, they are related to the variances associated to the singular vectors. But then, does it mean that the most probable value of the quantities at time $t_{M+1}$ is due to the least singular values (least variances)?   3. How do I use the singular values and vectors of PCA to predict the value of the quantitites at time $t_{M+1}$? Could you let me understand the idea behind? Thank you very much in advance.


title: Time-series data pre-aggregated into non-stationary rolling 12-month periods: are there special considerations for modeling?

text: I'm exploring the use of changepoint detection or other methods (am slowly becoming aware of wavelet transformation, etc. but have tons to learn in this area) to identify key shifts in health care performance patterns over time. However, many of the metrics I'm seeking to analyze (e.g., health care quality metrics) are both generally calculated and more reasonably interpreted as rolling-12 month aggregates. For example, it's important to me to track on a monthly basis the proportion of patients who are up-to-date on a certain lab test within the 12-month period ending that month, but I'm not particularly concerned with how many of these tests occurred specifically in August. So it's something sort of like having a moving average to work with as a raw starting point. That said, there are also reasons why this rolling-12 aggregation does not result in a stationary process either. My thought was to account for the data structure and seasonality by modeling it as a function of a 1-month and 12-month lag. Is this the proper way to think about this data? Is there anything else or a better approach I should be doing/considering? Again, my general goal is surveillance of the general trend as well as breaks -- so if it affects the answer, I'm looking at this in the context of using the R strucchange package, CUSUM statistics, or some other approach to identify good and bad anomalies. Thanks, Shelby


title: A problem with two ways of calculating AIC

text: I find that there are two ways to calculate AIC: AIC = -2ln(likelihood)+ 2K and   AIC = n*ln(RSS/n)+2K I have:               crf <- c(0.3333333, 0.5000000, 0.6666667, 0.6666667, 0.6666667, 0.8333333, 0.1666667, 0.3333333, 0.5000000, 0.5000000, 0.8333333, 0.5000000,0.6666667, 0.5000000, 0.6666667, 1.0000000)          co3 <- c(218.20, 243.84, 267.97, 286.31, 315.01, 315.01, 241.09, 242.52, 243.84, 245.04, 246.10, 284.15, 285.79, 287.31, 288.67, 289.49)          n <- length(crf)     model <- lm(crf ~ co3) #y=ax+b     aic1 <- AIC(model, k=2)     RSS <- sum(resid(model)^2)     aic2 <- n*log(RSS/n) + 2*3      I have `aic1=-7.220462` and `aic2= -52.62649` but if I will do like this:               aic3 <- n + n*log(2*pi) + n*log(RSS/n) + 2*3      I have `aic3=-7.22046`. I don't understand why, everywhere it is written that AIC = n*ln(RSS/n)+2K, what is this term n + n*log(2*pi)?


title: Fieller's Theorem Degrees of Freedom

text: In a previous question I was trying to figure out how to compute the confidence interval of a ratio. Having come to the conclusion that the correct way to do this is Fieller's Theorem, I realize I'm very unsure about my interpretation of the degrees of freedom. I'm trying to compute confidence intervals for a speedup -- computed as a ratio: $$ \frac{\text{time mean}_{\text{1 thread}}}{\text{time mean}_{n \text{ threads}}} $$ Each thread count has $k$ samples. I want to say I need 1 degree of freedom here... but I also have been known to be wrong in these things before...


title: A simple question in regards to a academic article

text: **Introduction:** I am currently doing some reading on psycological research and in this context I have some questions which may seem very simple. Nevertheless, it is not very transparent from the articles how the calculations are made so please feel free to answer as broad or as specific as you feel. **Question:** What confuses me is a statement like this: b= -.09, t(82) = -2.14, p<.05 or this: b= .34, t(68.87) = 2.03, p<.05 Should I consider this as the General Linear Model? Any input is appreciated. Thanks! See this example


title: In the survival analysis package of R, `survival`, where do we specify the hazard functional form?

text: I am trying to fit either the `coxph` or `survreg` function with the `Survival` package. I would like to specify my Hazard functional form as: $H(t) = exp\{ (\beta_0) + (\beta_1)*X_t +(\beta_2)*log(X_t)\}$. Is the way to specify it into my model as just straight insertion into the formula portion of `coxph`? As in:               coxph(Surv(time) ~ X_t + log(X_t), data = data.frame)


title: Negative binomial - likelihood ratio and likelihood ratio chi-squared

text: I ran a hierarchical negative binomial regression analysis, and got information relative to the log likelihood and likelihood ratio chi-squared in the output. I have the following questions regarding the goodness of fit:   1. What model is preferred? I heard the one with the log likelihood closer to zero in absolute form; is this correct?   2. Is the Likelihood ratio chi-squared the same as what some call the Wald chi-squared, to test the overall model, compared to previous ones? (i.e. null model)   3. Are there any "typical" values for the Log likelihood and the likelihood ratio chi-squared? I got these results for my 3 models:  $$ \begin{array}{|c|c|c|c|} \hline & \mbox{Model 1} & \mbox{Model 2} & \mbox{Model 3} \\ \hline \mbox{Log likelihood} & -3517.05 & -3278.52 & -3275.0 \\ \mbox{Likelihood ratio chi-squared} & 3727.3 & 4331.07 & 4338.3 \\ \hline \end{array} $$


title: How to prepare variance matrix in Statistica (or different tool)?

text: I have data which describe measurments in four different point for few days. Four water parameters have been measured. In other words, data contains following columns:   1. point name - we have four points names, of course this value is not unique for each row (as we have few measurment for few days)   2. First parameter value   3. Second parameter value   4. Third parameter value What I want to achieve is matrix with 16 rows and 16 columns (4 parameters multiplied by 4 points results in 16). On the diagonal of this matrix I will have ones. I know how to receive matrix four x four, but I haven't succeed in 16 x 16. Could you guide me how I can do that in Statistica or other modern statistical package (for example R!), please?


title: Assessing effect of intervention where aptitude tests used in pretest and posttest have some overlapping and some distinct items

text: In my research there is a pretest, and a treatment followed by a posttest. The pre and the posttests are similar but not exactly the same. A few questions were changed between the pre and posttests to avoid sensitizing the students to the questions in the pre-test, which might affect the internal validity of the study. Because I changed the post test, I got into the instrumentation threat, because now the pre and post tests are not exactly the same. In order to eliminate the instrumentation threat, I plan to compare the pre and the posttests taking only the common questions from both tests.   * Is this permissible statistically?    * Is there any literature that supports or talks about unequivalent pre and post tests?


title: Understanding the variance of random effects in lmer() models

text: I'm having trouble understanding the output of my `lmer()` model. It is a simple model of an outcome variable (Support) with varying State intercepts / State random effects:               mlm1 <- lmer(Support ~ (1 | State))      The results of `summary(mlm1)` are:               Linear mixed model fit by REML      Formula: Support ~ (1 | State)         AIC   BIC logLik deviance REMLdev      12088 12107  -6041    12076   12082     Random effects:      Groups   Name        Variance  Std.Dev.      State    (Intercept) 0.0063695 0.079809      Residual             1.1114756 1.054265     Number of obs: 4097, groups: State, 48          Fixed effects:                 Estimate Std. Error t value     (Intercept)  0.13218    0.02159   6.123      I take it that the variance of the varying state intercepts / random effects is `0.0063695`. But when I extract the vector of these state random effects and calculate the variance               var(ranef(mlm1)$State)      The result is: `0.001800869`, considerably smaller than the variance reported by `summary()`. As far as I understand it, the model I have specified can be written: $y_i = \alpha_0 + \alpha_s + \epsilon_i, \text{ for } i = \{1, 2, ..., 4097\}$ $\alpha_s \sim N(0, \sigma^2_\alpha), \text{ for } s = \{1, 2, ..., 48\}$ If this is correct, then the variance of the random effects ($\alpha_s$) should be $\sigma^2_\alpha$. Yet these are not actually equivalent in my `lmer()` fit.


title: Multiple time-series -- measuring medication weight gain

text: First - very new to Statistics; about half through a basic biostats book and an R book. I have a set of data where I"m trying to see if there is a correlation between a medication and weight gain. The dataset:   1. ~2200 patients   2. A list of their primary diagnoses, weight, date's of visit, and whether they were on medication at that time.    3. Each patient may be on or off the medication at various points, though I have limited the data set to patients receiving at least 90 days of medication.    4. Different patients are on / off the medication for varying # of days. Also, each patient has recorded visit for a varying amount of total days (though at least 90 based on 90 day min. medication limit from #2).  ![General Outlay](http://i.stack.imgur.com/Tx2B2.jpg) My main goals are:   1. Is this medication correlated with weight gain?    2. If so / not - does it depend on a certain threshold length of constant use?    3. (Eventually) do either of the above two correlations change based on the patients diagnoses (Some pt's will be grouped together)  My question is - really, where should I begin? What types of analyses should I be doing? I'm willing to put in the reading / work to figure it out; but I'm not quite sure where to look. As stated in intro, I have - some - experience with R, and am expecting to carry out my in depth analysis in that. UPDATE: So far, I have been looking at Average of change in weight / days on medication; so basically giving each patient a Wt. Change per day of med. Im planning on adding wt. change per day while off medication and all patient variables, but this illuminates what I need something more advanced for -- How do I account for trends occurring independent of medication. A thought - say a patient started a diet, and has been steadily losing 1lb per week prior to medicaiton. If that medication is started, it is possible the overall trend will overshadow the medication effect. In this case - is there a practical way to account for long term trends, without examining each patients data individually?


title: QQ plot with one dataset

text: I have been asked to "construct a (normal) QQ plot" for a set of numbers. However it is not a set of pairs or a pair of sets, just a single set of numbers. If we go ahead and assume the second probability distribution is the normal distribution, how would I go about satisfying this request?


title: SOM clustering for nominal/circular variables

text: Just wondering if anyone is familiar with clustering nominal inputs. I've been looking at SOM as a solution but apparently it only works with numerical features. Are there any extensions for categorical features? Specifically I was wondering about 'Days of the Week' as a possible features. Of course it is possible to convert it into a numerical feature (i.e. Mon - Sun corresponding to nos 1-7) however then the Euclidean distance between Sun and Mon (1&7) would not be the same as the distance from Mon to Tues (1&2). Any suggestions or ideas would be much appreciated.


title: How to design a study and test the effect of dosage level on cure probability?

text: My goal is to quantify dosage impact on cure probability for different patients. Let's suppose I have N patients with their charateristics such as age, gender, weight.... Also let's assume there is a known range for dosage, say from 0 to 1. Here is what I'm currently doing:   1. Randomly split all patients into 3 groups   2. Group 1 is given dosage of 0.1, group 2 dosage of 0.5 and group 3 dosage of 0.9. Please note size of group is different. Presumably group 1 is largest because I need decent number of cured patients and I know low dosage leads to low cure probability. Also note that I can not test more dosage because N is limited.   3. Get outcome (0/1) for all patients   4. Train logistic regression model for each group   5. If I want to get cure probability for any particular patient with untested dosage (say 0.7), I just calculate 3 probabilities based on models in step 4 and do some simple interpolation. My question is if there is any better way to do this. I'm not feeling very comfortable about the interpolation part.


title: Van Westendorp price modeling in R

text: I want to do price modeling using the Van Westendorp price sensitivity meter. Can anybody provide me resources describing how to do this using R?


title: Any suggestion for clustering method for unknown number of clusters and non-Euclidean distance?

text: I need some suggestion for clustering (unsupervised classification) method for a consulting project. I am looking for a method that hopefully has the following properties:   1. The subject of my study has three properties. One is represented by a (non-Euclidean) distance matrix and the other two are in the form of vectors in Euclidean space. The distance matrix comes from sequences and can be in the form of percent of dissimilarity or other measurement of distance of sequences. The algorithm should be able to take both vectors in euclidean space and non-euclidean distance as input. For example, K-medoids can work with a distance matrix but K-means can not.   2. I would like the algorithm to select the number of clusters and the weight for three properties automatically (with prior knowledge and constraint).   3. I have information of previously identified “centers of clusters”. I would like to incorporate it as prior or initial values.   4. As a statistician, I would prefer the method to have a clear likelihood or loss function. The closest thing I can think of is fitting a mixture model in Bayesian framework using reverse jump MCMC to determine the number of clusters. The vectors in R^d can be easily formulated into a normal likelihood but how to deal with the distance matrix is unclear to me. I can restrict the mean of normal likelihood to be at each of the observation of get the MCMC running but that does not have a clear mathematical / statistical meaning. Does anyone have experience with a similar problem? Suggestion to references will be highly appreciated!


title: Combining metrics for score card

text: We are building a score card system to measure the performance of various groups in the organization. The scores are generated by true / false metrics on individual transactions that each group completes. Transactions are qualified for measurement by some metrics but not others. The number of transactions also varies between groups. Here is a sample of what the data looks like:               Group  Transaction   Metric1Score  Metric1Qualified   Metric2Score  Metric2Qualified  ...      G1     Tran1         1             1                  0             1     G1     Tran2         1             1                  1             1     G2     Tran3         0             0                  0             1     ...      We'd like to combine all of these transaction metrics into one score per group and be able to compare scores between groups. Does anyone know how to get started on this or have any places to refer to?


title: Analysis of negative binomial distributions using glm.nb() in R

text: I have observations of four different groups of people. The dependent variable are count data (medical emergencies in the past 2 months). For each group, the dependent variable follows a negative binomial distribution with a maximum at 0. Now I would like to examine, if the factor "group" (meant to be categorical) is a significant factor for the number of reported medical emergencies. What I have done so far: I conducted a glm using `glm.nb` in R and examined the result of `glm.nb` using `anova()`. I would really appreciate, if someone could confirm for me, if this procedure is feasible for a negative binomial distributed response variable and a categorical factor. My second question: `anova()` with the result of the `glm.nb` produces the warning: `tests made without re-estimating 'theta'`. As far as I understand, theta is a dispersion parameter for the negative binomial distribution. However, the distributions of observations in each of the four groups have very different dispersions. Does R calculate a "mean" theta for all four distributions?


title: Why does the scikit-learn bootstrap function resample the test set?

text: When using bootstrapping for model evaluation, I always thought the out-of-bag samples were directly used as a test set. However, this appears not to be the case for the scikit-learn bootstrap function, which seems to build the test set from drawing with replacement from the out-of-bag data subset. What is the statistical reasoning behind this? Are there specific scenarios where this technique is better than just evaluating on the out-of-bag-sample or vice versa?


title: Efficient calculation of matrix inverse in R

text: I need to calculate matrix inverse and have been using `solve` function. While it works well on small matrices, `solve` tends to be very slow on large matrices. I was wondering if there is any other function or combination of functions (through SVD, QR, LU, or other decomposition functions) that can give me faster results.


title: How to estimate missing data points from a moving average curve?

text: The table below is extracted from Wikipedia. It shows the Pollution Standards Index (PSI) readings of Singapore. ![enter image description here](http://i.stack.imgur.com/Qgmmv.png) I am trying to estimate the four missing data points (2-5am) on day **20th June**. I did this by first plotting a graph of the 3-hourly readings. Then I try to derive the 1-hourly readings with the formula:               PSI(1hrly)n = 3 * PSI(3hrly)n - PSI(1hrly)n-1 - PSI(1hrly)n-2      This derived data is then used to plot a graph of the 1-hourly readings. In theory, the 3-hourly curve functions like a moving average curve that lags the 1-hourly curve. There shouldn't be too much deviation. I tried a few values but I am getting wildly erratic results. The curve below is plotted based on values 200,200,140 and 200. ![enter image description here](http://i.stack.imgur.com/cB9oM.png) Is there a way to get a rough estimate of the missing data?


title: Contours containing a given fraction of $(x,y)$ points

text: I often need make $(x,y)$ scatter plots that have many ($>10^5$) points. I've experimented with different ways of representing this many points that capture the distribution while getting around the messiness of actually putting a million points in an image. I've done the obvious things like thinning points, and plotting the local density of points as either contours or a colour map. I've found that thinning tends to be ineffective as the distributions I deal with are often exponentially peaked, so I just get the same problem in a narrower area, and local density representations are sensitive to the details of smoothing/sampling used to get the density field. Lately I've been thinking about trying a different strategy, but I haven't been able to come up with an algorithm to make it work (which is where this question comes in). This sort of plot, which I've heard cosmologists colloquially call a 'banana plot', is the inspiration, and I think is fairly common in Bayesian analysis: ![likelihood contours](http://inspirehep.net/record/1082797/files/FINAL_FIGURES_b11_spt_lcdm.png) Here the contours show the $1\sigma$ and $2\sigma$ confidence intervals for a combination of two parameters. In this application a 2D probability density function is generated, so making the plot is easy. I want to do something similar. I'd like to draw contours enclosing 99%, 95% and 68% (and more generally any percentage) of the data points in an $(x,y)$ dataset. Obviously such a contour is not unique, so an extra constraint is needed. I think something like minimising the area or maximising the average density of points inside the contour should give the desired effect. Does anyone know an algorithm that would produce such a contour? I eventually want to make a plot using matplotlib, so bonus points if you happen to know such an algorithm that is implemented in python, but I'll do the implementation myself if I have to. Not terribly familiar with the tags here on CV.SE, so if I'm missing any obvious ones please suggest/edit them in.


title: Classifying IP addresses with a decision tree

text: I am working on a classification problem using IP addresses as input and I am trying to find which IP addresses or subnets are likely to belong to a spammer. I have data consisting of the four octets in the IP address as well as a label: SPAM or OK. This feels like a good problem for a decision tree, but I am finding that most decision tree algorithms consider the order of the variables interchangeable. For example, `ctree` from `party` might output a rule interpreted as `IF octet1 == 192 && octet4 == 190 THEN label => SPAM` but with IP addresses, the order of the variables matters. That is, octet 2 must be considered after octet 1, and octet 3 must be considered after octet 2 and so on to yield rules like `IF octet1 == 206 THEN label => OK` `IF octet1 == 193 && octet2 == 64 && octet3 == 11 THEN label => SPAM` What type of decision tree model is this and what algorithm/tool could help here? Is there some sort of variation I can use? I prefer to stick with R.


title: Copula with 4 Poisson and 1 Bernoulli margins in R

text: I want to joint 4 Poisson variables and 1 Bernoulli variable with D-vine copula. How can I do this with R?


title: Inconsistency in output from my implementation of noncentral t CDF and R's pt()

text: I am trying to implement a noncentral _t_ CDF as expressed by Guenther (1978), Lenth (1989), and the Wikipedia article on the non-central _t_ in R. I have got my algorithm half working: when the signs of _t_ and Delta are equal, I get identical results to R's pt(), but when the signs are opposite, my results differ non-trivially. Can you help me figure out where I am going off? Here (minus argument validation bits) is my code:               PNCT <- function(t,Delta,df) {         tolerance <- 1e-12         Ftilde <- function(t,Delta,df) {             Delta <- sign(t)*Delta #switch sign if t<0             lambda <- Delta * Delta             j <- 0             SUM <- 0             SUMstep <- 1             y <- (t^2)/((t^2)+df)             while (SUMstep > tolerance) {                 pj <- (1/factorial(j)) * exp(-0.5*lambda) * ((0.5*lambda)^j)                 qj <- ( Delta/( sqrt(2)*gamma(j + 1.5) ) ) * exp(-0.5*(lambda)) * ((0.5*lambda)^j)                 SUMstep <- (pj*pbeta(y, j+0.5, df/2)) + (qj*pbeta(y, j+1, df/2))                 SUM <- SUM + SUMstep                 j <- j+1                 }             return(pnorm(-Delta) + (1/2)*SUM)             }         if (t >= 0) return(Ftilde(t,Delta,df))         else return(1 - Ftilde(t,Delta,df))         }      And here are examples showing the half-right performance of my attempt:               PNCT(t=2.01505, Delta=2.20582, df=5)     [1] 0.3999992     pt(q=2.01505, ncp=2.20582, df=5)     [1] 0.3999992          PNCT(t=-2.01505, Delta=2.20582, df=5)     [1] 0.0339698     pt(q=-2.01505, ncp=2.20582, df=5)     [1] 0.000155399          PNCT(t=2.01505, Delta=-2.20582, df=5)     [1] 0.9660302     pt(q=2.01505, ncp=-2.20582, df=5)     [1] 0.9998446          PNCT(t=-2.01505, Delta=-2.20582, df=5)     [1] 0.6000008     pt(q=-2.01505, ncp=-2.20582, df=5)     [1] 0.6000008      PS I have elsewhere attempted to implement Lenth's algorithm, both by following his published FORTRAN code, and by attempting to replicate implementations in other languages, including R's internal C version, but with even less success. For this question, I am interested in how I have gaffed the above code.


title: Using non-stationary time series data in OLS regression

text: I am using 1983-2008 annual data to test if both gini coefficients and gross national saving in China and the US can affect the US current account balance. The data seem to be non-stationary, but I am a beginner and only knows the basic multiple regression model and autoregressive distributed lag model, can I still use them to these model to test the effects? I know the models would be biased and not accurate, but do they give some useful information. My chosen control variables are real GDP, interest rate, dollar index and maybe some other national income components.


title: Multi Class vs 2 class Naive Bayes

text: I was wondering what are the implications of using a **multi-class Naive Bayes** versus a **2 class Naive Bayes** (for one against everything). Which technique performs **better**? I've previously came across an instance where **2 class LDA classifier** performed **better** than a **multi-class LDA**. Hence, this question. I'm currently using **Scikit-Learn** machine learning package in python to use Naive Bayes classifier.


title: Degrees of freedom correction to pooled regression of demeaned data in R

text: Is there a canonical way to do degrees of freedom (DF) correction to `ols` objects from the `rms` package (and to `lm` objects in general)? I am demeaning my data then running a pooled regression instead of using indicator variables and realizing there's a lot that I don't understand.               library(rms)     library(lmtest)     library(sandwich)          # generate data     set.seed(2001)     dat1 <- data.frame(firm=rep(letters[1:5], each=5),                        year=rep(1:5, times=5) + 1980,                        x1=rnorm(5*5),                        x2=rnorm(5*5),                        x3=rnorm(5*5),                        x4=rnorm(5*5),                        x5=rnorm(5*5),                        stringsAsFactors=FALSE)          # demean on "firm" to obtain within estimator     myDemean <- function(x, use=NULL, skip=NULL) {         innerFun <- function(x, y) {             if (is.numeric(x)) {                 x - ave(x, y)             } else {                 x             }         }         notUse <- !(names(x) %in% union(use, skip))         x[notUse] <- lapply(x[notUse], FUN=innerFun, y=x[use])         x     }     dat2 <- myDemean(dat1, use="firm", skip="year")          # fit models     model1 <- ols(x1 ~ x2 + x3 + x4 + x5 + factor(firm), data=dat1)     model2 <- ols(x1 ~ x2 + x3 + x4 + x5, data=dat2)     coeftest(model1)     coeftest(model2) # wrong SE without df correction          # adjust DF in `coeftest`     oldDF <- model2$df.residual         newDF <- oldDF - (length(unique(dat1$firm)) - 1)     coeftest(model2, df=newDF) # still wrong SE          # adjust DF in model object     model2$df.residual <- newDF         coeftest(model2) # still wrong SE          model2$df.residual # although object updated to DF=16          # can get correct SE with adjusted Var-Covar matrix     newCV <- model2$var * oldDF / newDF         newCV[2:5, 2:5]         model1$var[2:5, 2:5]     coeftest(model2, vcov=newCV, df=newDF)      coeftest(model2, vcov=newCV) # but same p-values with/without DF correction      Is there a "correct" way to make the DF correction? Or should I just do the last one and stop worrying?


title: GLM with unequal sample size

text: I am trying to figure out how to run a GLM with a poisson distribution in R. The data has 4 treatments: $\text{A}(n=16)$, $\text{B}(n=17)$, $\text{C}(n=16)$, $\text{D}(n=20)$, and 2 time periods. This data is count/area so there are decimals. Ultimately, I would like to know if there are differences between the treatments. What is the appropriate way to format these data so that I can run an ANOVA using GLMs with a poisson distribution in R? Below is an example of the data:               time1_A = c(0,0,0,0,0,0,0,0,39.33,0,0,0,76.22,0,0,0)     time1_B = c(0,0,0,0,0,0,0,0,0,0,0,0,39.82,0,0,0,0)     time1_C = c(0,0,63.22,0,0,0,0,0,39.33,0,0,0,45.22,0,0,0)     time1_D = c(0,0,0,0,0,0,0,0,79.83,0,0,0,0,0,0,0,0,0,0,0)     time2_A = c(0,0,39.33,0,0,39.33,0,0,0,0,39.33,39.33,79.62,39.33,0,39.33)     time2_B = c(0,0,0,0,39.33,0,39.33,39.33,0,11.94,0,0,39.33,39.33,79.62,0,15.92)     time2_C = c(39.33,0,39.33,0,199.6,119.4,79.6,0,39.33,318.7,199.6,0,39.33,278.6,278.6,238.5)     time2_D = c(0,0,0,0,0,0,0,0,0,0,39.33,0,0,0,0,0,0,0,0,0)


title: effect of parameters on ARMA process

text: > **Possible Duplicate:**   >  Effect of moving average parameter on variability and variance of demand What is the effect of Theta on MA(1) process, when we move from zero to 1 and from zero to -1, how the process will be affected? the same question for AR(1)


title: Two or more time series. What is the best way to test whether one of them is leading and by what time period?

text: I am trying to prove/disprove that one time series is leading trend for the other ones. Two time series are (probably) independent and the movements are caused by some (let's assume unknown) common factors. What method would be most appropriate? I also want to find the leading period? Data is quarterly.


title: Comparing smoothing splines vs loess for smoothing?

text: I wish to better understand the pros/cons for using either loess or a smoothing splines for smoothing some curve. Another variation of my question is if there is a way to construct a smoothing spline in a way that will yield the same results as using loess. Any reference or insight are welcomed. Thanks.


title: Factoring of conditional probability

text: I have been watching Tom Mitchell's lecture on Bayes Nets: http://cc- web.isri.cmu.edu/CourseCast/Viewer/Default.aspx?id=bc507778-7a18-4121-b345- 83d9bab72f55 He states the following by way of the chain rule. I'm being dense I am sure but don't recognize this quickly. Is it a straight application of the chain rule? $$P(x_3,x_4\ |\ Y) = P(x_3\ |\ Y, X_4 )P(X_4\ |\ Y)$$


title: Calculating the $p$-value of an $F$- statistic

text: I am trying to implement an algorithm for calculating $p$-values of $F$-tests and I need this method to be highly precise. It is easy to implement this with $z$- or $t$- tables, however I don't know how to do this with $F$-values. I have seen some online calculators that do this job, what method do they use? Interpolation of couple of $F$-tables would surely be unprecise. Any help would be appreciated. Thanks in advance


title: Stata's tsset in panel with two time dimensions

text: I have a panel of bond spreads. Spreads are from a specific issuer (firm) and maturity (eg 2, 4, etc years). The problem is that Stata's tsset (or xtreg) command takes two variables to set the data structure. In my case no two variables uniquely identifies each observation. If I do: tsset issuer_id Date I get the error: repeated time values within panel My data looks like the following:               Issuer   Date   Spread Maturity IndependtVar1 IndependentVar2 ...     Issuer1  apr12  2.2    2           ...     Issuer1  apr12  3.3    4     ...     Issuer2  apr12  2.5    2     Issuer2  apr12  2.8    4     ...      Is there any way around this? (except for converting to a wide-format data, eg a separate spread variable for each maturity) Thank you


title: Standard deviation of several measurements with uncertainties

text: I have two 2 hours of GPS data with a sampling rate of 1 Hz (7200 measurements). The data are given in the form $(X, X_\sigma, Y, Y_\sigma, Z, Z_\sigma)$, where $N_\sigma$ is the measurement uncertainty. When I take the mean of all measurements (e.g. the average Z value of those two hours), what is its standard deviation? I can of course calculate the standard deviation from the Z values, but then I neglect the fact that there are known measurement uncertainties... Edit: The data is all from the same station, and all coordinates are remeasured every second. Due to satellite constellations etc., every measurement has a different uncertainity. The purpose of my analysis is to find the displacement due to the an external event, (i.e. an earthquake). I would like to take the mean for 7200 measurements (2h) before the earthquake and another mean for 2h after the earthquake, and then calculate the resulting difference (in height for example). In order to specifiy the standard deviation of this difference, I need to know the standard deviation of the two means.


title: Sampling frame of a voluntary survey

text: If one is conducting a survey in which all the respondents are volunteers, would there be a true sampling frame? I see three possibilities:   1. We could view this as a situation where there is no set list of subjects from which those that end up being included in the sample are chosen -- so therefore there is no sampling frame   2. The sampling frame could be every possible individual in the world that has access to a computer (it it's online) or exits in physical space and who might encounter the person administering the survey- say anyone in the world who might live in or might be visiting the city and might happen to walk by the corner where the person who is administering the survey might be standing. This would just mean that the vast majority of the sampling frame would have an extremely low inclusion probability   3. The sampling frame is all people who have a reasonable opportunity to take the survey- so if it were an internet survey, all of those who were asked to take the survey by email or had seen a page that had a link to the survey. Similarly if someone were administering a survey in say a mall, then the sampling frame would be all people who were shopping in that mall in a given time.


title: Specifying negative costs (benefits) in rpart's loss matrix

text: How do we specify negative costs in rpart? The documentation says the diagonals of the loss matrix should be zero. Is there an alternative to specify the benefits of correct classification (that is, the negative cost)?


title: Which statistical method should I use?

text: I wish to analyse the following : Independent Variable (IV): User of online financial reports Perceived usefulness and Perceived Quality (Sub to Relevance, Reliability, Understandability, Comparability and Timeliness) (all using 5-point Likert scale.) towards Internet Financial Reporting in Malaysia. Dependent Variable (DV): Intention to use Internet Financial Reporting information as a decision-making instrument (5-point Likert scale) Questions:   1. Which statistical method should I use to analyse my data? Previously I was thinking to use multiple regression and correlation to analyse it, but I just found that I'm not fulfilling the assumptions, because the DV should be measured on a continuous scale.   2. What are the minimum numbers of observations I need? I read before that it is at least 10 for each IV. Does that mean I only need 60 observations?  I'm using SPSS to analyse the data.


title: Entropy weighted Naive Bayes performs poorer than regular Naive Bayes?

text: I have a text classification problem, where there are many different classes, and the text to be classified is very short (about 1 sentence each):               this is sentence one, label1     this is sentence two, label2     ...      I am using the words in the the text as features with a naive bayes classifier. When computing: ![enter image description here](http://i.stack.imgur.com/zuPGi.png) the naive bayes classifier just tabulates counts. However, I wanted to weight the count for each word: ![enter image description here](http://i.stack.imgur.com/1Ym2z.png) Where the weight for a word/class pair $w_i, c_j$ is determined as follows:   1. Take the subset of the training data where $w_i$ occurs   2. Using this subset, compute the entropy on the distribution over class labels. The higher this entropy value the more "disordered" the class labels are for the given word $w_i$ (meaning this is a bad word with respect to classification), and vice versa. Normalize this by the maximum possible entropy.   3. Take the inverse of the above (now words with a low entropy, good words, have higher value)   4. Multiply the above by the fraction of times $c_j$ occurs in the subset mentioned in (1). The reason for this step, is because the quantity computed in (3), inverse normalized entropy, will be high if the word associates with a certain class, no matter what class. However, we don't want the weight to be high for a word/class pair if the word has a high inverse entropy which is high due to some other class we are not concerned with. # What went wrong? I felt that this weighting mechanism would improve results: the weight for a word/class pair is high when a word is informative of that class. However, after implementing the above, my results were worse than the normal naive bayes classifier with no weights! I am not sure why, any insight would be highly appreciated! (maybe someone with sharper math skills than me can point out a potential flaw?)


title: Using the rugarch package when sample size is small

text: I have a question about the `rugarch` package. My sample size is 43 and I have a problem to model a garch whose mean equation includes an exogenous model; otherwise my mean equation is linear regression, that is, $y_i=a+bx_i+e_i$ that $i=1, \ldots, 43$ and $e_i$ follow a $\text{garch}(1,1)$. But when I run it simultaneously R says `you need at least 100 points`. I don't know what to do.


title: How to quantify mis-specification bias and compare against smoothing bias for a non-parametric estimate of a randomly allocated continuous treatment?

text: Suppose that there is a data-generating process $$ y = \alpha + g(x) + \epsilon $$ which is to say that an outcome is some function of $x$. Suppose that $x$ is randomly assigned, so $\text{Cov}(x,\epsilon)=0$. But the functional form mappying $x$ to $y$ is unknown. Say that $y$ is continuous and the variance of $\epsilon$ is normal and iid. Think of $\epsilon$ not as measurement error, but rather as everything else that determines $y$ besides $x$. So, $g(x)$ represents the causal effect of $x$ on $y$. Say that you've got reason to believe that a straight line is a bad model for the data. Your instinct might be to use a gam:               library(mgcv)     model = gam(y~s(x))      But smoothing is model selection (via a ridge penalty), and model selection biases coefficients. This is no good, because you want unbiased causal effects. You can represent the various splines that make up a gam as $\beta_j$ for a given spline basis. Given a non-zero smoothing parameter $\lambda$, $E[\hat\beta_j] \neq \beta_j$. The bias-variance tradeoff for a smoothing spline is often presented as $$ PSS(\hat f, \lambda) = \displaystyle\sum_{i=1}^N\left( y_i-\hat f(x_i)\right)^2 - \lambda \int\left( \hat f''(x)\right)^2 $$ which is to say that as the smoothness increases, the bias increases. **This is the problem -- bias away from what?** The bias is away from the true values of $\beta_j$. I don't care about $\beta_j$, I care about $g(x)$, which is roughly represented by a set of penalized splines. So if I were to take the ridiculous position that my set of unpenalized splines (chosen by allocating them evenly across the quantiles of $x$) were the $true$ model, then indeed, my estimates of $g_{smooth}(x)$ would be "biased" away from my estimates of $g_{unpenalized}(x)$. So I want to be able to quantify something like $$ E[\hat g_\hat{\lambda}(x)] = g(x) + \text{smoothing bias} - \text{misspecification bias} $$ so that I can discriminate whether bias is minimized using a gam or minimized using a parametric approximation. Intuitively, if the true $g(x)$ looks roughly sigmoidal and I fit it with a straight line, then smoothing bias would be zero and mis-specification bias will be high. If I use a gam, then smoothing bias will be positive and mis- specification bias will be lower (though not zero, because the basis dimension is finite and $\lambda$ is estimated). Basically, I want to be able to argue that the choice of a smooth functional form is more or less biased than a linear regression coefficient for a completely exogenous, continuously-valued treatment. So how do I quantify mis- specification bias? One of the best forms of help here would be citations of papers to read that have thought about this. I'm not sure what terms to search for when scouring the literature. Note: I'll have control variables as well as $x$, so lowess is out. In addition, I'd like to look at heterogeneous treatment effects via interactions with $x$. So I'd like to see if its not impossible to use `gam` EDIT: no ideas, not even any comments?


title: Autocorrelation and partial autocorrelation interpretation

text: I am learning about ACF and PACF graphs. I am not sure I understand how to interpret the one I got for my data. I have searched google for some ACF and PACF examples, and I found some samples of different processes, however, the one I am getting doesn't look similar to any. Does this means there is no seasonality, trend and other processes? I have also created graphs for unstandardized residuals of my model, do not understand what those means? Does it somehow relate to white noise? ![ACF & PACF of y and residuals](http://i.stack.imgur.com/KK7vo.png)


title: Determine presence of "plateau" phase with lme4

text: I am trying to determine whether the data depicted in the following figure shows a plateau phase for higher scramblings. ![enter image description here](http://i.stack.imgur.com/QaBpZ.png) As I have often been encouraged to use linear models instead of anova - I have decided to tackle this issue with lme4. The output I get for my latex document from the texreg lme4 parser is as follows:                                    Model 1       (Intercept) : 2.03    [1.78; 2.28]*   COIsc-14 : -0.23    [-0.38; -0.08]*   COIsc-18 : -0.23    [-0.38; -0.08]*   COIsc-22 : -0.19    [-0.34; -0.04]*   COIsc-26 : -0.32    [-0.47; -0.17]*   COIsc-6 : 0.45    [0.30; 0.60]*   AIC : 1697.17   BIC : 1735.04   Log Likelihood : -840.59   Deviance : 1681.17   Num. obs. : 840   Num. groups: ID : 7   Variance: ID.(Intercept) : 0.08   Variance: Residual : 0.41   *0 outside the confidence interval                                    Model 1       My confidence interval is 0.95 . Can I make the point I am trying to make by the fact that the 14, 18, 22, and 26 values are all in each other's confidence intervals? Also, what information do I get from               Variance: ID.(Intercept) : 0.08           Variance: Residual : 0.41       ? Does that mean that 0.51 of my variance is explained by my categories, 0.08 by my IDs (which I define as random) and 0.41 by unidentified sources? 0.41 sounds like a lot, is this a good model then?


title: How to convert from VAR (vector-auto-regression) on differences to original series. How to simulate orignal series from VAR on differences

text: I have been using the VAR {vars} program to find a fit for the following bi- variate time series (subset):               bivariateTS <- structure(c(0.950415958293559, 0.96077848972081, 0.964348957109053,      0.967852884998915, 0.967773510751625, 0.970342843688257, 0.97613937178359,      0.980118627997436, 0.987059493773907, 0.99536830931504, 1.00622672085718,      1.01198013845981, 1.01866618122606, 1.02039666233818, 1.02400374633675,      1.02493877778841, 1.02865185614599, 1.03337978432753, 1.03244791539977,      1.03237166203917, 1.02853704067597, 1.0271704031946, 1.02588255626357,      1.02300209859012, 1.02397946952311, 1.02052043198189, 1.01593679111077,      1.01535467405129, 1.01368552421158, 1.00833115261092, 1.00495328099247,      1.00334161454411, 1.00029163432818, 0.999268774926758, 0.998174371988049,      0.999643729403106, 1.00255061235855, 1.0036328278581, 1.00343133578339,      1.00020064935273, 0.478881778679413, 0.404641679179503, 0.622804024030457,      0.474677494404522, 0.433851413612414, 0.544920292447026, 0.561093836992123,      0.563908924197768, 0.395583533713216, 0.273458857040352, 0.284537535240248,      0.277701982320208, 0.451280071722745, 0.401744885050984, 0.533760449322162,      0.488858723259857, 0.426799781286085, 0.503196832449693, 0.526637755320189,      0.599897302279671, 0.519986312138689, 0.460526054741527, 0.550734747489121,      0.487644621564283, 0.524566370547446, 0.700614666580712, 0.620626909353966,      0.661635856684872, 0.631149920530168, 0.668383815717266, 0.646705357249337,      0.685555413824799, 0.650175601510345, 0.679512593649786, 0.67242973935153,      0.599346651479619, 0.630852735978974, 0.657541680143728, 0.694629910826878,      0.654572164624051), .Dim = c(40L, 2L), .Dimnames = list(NULL, c("a", "b")))      I have used the "usual" techniques to show that both of the series are I(1), (integrated of order 1) AND that there is no co-integration. In addition I used the VARselect program to estimate the number of lags. To estimate the relationship between these variables the standard recipe (as far as I know) is to take first differences of the individual series and apply the VAR program to the new bivariate difference series. The VAR program seems to derive a reasonable set of coefficients. My Question:   **How do I simulate future paths for the original pair of (undifferenced) series?** * * * _The rest of this note is just how I have done this up until now..._ The method I have been using is to take a simulated path of the differences (note that in my case there are 2 simulated paths) say               > x.d     [1]   9  90 900      and then add back the starting value (1 in this case)               > c(1, 1+cumsum(x.d))     [1]    1   10  100 1000       A kind of re-integration of the series. In a real case I would use the final value(s) of `bivariateTS (1.0002006, 0.6545722 )` in place of the "1" in the cumsum above. My problem:   **The paths that I get make no sense.** Here is the function I use (an alteration of Professor Pfaff's predict.varest) to simulate var paths. My apologies if it is bad form to use alter someone else's code for other purposes.               # Modified VAR simulation function jigged up from the predict function of vars package     simVARpath<-function(object, n.ahead, mult = 1) {         K <- object$K             p <- object$p         obs <- object$obs             type <- object$type         data.all <- object$datamat             ynames <- colnames(object$y)         n.ahead <- as.integer(n.ahead)         Z <- object$datamat[, -c(1:K)]             B <- Bcoef(object)             if (type == "const") {                 Zdet <- matrix(rep(1, n.ahead), nrow = n.ahead, ncol = 1)                 colnames(Zdet) <- "const"             }             else if (type == "trend") {                 trdstart <- nrow(Z) + 1 + p                 Zdet <- matrix(seq(trdstart, length = n.ahead), nrow = n.ahead,                      ncol = 1)                 colnames(Zdet) <- "trend"             }             else if (type == "both") {                 trdstart <- nrow(Z) + 1 + p                 Zdet <- matrix(c(rep(1, n.ahead), seq(trdstart, length = n.ahead)),                      nrow = n.ahead, ncol = 2)                 colnames(Zdet) <- c("const", "trend")             }             else if (type == "none") {                 Zdet <- NULL             }             if (!is.null(eval(object$call$season))) {                 season <- eval(object$call$season)                 seas.names <- paste("sd", 1:(season - 1), sep = "")                 cycle <- tail(data.all[, seas.names], season)                 seasonal <- as.matrix(cycle, nrow = season, ncol = season -                      1)                 if (nrow(seasonal) >= n.ahead) {                     seasonal <- as.matrix(cycle[1:n.ahead, ], nrow = n.ahead,                          ncol = season - 1)                 }                 else {                     while (nrow(seasonal) < n.ahead) {                         seasonal <- rbind(seasonal, cycle)                     }                     seasonal <- seasonal[1:n.ahead, ]                 }                 rownames(seasonal) <- seq(nrow(data.all) + 1, length = n.ahead)                 if (!is.null(Zdet)) {                     Zdet <- as.matrix(cbind(Zdet, seasonal))                 }                 else {                     Zdet <- as.matrix(seasonal)                 }             }             if (!is.null(eval(object$call$exogen))) {                 if (is.null(dumvar)) {                     stop("\nNo matrix for dumvar supplied,                            but object varest contains   exogenous variables.\n")                 }                 if (!all(colnames(dumvar) %in% colnames(data.all))) {                     stop("\nColumn names of dumvar do not coincide with exogen.\n")                 }                 if (!identical(nrow(dumvar), n.ahead)) {                     stop("\nRow number of dumvar is unequal to n.ahead.\n")                 }                 if (!is.null(Zdet)) {                     Zdet <- as.matrix(cbind(Zdet, dumvar))                 }                 else {                     Zdet <- as.matrix(dumvar)                 }             }             Zy <- as.matrix(object$datamat[, 1:(K * (p + 1))])              cov<-summary(object)$covres         forecast <- matrix(NA, ncol = K, nrow = n.ahead)         lasty <- c(Zy[nrow(Zy), ])         for (i in 1:n.ahead) {             lasty <- lasty[1:(K * p)]             Z <- c(lasty, Zdet[i, ])             rv<-mult*matrix(rmvsnorm(1, dim = dim(cov)[1], Omega = cov), ncol = 1)             forecast[i, ] <- B %*% Z + rv             temp <- forecast[i, ]             lasty <- c(temp, lasty)         }         forecast     }


title: How to detect polarized user opinions (high and low star ratings)

text: If I have a star rating system where users can express their preference for a product or item, how can I detect statistically if the votes are highly "divided". Meaning, even if the average is 3 out of 5, for a given product, how can I detect if that is a 1-5 split versus a consensus 3, using just the data (no graphical methods)


title: How to correct for heteroskedasticity in fixed effects panel regression with correction for clustered standard error?

text: I am a student at RSM and I have a question regarding my regression analysis for my thesis as I have encountered issues I do not know how to deal with. I have performance data (dependent variable) of companies in 13 different countries, all of which are the ones I wish to regress against financial liberalization factors (independent variables). My time period run from 1980 to today. Each company has multiple observations over different time periods, after correcting for firm fixed effects as I wish to run a regression with _all_ my data for all countries in one. I have also included the clustering standard errors for the companies as the observations for one company are clustered. However, the issue now is that all my models appear to highly suffer from heteroskedasticity. After reading numerous econometric papers which are too technical for me to understand, I am clueless as to how I should proceed. Could you please advice me on how to fix this problem?


title: How to define CDF of bivariate normal distribution in rjags?

text: I am trying to define my own likelihood in rjags. I want to define the CDF of bivariate normal distribution under rjags. Does anyone know how to define the CDf of bivariate normal in rjags?


title: Is this a valid method to control the FWER?

text: I have a huge number, say $M$, of hypotheses that are potentially correlated. I have a dataset $D$ of random samples from an unknown distribution and I want to do test the hypotheses for significance using $D$, while controlling the FamilyWise Error Rate (FWER) to be, say, at most $\delta$. Some people have suggested me to use the following method, but I don't think it is statistically valid, and they can't seem to be able to argue about its being valid, so I thought about asking here. It is basically the holdout method but used in a multiple hypotheses testing setting. Take the dataset $D$, split it randomly into two datasets $D_{exp}$ and $D_{eval}$. Use $D_{exp}$ to filter out some of the $M$ hypotheses, for example testing them using critical value $\delta$ (i.e. with no multiple hypotheses correction). Suppose that only $k$ of them pass the test on $D_{exp}$. Test these $k$ on $D_{eval}$, using $\delta/k$ as critical value (or equivalently, using the Holm-Bonferroni method on those $k$ hypotheses), and accept those that pass this test. The idea is that $k\ll M$ and therefore this method has more power than applying the Bonferroni inequality directly (which is actually computationally infeasible in my case, as $M$ is potentially of the order of $10^{300}$.). Does this method actually control the FWER? If so why, and is there a reference for this, possibly with a proof? I believe it's of key importance that I'm in a multiple hypotheses testing setting, with potentially correlated hypotheses. Thanks in advance.