title: How does the confidence interval for the parameter p of a bernoulli trial vary with the value of p itself?

text: This is a very basic question (I'm currently studying undergrad level statistics), but I was hoping for some clarification regarding an assertion I read in a newspaper article earlier today. The author asserts that > evidence about the risks (of childbirth) has been hard to come by and > difficult to interpret. This is partly because the overall risks of maternal > and neonatal death are now very small (about five per 100,000 women die in > childbirth and four per 1,000 babies), so large numbers of mums are needed > to assess relative risks. Now, intuitively, this seems uncontroversial - if an event occurs rarely, then you'd expect more trials to be needed to achieve a good estimate of the likelihood of the event, when compared to one that occurs more frequently. However, I'm having trouble seeing how the mathematics bares this out. If we take the example in the article, we can model the probability of the death of a mother during childbirth as a Bernoulli trial with an estimate of the parameter $p$ given by $\hat{p}=\frac{5}{1000}$. From this, we can construct a 95% confidence interval for the true value of p: $$ p^{\pm} = \hat{p} \pm 1.96 \frac{p(1-p)}{\sqrt{n}} $$ However, if we take the limit as $p \to 0$ we get $$ \begin{aligned} \lim_{p \to 0} p^{\pm} = \hat{p} \pm 1.96 \frac {0(1-0)}{\sqrt{n}}\\ = \hat{p} \pm 1.96 \frac {0}{\sqrt{n}}\\ = \hat{p} \pm 0 \end{aligned} $$ Therefore, it appears that our estimate of $p$ in fact becomes _more_ accurate as $p$ gets smaller, regardless of our value of $n$. This seems pretty counterintuituve to me - am I going wrong somewhere, and if so, where? Many thanks, Tim


title: What statistics algorithm should be used to identify why something is increasing?

text: I work at a hospital and have been asked to use statistical algorithms to identify an increase in Census which is really the number of patient days in a hospital. It has trended up about 15 percent in the past three months and we have not been able to identify why using basic descriptive statistics. I have dozens of fields to choose from such as PatientAge, DiagnosisCode, Physician who admitted patient, InsuranceType, ProcedurePerformed etc. There is an integer field called units and the overall units has increased which means there are more patients coming to the hospital. I want to analyze all of these fields and identify which field or combination of fields is producing an increase in Units. I tried simple correlations but nothing is standing out. Could I use regression or decision trees to predict an increase in something? Thanks.


title: Why could mean centering change results

text: I used centering for my variables due to multicolleanirity and surprisingly the results (from before to after centering) changed for two interacted variables; one from significantly negative to significantly positive and the other from significantly positive to significantly negative. I read that results normally change for the main effect, but not for the interaction. > In my case the interaction results were changed > What could be the reason > and is this normal and what should I do? I am doing panel data analysis with fixed effects and here is my model > Y it = β0 + β1 X1 it +β2 X2 it + β3 X3 it + β5 M it + β5 (X1 it × 2010 i) + > β6(X2 it×2010i)+ β7 (X3 it × 2010 i)+ Year dummies + λi+ϵ it M is the interaction variable and X are the independent variables ... when I run the model first, the VIF was too high and then whether I center or use z score the results change Appreciate any feedback > thank you so much


title: EM algorithm Practice Problem

text: This is a practice problem for a midterm exam. The problem is an EM algorithm example. I am having trouble with part (f). I list parts (a)-(e) for completion and in case I made a mistake earlier. Let $X_1,\ldots,X_n$ be independent exponential random variables with rate $\theta$. Unfortunately, the actual $X$ values are not observed, and we only observe whether the $X$ values fall within certain intervals. Let $G_{1j} = \mathbb{1}\left\{X_j < 1\right\}$, $G_{2j} = \mathbb{1}\left\{1< X_j<2\right\}$, and $G_{3j} = \mathbb{1}\left\{X_j > 2\right\}$ for $j=1,\ldots,n$. THe observed data consist of $(G_{1j},G_{2j},G_{3j})$. (a) Give the observed data likelihood: $\begin{align*} L(\theta | G) &= \prod_{j=1}^n \text{Pr}\left\{X_j < 1\right\}^{G_{1j}}\text{Pr}\left\{1<X_j<2 \right\}^{G_{2j}}\text{Pr}\left\{X_j >2\right\}^{G_{3j}}\\ &= \prod_{j=1}^n \left(1-e^{-\theta}\right)^{G_{1j}}\left(e^{-\theta}-e^{-2\theta}\right)^{G_{2j}}\left(e^{-2\theta}\right)^{G_{3j}} \end{align*}$ (b) Give the complete data likelihood $\begin{align*} L(\theta | X,G) &= \prod_{j=1}^n \left(\theta e^{-\theta x_j}\right)^{G_{1j}}\left(\theta e^{-\theta x_j}\right)^{G_{2j}}\left(\theta e^{-\theta x_j}\right)^{G_{3j}} \end{align*}$ (c) Derive the predictive density of the latent variable $f(x_j|G,\theta)$ $\begin{align*} f(x_j|G,\theta) &= \dfrac{f_{X,G}(x_j, g)}{f_G(g)}\\ &= \dfrac{ \theta e^{-\theta x_j}\mathbb{1}\left\{x_j \in \text{region r s.t. } G_{rj}=1\right\}}{\left(1-e^{-\theta}\right)^{g_{1j}}\left(e^{-\theta}-e^{-2\theta}\right)^{g_{2j}}\left(e^{-2\theta}\right)^{g_{3j}}} \end{align*}$ (d) E-step. Give the function $Q(\theta,\theta^i)$ $\begin{align*} Q(\theta,\theta^i) &= \text{E}_{X|G,\theta^i}\left[ \log{f(\mathbf{x}|G,\theta)}\right]\\ &= n\log{\theta} - \theta\sum_{j=1}^n\text{E}\left[X_j|G,\theta^i\right] - N_1\log{(1-e^{-\theta})} - N_2\log{(e^{-\theta}-e^{-2\theta})} - N_3\log{e^{-2\theta}}\\ &= n\log{\theta} - \theta\sum_{j=1}^n\text{E}\left[X_j|G,\theta^i\right] - N_1\log{(1-e^{-\theta})} - N_2\log{(e^{-\theta}(1-e^{-\theta}))} + 2\theta N_3\\ &= n\log{\theta} - \theta\sum_{j=1}^n\text{E}\left[X_j|G,\theta^i\right] - N_1\log{(1-e^{-\theta})} + \theta N_2 -N_2\log{(1-e^{-\theta})} + 2\theta N_3 \end{align*}$ where $N_1=\sum_{j=1}^n g_{1j}, N_2=\sum_{j=1}^n g_{2j}, N_3=\sum_{j=1}^n g_{3j}$ (e) Give expressions for $\text{E}\left[X_j|G_{rj}=1,\theta^i\right]$ for $r=1,2,3$. I will list my results which I am pretty sure are right but the derivations would be a bit long for this already looong question: $\begin{align*} \text{E}\left[X_j|G_{1j}=1,\theta^i\right] &= \left(\dfrac{1}{1-e^{-\theta^i}}\right)\left(\dfrac{1}{\theta^i}-e^{-\theta^i}(1+1/\theta^i)\right)\\ \text{E}\left[X_j|G_{2j}=1,\theta^i\right] &= \left(\dfrac{1}{e^{-\theta^i}-e^{-2\theta^i}}\right)\left(e^{-\theta^i}(1+1/\theta^i)-e^{-2\theta^i}(2+1/\theta^i)\right)\\ \text{E}\left[X_j|G_{3j}=1,\theta^i\right] &= \left(\dfrac{1}{e^{-2\theta^i}}\right)\left(e^{-2\theta^i}(2+1/\theta^i)\right) \end{align*}$ This is the part that I am stuck on, and it might be because of an earlier mistake: (f) M-Step. Find the $\theta$ that maximizes $Q(\theta,\theta^i)$ From law of total expectation we have $\begin{align*} \text{E}\left[X_j|G,\theta^i\right] &= \left(\dfrac{1}{\theta^i}-e^{-\theta^i}(1+1/\theta^i)\right) + \left(e^{-\theta^i}(1+1/\theta^i)-e^{-2\theta^i}(2+1/\theta^i)\right) + \left(e^{-2\theta^i}(2+1/\theta^i)\right)\\ &= 1/\theta^i \end{align*}$ Therefor $\begin{align*} Q(\theta,\theta^i) &= n\log{\theta} - \theta\sum_{j=1}^n\text{E}\left[X_j|G,\theta^i\right] - N_1\log{(1-e^{-\theta})} + \theta N_2 -N_2\log{(1-e^{-\theta})} + 2\theta N_3\\ &= n\log{\theta} - \theta\dfrac{n}{\theta^i} - N_1\log{(1-e^{-\theta})} + \theta N_2 -N_2\log{(1-e^{-\theta})} + 2\theta N_3\\ \dfrac{\partial Q(\theta,\theta^i)}{\partial \theta} &= \dfrac{n}{\theta} - \dfrac{n}{\theta^i} - \dfrac{(N_1+N_2)e^{-\theta}}{1-e^{-\theta}} + N_2+2N_3 \end{align*}$ Next I should set this equal to zero and solve for $\theta$, but I have tried this for a very long time and I cannot seem to solve for $\theta$!


title: Why does repeated-measures ANOVA assume sphericity?

text: I mean by sphericity the assumption that the variance of all pairwise differences should be the same between groups. In particular, I don't understand why this should be the assumption and not that the variances of the group scores themselves be the same.


title: Is $F(E[Y_n]) \approx E[F(Y_n)]$ a reasonable approximation?

text: Studying the asymptotic distribution of order statistics I came across this approximation: $$F \left( E \left[ Y_n^{\left(n \right)} \right] \right) \approx E \left[ F \left( (Y_n^{\left( n \right)} \right) \right] $$ where $Y_n^{\left( n \right) }$ denotes the sequence of the nth order statistics and $F(.)$ is the distribution function of the underlying distribution. I understand that because of the Probability Integral Transform, $F \left( Y_n^{ \left(n \right)} \right)$ is the largest observation from a sample of size $n$ from a uniform distribution over $(0,1)$. Intuitively it makes sense to claim that the area to the left of the mean of $Y_n^{ \left( n \right) } $ is approximately equal to the mean of the area to the left of $Y_n^{ \left( n \right) } $. Is that really the case though? Also does this approximation hold for other order stastistics as well? All help is greatly appreciated, thank you.


title: Need help on Cramér-Chernoff method in concentration inequalities

text: I am reading up on the Cramér-Chernoff method in concentration inequalities. The idea is to use Markov's inequality and the monotonic transformation $\phi(t) = e^{\lambda t}$ where $\lambda \geq 0$. Consider a non-negative random variable $Z$ with finite expectation. Then Markov's inequality implies $$ (1) \qquad \Pr(Z \geq t) = \Pr(\phi(Z) \geq \phi(t)) \leq \frac{\mathbb{E}e^{\lambda Z}}{e^{\lambda t}} $$ In the Cramér- Chernoff method, we continue by letting $$ \qquad \psi_Z(\lambda) = \log \mathbb{E} e^{\lambda Z} \quad \text{ for all } \lambda \geq 0, $$ then standard texts claim that we have $$ (2) \qquad \Pr(Z \geq t) \leq \exp \left( \inf_{\lambda \geq 0} ( \psi_Z(\lambda) - \lambda t ) \right). $$ For example, see page 7/22 of the lecture note http://users.math.uni- potsdam.de/~blanchard/lectures/lect_2.pdf. My question is: how we can obtain the last inequality? I understand that since (1) is true for all $\lambda \geq 0$, we can have $$ \Pr(Z \geq t) \leq \inf_{\lambda \geq 0} \frac{\mathbb{E}e^{\lambda Z}}{e^{\lambda t}} = \inf_{\lambda \geq 0} \left( \mathbb{E} \exp (\lambda Z - \lambda t) \right), $$ but I am confused with how this can lead to (2). I'd appreciate for any help. Thank you in advance.


title: Calculating Perplexity

text: In the Coursera NLP course , Dan Jurafsky calculates the following perplexity: Operator(1 in 4) Sales(1 in 4) Technical Support(1 in 4) 30,000 names(1 in 120,000 each) He says the Perplexity is 53. The lecture can be found here: https://class.coursera.org/nlp/lecture/129 @ 8:40. I've been trying to figure out how it's calculated but I just can't get it. I've tried to switch to log's because of the 1/120000 but I still can't seem to figure it out.


title: Which distributions are parameterization invariant when based on the Jeffreys prior?

text: I understand that the Jeffreys prior provides a method for constructing a prior distribution over parameters for a given model (likelihood function) such that the prior distribution is "invariant under reparameterization." I understand this invariance to mean that the Jeffreys prior for a given set of parameters can be converted to a prior distribution over a second set of parameters (via the standard change of variables method for probability distributions), and the resulting prior will match the Jeffreys prior for the second set of parameters. Does a similar kind of invariance exist for:   1. the posterior distribution based on the Jeffreys prior? (i.e. does the posterior derived from a Jeffreys prior possess the same invariance properties as the Jeffreys prior?)   2. the prior predictive distribution based on the Jeffreys prior? (i.e. does the prior predictive distribution derived from a set of parameters and the corresponding Jeffreys prior match the prior predictive distribution derived from a second set of parameters and the corresponding Jeffreys prior?)   3. the posterior predictive distribution based on the Jeffreys prior? (similar to #2)   4. the MAP (maximum _a posteriori_ ) parameter estimate based on the Jeffreys prior? (i.e. does the distribution over the data space based on one set of MAP parameters match the distribution based on the MAP parameters for a different parameterization, where both parameter posterior distributions are based on the corresponding Jeffreys priors)


title: How to find 3D normal distribution function from a sample of 3D data?

text: I have a series of three dimensional data like so:                   1    2     3     4    5          1   0    0     0     0    0          2   0   0.3   0.5   0.3   0          3   0   0.5   1.0   0.5   0          4   0   0.3   0.5   0.3   0          5   0    0     0     0    0      Which is clearly a normal distribution. I want some way to give this as input, and to get as output the mean and standard deviation. Is there a program which can do this for me? Preferably an online solution, if not, then a free program for linux. However, if necessary I am willing to code it myself. Finally, as a last resort I might be able to get the program for free from my school. Note: I was considering placing this in the Mathematics site, but I thought it bordered more on data. Feel free to move it if necessary.


title: Minimum number of studies for meta-analysis?

text: Is there any reference or guideline regarding the minimum number of studies for conducting a meta-analysis?


title: What is the process for computing the power of a difference in mean test?

text: Suppose I have two samples with the following statistics and I want to determine if there is a difference in their means. First:                  Min. 1st Qu.  Median    Mean 3rd Qu.    Max.          0.0    39.0    39.0    38.6    40.0    48.0       Others:                  Min. 1st Qu.  Median    Mean 3rd Qu.    Max.         4.00   39.00   39.00   38.51   39.00   50.00       Performing a `t.test` in R I get: > t.test(first, other)                       Welch Two Sample t-test          data:  first and other     t = 1.3771, df = 8981.579, p-value = 0.1685     alternative hypothesis: true difference in means is not equal to 0     95 percent confidence interval:      -0.03304794  0.18912248     sample estimates:     mean of x mean of y 38.60095  38.52291       If our $\alpha$ was 0.05 we would not reject the null hypothesis. But suppose we assume that the means are not the same and the actual difference in the means is 0.078. How would I go about computing the power of this test given an $\alpha$ of 0.05? I believe I can approximate power by sampling from the distributions and comparing the means and counting the number of times the null hypothesis would be rejected with a given $\alpha$ but I want to confirm that my simulation is correct. Following @whuber's suggestion from this post I tried computing the power using the following:               first.mean = 38.60095     first.sd = 2.791901     first.count = 4413          other.mean = 38.51326     other.sd = 2.634525     other.count = 4735          z = replicate(10^3, {       first = rnorm(first.count, mean=first.mean, sd=first.sd)         other = rnorm(other.count, mean=other.mean, sd=other.sd)         t.test(first, other)$p.value     })     length(z[z < 0.05])/length(z) # mean(z < 0.05)      I get 0.30 (thanks for catching the typo @John). Is there a way to confirm this analytically? This is from Allen Downey's ThinkStats book chapter 7 exercise 7. I took a stab at simulating the power and now want to confirm my result analytically if possible. Thanks!


title: Inconsistency between RMSE and 95% CI Coverage

text: I am running simulations to compare different weighting methods to estimate the mean of y (with missing values). I use bias, RMSE and 95% CI Coverage as my performance metrics. However, looking at RMSE and coverage, I arrive at different conclusions. For example, for method A and B, bias(x10^3) are 75 and -134 respectively, the corresponding RMSEs(x10^3) are 158 and 159, while coverages are 85% and 34%. Given that the variance of method A is higher than method B, and the weighting model is non-linear, is it normal to have such RMSE-coverage inconsistency? I have been searching for references which point out to this inconsistency. But all I could find was this statement: "In selecting a best single model or model ensemble, the three types of assessments may have to be balanced. The model with the lowest root mean squared error (RMSE) may do worse on trend or have coverage that is too high or too low..." http://www.pophealthmetrics.com/content/pdf/1478-7954-10-1.pdf And if this is the case, can we think of Coverage as a more important performance metric?


title: Bayes decision theory: Classification error probability

text: In Bayesian decision theory: Given $\omega_1$ and $\omega_2$ as two classes for classification, $P\left( \omega_1 \right)$ and $P\left( \omega_2\right)$ their prior probabilities, $x$ the feature vector representing an unknown pattern, $P\left( \omega_1 | x\right)$ and $P\left( \omega_2 | x\right)$ as posteriori probabilities, $p\left( x |\omega_i\right)$ the likelihood function of $\omega_i$ with respect to $x$; and $R_i$ the region of feature space where decision is in favor of $\omega_i$ why is it that the minimization of the probability error $P_e$ , given by: \begin{align*} P_e &= P\left( x \in R_1, \omega_2 \right) + P\left( x \in R_2, \omega_1 \right)\\ &=P\left(\omega_2 \right) \int \limits_{R_1}p\left( x |\omega_2 \right) dx + P\left(\omega_1 \right) \int \limits_{R_2}p\left( x |\omega_1 \right)dx\\ &= \int \limits_{R_1}P\left( \omega_2 | x\right)p(x) dx + \int \limits_{R_2}P\left(\omega_1 | x \right)p(x)dx \end{align*} by choosing regions $R_1$ and $R_2$ of feature space so that: \begin{align} R_1 &: P\left(\omega_1 | x \right) > P\left( \omega_2 | x\right)\\ R_2 &: P\left(\omega_2 | x \right) > P\left( \omega_1 | x\right) \end{align} is said to be **not always the best** for minimizing $P_e$ ?


title: Graduate research methods: what statistical test should I use?

text: I hope someone her can help me because I'm so lost. I'm trying to figure out what statistical test I should use for my study: I want to look at how promotions, stadium age, and environmental sustainability impact attendance at sporting events. What statistical test should I use? If someone could help me I'd really appreciate it


title: How do I go about conducting model diagnostics on WLS?

text: I'm familiar with the diagnostics required for OLS, however I'm in new territory with a model I'm fitting to data in R, using Poisson regression with GLM. What are the standard methods in evaluating a WLS model?


title: Where can I find good publicly available data that I could use to teach z-scores to my college students?

text: I am sick of using the examples in the book. Is there an easy place to find data for which z-score/percentile/normal distribution stuff would be easy to see?


title: When should I *not* use R's nlm function for MLE?

text: I've run across a couple guides suggesting that I use R's nlm for maximum likelihood estimation. But none of them (including R's documentation) gives much theoretical guidance for when to use or not use the function. As far as I can tell, nlm is just doing gradient descent along the lines of Newton's method. Are there principles for when it's reasonable to use this approach? What alternatives are available? Also, are there limits on the size of the arrays, etc. one can pass to nlm?


title: Getting the equation from R's lm when using a product

text: Similar to this question: How to translate the results from lm() to an equation? in which the top voted answer said how to get the form of an equation from `lm(y ~ x)` and equivalently for `lm(z ~ y + x)` and other sums, I'm wondering: how can the equation be obtained from the form of `lm(z ~ y*x)`? I have:               > summary(lm(log(z) ~ x*y))          Call:     lm(formula = log(z) ~ x * y)          Residuals:           Min        1Q    Median        3Q       Max      -0.181142 -0.073755  0.000481  0.082088  0.200902           Coefficients:                  Estimate Std. Error  t value Pr(>|t|)         (Intercept)  -9.85368    0.09304 -105.906  < 2e-16 ***     x           -97.41166    6.28269  -15.505  < 2e-16 ***     y            -2.26398    0.14243  -15.895  < 2e-16 ***     x:y          91.69016    9.77390    9.381 6.95e-14 ***      It returns 4 coefficients -- `(Intercept)`, `x`, `y`, and `x:y` -- but I'm not sure how to put them together to get the final equation. Is it simply that `x:y` term multiplied by `x*y` plus the intercept? That is, in this case $\log z = 91.69 xy -9.853$?


title: Learning a univariate transform (kernel?) for novelty detection

text: I have 150 observations, 500 features, and I am interested in novelty detection (outlier detection): given a new observation (let's say 'patient') I want to know if it is different from the previous ones (let's call it 'control'). If I had a lot of data, I would probably be using statistical testing at the univariate parameter level, but, because of multiple testing issues, I end up exploring the tails of the control distribution to achieve significance, and I do not have enough data for non parametric tests for such small p-values. I am doing one class SVMs that alleviate this issue by learning a global decision strategy. The limitations of this approach are   1. it is very 'blackboxy'   2. it works poorly if the data is very 'anisotropic', i.e. the marginal distributions of the control are very dissimilar in different directions. A trick to work around problem 2 is to center an norm the univariate parameters (this is often called creating 'Z scores'). Ideally, one would like to whiten the data using the control covariance, but there is not enough data to compute it. The values fed in the OC-SVM then can be interpreted as a univariate test statistic (under a normal null distribution for the controls). In my case, I can see from the histograms that the control's distribution is heavy tailed. I would like to learn a univariate transform making it closer to a standard normal. By the way, I have no reference about such practices. I have learned them empirically, and from lab discussions. Any pointer would be welcome, even if they don't directly answer my question.


title: Monte Carlo estimation of convex hull overlap probability

text: This is a statistical version of my Math.SE  post. Given natural numbers $b$ and $r$, uniformly randomly choose $b+r$ points within a unit square. Call the $b$ points the _blue points_ and the $r$ points _red points_. Let $p(b,r)$ be the probability that $H_r$, the convex hull of red points, overlaps with $H_b$, that of the blue points. How do I efficiently estimate $p(b,r)$ via Monte Carlo (or perhaps better) methods? I can think of averaging the probability over a large number of randomly chosen test cases, but how can I sample the search space in a way that gives me error bounds?


title: Comparing two lists to measure for similarity

text: The problem I am trying to solve is finding the probability that two people are the same by cross-referencing the associates of the two people. For example, if person A is associated with the following people:   * Jeff   * Rick   * Jessica   * Mary Person B is associated with the following people:   * Ryan   * Mary   * Dennis   * Scott   * Jeff   * Sharon   * Rick   * Larry   * James So these two people have the following people in common:   * Mary   * Jeff   * Rick How would I go about figuring out the likelihood that Person A and Person B are the same based on the common relationship with the three people above? There are three factors I can see right now, but I don't know how to weigh any of them:   1. Ratio of common associates (doubled because seen from both sides) over the total number of associates   2. Ratio of common associates over the number of associates for Person A   3. Ratio of common associates over the number of associates for Person B I'm not a statistician, so I don't know if what I've presented is the correct way to solve the problem. Can anyone provide some guidance?


title: Gaps in time series and time series validity

text: After doing some reading on CrossValidated, I understood that we can use "imputation" techniques to fill in the gaps (if they are random). But I am not clear on following questions:   1. How many consecutive gaps may make data set invalid for forecasting?   2. How many total gaps in data set makes it as invalid. For example I have hourly data for week, which means 188 total points in my data set   1. Case 1: assume if we are missing 3 consecutive days of data is missing, can we still consider that data set as valid data set?   2. Case 2: assume overall 80 data points are missing out of 188 points, can we still consider data set as valid? I am using HoltWinters implementation in java for forecasting. Any help would be appreciated.


title: What kind of t-test?

text: I have 3 algorithms (A, B, C) which predict "author name" of each text article. I have tested all the algorithms on 3 different corpora (X, Y, Z) (different number of articles from 3 newspapers). I have got the following result. The result provides the number of "Correct" and "Incorrect" answers and the "Accuracy" (correct_answer/total) for each Algorithm on each Corpus:               Algorithm A =>    Corpus X: 80%,    Corpus Y: 93%,    Corpus Z: 94%     Algorithm B =>    Corpus X: 78%,    Corpus Y: 90%,    Corpus Z: 89%     Algorithm C =>    Corpus X: 90%,    Corpus Y: 83%,    Corpus Z: 97%      I have to compare the results by T-test and sign-test in order to define which one is the best algorithm. But I do not know which kind of T-test I have to use (independent, two-tailed, paired…). Also I have heard that I have to use one-way ANOVA instead of t-test.


title: GLM and categorical variable R - remove one category

text: I am currently running quasipoisson models with a continuous response variable and 13 covariates. I am using `glm()` and `summary()`. Most of my covariates are continuous but I also have the `YEAR` and `SEASON` as categorical variables. When I run the model it appears that the year 2008 is not significant while all the other years are. How can I remove only this year from my model? Do I have to remove the entire covariate?


title: Confusion related to derivation of the probability distribution

text: I was reading a paper where they showed that when $\mathbf{X}$ and $\mathbf{Z}$ are two multidimensional variables $$p(z_{nk}=1|\mathbf{X}, \mathbf{Z}_{\neg nk})\propto p(\mathbf{X}|\mathbf{Z})p(z_{nk}=1|\mathbf{Z}_{\neg nk}).$$ I didn't get how this expression was obtained. Suggestions?


title: Would it be incorrect to apply tests appropriate for normal distribution to data with the following distribution?

text: ![enter image description here](http://i.stack.imgur.com/yc6o9.png) Histogram showing the frequency for test scores (range 15-75) for a total of 91 students.


title: Convergence in distribution (central limit theorem)

text: > If $X_1, ... , X_n$ be iid exponential with mean $1/\lambda$. Let $S_n = X_1 > + ... + X_n$. > > a) Show that $S_n$ is $\Gamma(n, 1/\lambda)$. Each $X_i$ is $\Gamma(1, 1/\lambda)$ by the definition of an exponential distribution. We can look at the mgf of $S_n$. $$M(t) = E(e^{t \sum X_i}) = E(e^{tX_1})...E(e^{tX_n}) = (\frac{1}{1-t/\lambda})...(\frac{1}{1-t/\lambda}) = (\frac{1}{1-t/\lambda})^n,$$ which is the mgf of $\Gamma(n, 1/\lambda)$. So $S_n \sim \Gamma(n, 1/\lambda)$. > b) Show that $\frac{(S_n - n/\lambda)}{\sqrt{n}}$ converges in distribution > and identify the asymptotic distribution. By the Central Limit Theorem, we know that $\frac{(S_n - n\mu)}{\sqrt{n}(1/\lambda^2)} = \frac{(S_n/n - \mu)}{(1/\lambda^2)(1/\sqrt{n})}$ converges in distribution to $N(0,1)$. We see that $\frac{(S_n/n - 1/\lambda)}{(1/\lambda^2)(1/\sqrt{n})} = \frac{(S_n/\sqrt{n} - \sqrt{n}/\lambda)}{(1/\lambda^2)}(\frac{\sqrt{n}}{\sqrt{n}}) = \frac{S_n - n/\lambda}{(1/\lambda^2)(\sqrt{n})}$ So, $$\frac{S_n - n/\lambda}{(1/\lambda^2)(\sqrt{n})} \xrightarrow{D} N(0,1)$$ implies that $$ \frac{S_n - n/\lambda}{(\sqrt{n})} = (1/\lambda^2)\frac{S_n - n/\lambda}{(1/\lambda^2)(\sqrt{n})} \xrightarrow{D} N(0, 1/\lambda^4).$$ > c) Show that $\sqrt{n}(\sqrt{S_n/n} - \sqrt{1/\lambda})$ converges in > distribution and identify the asymptotic distribution. > > There is a theorem in our textbook ( _Introduction to Mathematical > Statistics_ by Hogg) which states, > > _Let $\{X_n\}$ be a sequence of random variables such that $\sqrt{n}(X_n - > \theta) \xrightarrow{D} N(0, \sigma^2)$. Suppose the function $g(x)$ is > differentiable at $\theta$ and $g'(\theta) \not= 0$. Then $\sqrt{n}(g(X_n) - > g(\theta)) \xrightarrow{D} N(0, \sigma^2(g'(\theta)^2)$_. From part b), we see that $$\frac{\sqrt{n}(S_n/n - 1/\lambda)}{(1/\lambda^2)} \rightarrow N(0,1).$$ So $$\sqrt{n}(S_n/n - 1/\lambda) \rightarrow N(0, 1/\lambda^4)$$. Now we apply the theorem. We have $g(x)=\sqrt{x}$. So $g'(x)=1/(2\sqrt{x})$, and $g'(\frac{1}{\lambda}) = \frac{\sqrt{\lambda}}{2}$. So by the theorem, $$\sqrt{n}(\sqrt{S_n/n} - \sqrt{1/\lambda}) \rightarrow N(0, (\frac{\sqrt{\lambda}}{2\lambda^4})).$$ Do you think my answers are correct?


title: Hypothesis Testing With Only Two Data Points

text: Let's say I have monthly data on cars on the road and number of car accidents. I want to know whether cars on the road is somehow related to the number of car accidents. Of course, I could conduct regression analysis or correlation analysis. However, when I have a set of 50 or so months with the number of cars on the road and number of accidents, how can I perform hypothesis testing to test whether a rise in cars on the road was related to an increase in automobile accidents. Tagged as R because I'll be using R though my question relates to 'how would I solve this basic problem"


title: Question about Hybrid Monte Carlo

text: In a Hamiltonian system, the Hamiltonian is always preserved, but in HMC algorithm, new state is accepted by probability $\min(1,\frac{\exp(-H_{new})}{\exp(-H)})$, I think increase or decrease in Hamiltonian is due to numerical error, why should we always prefer a state with low Hamiltonian while what we actually want is low potential energy?


title: How do I model the probability of two empirical distributions?

text: I have two distributions: A, and B. Each distribution is filled with the numbers 1.0-10.0. These distributions are NOT simple functions, like the gaussian, but are merely empirical counts. Essentially, I want to create a model for the probability that any given number is A. You can imagine that this is easy using histograms; we would create 1.0 sized bins, count the number of A in that bin Bin(A), count the number of B Bin(B) in that bin, and create a new bin for that range on a new histogram with the height value being the percentage Bin(A) / (Bin(A)+Bin(B)). My question is how to do this using continuous random distributions. Using either Python or R works fine. I feel as though there is something critical I am missing or failing to understand about this problem, because while it seems like a very trivial problem to me, I can find little information on how to solve it in either of those languages, both of which I am fairly experienced with


title: What test(s) can I use to compare 4 percentages for both pre- and post-treatment conditions?

text: Have 4 groups of results from 4 different treatments. Want to compare pre- and post treatment within groups and between groups


title: How to interpret Pearson's $r$, Spearman's $\rho$, Cramér's $\phi_{c}$ and Cramér's $V$ correctly?

text: Is it right that both Cramér's $\phi_{c}$ (and consequently Cramér's $V$ as well) are similar to the Pearson correlation coefficient, $r$, (and so, to Spearman's $\rho$) in its interpretation, as it is written here (http://en.wikipedia.org/wiki/Phi_coefficient)? Does a Cramér's $V$ value of 0.30 indicate a stronger relationship between two variables than a Pearson correlation coefficient of 0.20? Consequently: Can I compare the strength of the relationship between different variables? Let's say for example. The strength between income and consumer spending is stronger than the relationship between distance (in classes) and consumer spending or gender and consumer spending?


title: Using kknn regression on a time series

text: I'm trying to perform a knn regression on a sales time series using the kknn package. I can't see how you can restrict the neighbours to past sales only though. e.g.               Product, Price, Sale Date     A,       £100,  2012-01-04     A,       £100,  2012-01-05     A,       £200,  2012-01-06      If I want to predict for the second row, it should only consider the first row as a neighbour -- not the third. Is there a way to do this or am I going to have to write something custom? If it's going to have to be custom, what's the most efficient way of doing this in R? My dataset is 200,000 rows so I imagine it could be slow going using a loop. Any help you can offer is appreciated. Thanks


title: Open source implementation elastic net in C or C++

text: Can anyone provide or point me to a freely available implemention of Elastic Net in C or C++?


title: What are the dangers of violating the homoscedasticity assumption for linear regression?

text: As an example, consider the `ChickWeight` data set in R. The variance obviously grows over time, so if I use a simple linear regression like:               m <- lm(weight ~ Time*Diet, data=ChickWeight)      which aspects of the model will be questionable? Are the problems limited to extrapolating outside the Time range? How tolerant is linear regression to violation of this assumption (i.e., how heteroscedastic does it have to be to cause problems?).


title: Repeated categorical measures in a survey

text: Last year, we did a survey that asked people a question about their belief in climate change. The respondents could select 1 from a list of 5 categorical responses. In the fall, there was an event that we believe might affect respondents' climate change beliefs, so we're repeating the survey this year. I haven't done much in the way of longitudinal analysis: how can I compare the responses to see if the proportions have changed? Thanks!


title: Principal component analysis with Voxels (Vector Pixels) in Matlab

text: I have Vectors which have a certain position (Voxels) in an image. I would like to perform a pca to cluster out all Voxels which are correlating with each other. I have for example three Voxels in 1D:               syn_data_1 = [1;0;0;1;1];     syn_data_2 = [1;0;0;1;1];     syn_data_3 = [0;0;1;0;1];     syn_data = [syn_data_1, syn_data_2, syn_data_3]     %syn_data(:,1) is the Voxel in position 1 in 1D etc      Now Position 1 and 2 are strongly correlating, while 3 doesn't. It is possible to use `corr` to see which Voxels are correlating, but it is impossible to do it for all Voxels on a big data set in 3D. Is there a way to perform a pca on such data, so that I can see a clustering of Voxels which are similar? PS: Please don't be confused by the word Voxel, since in the end of the day I am just talking about Pixels which are having several properties represented by a vector. I am of course happy to provide further information if this can help to understand my question.


title: What theories should every statistician know?

text: I'm thinking of this from a very basic, minimal requirements perspective. What are the key theories an industry (not academic) statistician should know, understand and utilize on a regular basis? A big one that comes to mind is Law of large numbers. What are the most essential for applying statistical theory to data analysis?


title: Detecting increases in twitter list activity

text: I'm trying to find more information on how to measure and detect increases in activity in twitter list activity. I have access to all the tweets produced by a list of users. The idea is to try to predict or detect when there is some breaking news that everyone is talking about. A sensible thing is to plot number of tweets over time, and proportion of active users in a list (if 10 out of 30 users in a list tweet, proportion of active users = 0.33) The problem with this is that it is difficult to spot spikes in activity due to daily rhythm - people tend to sleep, and only tweet during the day... who knew? A quick and dirty example: (Cycle plot?) http://savedbythegoog.appspot.com/?id=6708dc9a2c90f9f937755542900a1b1957f20f3a Looking at a month of data is slightly better - but not really informative: http://savedbythegoog.appspot.com/?id=014b2246ba3d762b5201f0048812cd0675d5721a Those peaks could be annotated with the top hashtags to give some meaning to them. That's what I got so far. Are there any specific things i could read up on to get a better measure or make things more interpretable? Any other approaches or ideas more than welcome!


title: What's statistic should I use to determine if Time Series A (economic) is a leading indicator of Time Series B (Sales)

text: I'm trying to utilize economic and industry data to forecast sales of industrial companies. I have numerous economic and industry indices and as a first step I would like to indentify those indices that are good leading indicator for sales. I have already utilized graphing, but I was hoping there was a statistic that would also provide some insight. Some preliminary research has uncovered the Granger Causality statistic, but I was hoping to see if there were any other recommended alternatives or to confirm that Granger Causality is the current best practice. Any input or suggestions would be greatly appreciated.


title: R: Visualizing document clustering results

text: I have a k-means clustering result with 35 clusters, there are 5000 documents that each belong to one of the 35 cluster. I would like to visualize the results of the clustering algorithm on a scatter plot (or something similar) where each document is colored based on which cluster they belong to, and their distance on the visualization is proportional to their distance in similarity (i.e. the more similar they are, the closer they appear on the visualization). Ideally, it would also be nice to see the top 10 words that belong to the clusters. I am attaching my code for the clustering algorithm, it deals with data from a database.               myCorpus <- Corpus(VectorSource(userbios$bio))         docs <- userbios$twitter_id     # convert to lower case     myCorpus <- tm_map(myCorpus, tolower)     # remove punctuation     myCorpus <- tm_map(myCorpus, removePunctuation)     # remove numbers     myCorpus <- tm_map(myCorpus, removeNumbers)     # remove URLs     removeURL <- function(x) gsub("http[[:alnum:]]*", "", x)     myCorpus <- tm_map(myCorpus, removeURL)     # add one extra stop words:  "via"     myStopwords <- c(stopwords('english'), "twitter", "tweets", "tweet", "tweeting", "account")               # remove stopwords from corpus     myCorpus <- tm_map(myCorpus, removeWords, myStopwords)               myTdm <- TermDocumentMatrix(myCorpus, control = list(wordLengths=c(1,Inf), weighting=weightTfIdf))     # remove sparse terms     myTdm2 <- removeSparseTerms(myTdm, sparse=0.90)          m2 <- as.matrix(myTdm2)     #cluster terms     distMatrix <- dist(scale(m2))     fit <- hclust(distMatrix, method="ward")     # transpose the matrix to cluster documents (tweets)      m3 <- t(m2)          # k-means clustering      k <- 35     kmeansResult <- kmeans(m3, k)     #cluster centers     round(kmeansResult$centers, digits=3)         for (i in 1:k) {           cat(paste("cluster ", i, ": ", sep=""))           s <- sort(kmeansResult$centers[i,], decreasing=T)           cat(names(s)[1:15], "\n")           # print the tweets of every cluster + #            print(docs[which(kmeansResult$cluster==i)])     }


title: appropriate hypothesis test for non-negative quantity

text: This question is intentionally written in a very abstract way, because I am interested in a generic solution, or for pointers to different possible concrete cases in different fields: I am measuring a particular quantity in a sample of subjects taken from a given population, and I want to apply a null hypothesis test to check whether there is evidence for "an effect" in that population. Under most circumstances, I would apply a one-sample t-test to test whether the population _mean_ is different from zero. In this case, however, I know that the underlying quantity can only be zero or positive, but _never negative_. In my view this invalidates the null hypothesis of a population mean of zero, because as soon as there is any variation across subjects in the population that mean has to be non-zero anyway. Or to put it differently: A t-test applied in such a situation effectively tests whether the _variance_ in the population is non-zero, which however we usually assume to be the case anyway, because no person is exactly like the other. The situation is complicated by the fact that I can only obtain estimates of that underlying non-negative quantity, where the sampling distribution can be safely assumed to be normal, so that in contrast to the true value these estimates can and will often be below zero. Two questions:   * _Is this a known situation with a generic approach?_ I couldn't find anything in standard textbooks on hypothesis testing.   * _Are you aware of instances of this situation?_ It might help me get to the core of the problem if I can compare my concrete case to other concrete cases in different disciplines.


title: Two player dice game probability

text: $A$ and $B$ play a dice game where a player wins if their score is higher. $B$ wins if their score is equal. What is the probability of $A$ winning if both the players roll their dice $n$ times and choose their best scores (out of their $n$ throws)?


title: Correct use of terms (about sampling distributions)

text: I am taking random numbers from an uniform distribution. After I have taken some number of samples, I plot the probability density function. Now, I'm wondering how should I call the plot. Is it probability density function of uniform sampling distribution? Due to small number of samples it does not look like exactly uniform distribution, so I don't know what would be the correct term. Near-uniform? Edit:               x = rand(1000, 1);     [f, t] = ksdensity(x);     figure;     plot(t, f);


title: Scaling count data that contains zeros

text: Imagine that I have 40 observations in region A, and 50 observations in region B. The total size of region A in square kilometers is 50km, and the total size of region B is 40km. Within each region, each observation has a count of events that ranges from 0 (fairly frequent) to 20. What is the best way to hypothesis test for a difference in mean events per square kilometer between region A and region B? I've considered the following options: 1) Using a negative binomial regression with the area of each observation as a covariate. 2) Scaling the dependent variable (event counts) somehow. 3) Subsetting A and B via matching to achieve an equivalent total area on each side. Thanks in advance.


title: Help for planning a neuroscience experiment

text: I'm new here. I am planning a neuroscience experiment. I will be measuring brain signals from about twenty subjects. I will present the subjects four different kinds of stimuli. After all the data processing, I will have four power spectrums for each participant (one per stimulus type). That is, I will see the magnitude (power) of each frequency (something like this but calculated from the brain data). There is also a fifth stimulus type, which is rest (no stimulus; used as as control measure or "baseline" in neuroscience jargon). From each of the power spectrums, I will pick two frequency bands and calculate the average power in each band. The first band is 8-12 Hz ("alpha" in neuroscience jargon) and the other 16-22 Hz ("beta" in neuroscience jargon). I am expecting to observe four phenomena:   * For stimulus type 1, alpha is larger than beta. Both are larger than baseline.   * For stimulus type 2, beta is larger than alpha. Both are larger than baseline.   * For stimulus type 3, alpha is larger than baseline activity. Beta is equal to baseline.   * For stimulus type 4, beta is larger than baseline activity. Alpha is equal to baseline. How should I do my statistical tests? Is my experimental design sufficient? Note that the power of beta is always less than the power alpha in the rest state. I am asking this because our field appears to have some problems regarding the rigour of statistical analyses (see e.g. http://www.nature.com/neuro/journal/v14/n9/full/nn.2886.html).


title: How can I analyze my incoming email?

text: I would like to analyze the email I receive in my Gmail inbox in order to systematically come up with effective Gmail filters for the most common types of email. I am prepared to manually curate and classify a very large volume of emails (perhaps several 1000) to supply training data, if a machine learning strategy requires it. What I am hoping to get from the analysis: You may be familiar with email management techniques such as Inbox Zero. Generally, they rely on triage of emails. Because I get a very large amount of irrelevant email (mailing lists which are mostly uninteresting, coworkers who unnecessarily CC unimportant email to everyone, emails which the sender thought I would want to read but was mistaken) triage takes too much time. Therefore, I want to see if statistical analysis of the text, type of attachment, sender and recipients can be used to automate triage of the "worst offender" classes of email. Some examples of actual data I would like to get:   * Word clouds for emails I consider important and/or urgent vs. email that I will ignore.   * Bayesian rules such as "emails containing a hyperlink or attachment are more important than usual if coming from person X, but not more important if from person Y, and less important from person Z" (perhaps X is my boss who sends me important documents, while Z is a coworker who only ever uses attachments to share cat pictures)   * A rule that could filter emails such as announcement lists, by checking if it contains an announcement I would like to see (based on keywords and/or scheduled time of announcement) I know there is a wealth of know-how in the tech industry for statistical analysis of email. Google, for instance, already does almost what I want with their new inbox tabs feature. Unfortunately their implementation is too opaque and gives too little control to the user, so I would prefer to do my own. How can I do this analysis? Can I use any existing tools, or must I program my own?


title: Interpreting a logistic regression model with multiple predictors

text: I performed multivariate logistic regression with the dependent variable `Y` being death at a nursing home within a certain period of entry and got the following results (note if the variables starts in `A` it is a continuous value while those starting in `B` are categorical):               Call:     glm(Y ~ A1 + B2 + B3 + B4 + B5 + A6 + A7 + A8 + A9, data=mydata, family=binomial)     Deviance Residuals:          Min       1Q   Median       3Q      Max       -1.0728  -0.2167  -0.1588  -0.1193   3.7788            Coefficients:                  Estimate Std. Error z value Pr(>|z|)         (Intercept)  20.048631  6.036637   3.321 0.000896 ***     A1           0.051167   0.016942   3.020 0.002527 **      B2          -0.664940   0.304299  -2.185 0.028878 *       B3          -2.825281   0.633072  -4.463 8.09e-06 ***     B4          -2.547931   0.957784  -2.660 0.007809 **      B5          -2.862460   1.385118  -2.067 0.038774 *       A6          -0.129808   0.041286  -3.144 0.001666 **      A7           0.020016   0.009456   2.117 0.034276 *       A8          -0.707924   0.253396  -2.794 0.005210 **      A9           0.003453   0.001549   2.229 0.025837 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 485.10  on 2206  degrees of freedom     Residual deviance: 417.28  on 2197  degrees of freedom     AIC: 437.28          Number of Fisher Scoring iterations: 7           (Intercept)           A1           B2           B3           B4           B5           A6           A7           A8           A9      5.093426e+08 1.052499e+00 5.143045e-01 5.929197e-02 7.824340e-02 5.712806e-02 8.782641e-01 1.020218e+00 4.926657e-01 1.003459e+00                              2.5 %       97.5 %     (Intercept) 3.703525e+03 7.004944e+13     A1          1.018123e+00 1.088035e+00     B2          2.832698e-01 9.337710e-01     B3          1.714448e-02 2.050537e-01     B4          1.197238e-02 5.113460e-01     B5          3.782990e-03 8.627079e-01     A6          8.099945e-01 9.522876e-01     A7          1.001484e+00 1.039302e+00     A8          2.998207e-01 8.095488e-01     A9          1.000416e+00 1.006510e+00      As you can see, all of the variables are "significant" in that their p values are below the usual threshold of 0.05. However looking at the coefficients, I'm not quite sure what to make of these results. It seems that although these variables contribute to the model, looking at the odds ratios, they don't seem to really seem to have much predictive power. Of note, when I calculated the AUC, I got approximately 0.8. Can I say that this model is better at predicting against mortality (e.g. predicting that seniors will live past the prescribed period) compared to predicting for mortality?


title: Calculating TF-IDF in Matlab for Stop Words

text: I have been calculating tf-idf for removing stop words in my file. I am experiencing problems. First let us have a look at code I have written..               %MATLAB CODE%          a = load('myexcelfile');       wcf     =   cell2mat(a.raw(2:end-2,2:7)); % obtaining frequencies in matrix form       %raw is a cell present within my excel file % end-2 I have done because last 2 were NaN               words = a.raw(2:end-2,1); % Obtaining words in array format          %CALCULATING IDF HERE%          m=numel(words);     for i=1:size(wcf,2)     nz = nnz( wcf(:, i));          % if not zero, assign a weight:     if nz     I(i) = log10( m / nz );     end     end          %CALCULATING TF HERE%          tf_all= wcf(:,1:6);      %Applying Formula%               for j=1:size(wcf, 2)     Y(:, j) = tf_all(:, j)*I(j);     end      Please tell me that Am I doing correct? and if yes, then Matrix I get comprises of values where I am unable to screen out stop words i.e. is are am the about etc. Actually I get weight low for frequent words but which threshold or on what basis should I use? Please guide and correct me anywhere Thanks in advance


title: Beta-Mixing Time Series

text: There are plenty of resources on how to compute $\beta$-mixing coefficients for a time series and to check if a time series is $\beta$-mixing or not. However, I am struggling to actually find concrete examples of time series I can generate on my own that are defined to be $\beta$-mixing or not. It would be great if someone could provide simple examples of both. Even better would be to provide rules that separate the generation of $\beta$-mixing time series from those that are not. Thank you.


title: Choosing a regression model

text: How can one objectively (read "algorithmically") select an appropriate model for doing a simple linear least-squares regression with two variables? For example, say the data seem to show a quadratic trend, and a parabola is generated which fits the data quite well. How do we justify making this the regression? Or how do we eliminate the possibility of there existing a better model? What I'm really worried about is this: we could just keep adding polynomial terms until we had a perfect fit for the data set (an interpolation of the points), with no error whatsoever. But this would be useless as far as _predicting_ or extrapolating, because there would be no reason to think that the "model" was actually appropriate. So how does one balance the needs of accuracy and intuitive appeal? (Also, please alert me if this has been asked before, I assumed it would have been but did not find anything.)


title: Should the outcome variable be measured at least twice for a longitudinal study?

text: I am trying to find the association between BMI and onset age of a condition with linear regression model. I have multiple records of BMI measurement. But the outcome variable, onset age of condition, obviously, only once. It seems outcome variable should be measured at least twice for a longitudinal study. Is this correct? If yes, then how do I make use of this multiple BMI observation in my linear regression model?


title: Chow forecast test

text: I have a question about the chow forecast test. I don't understand how you can proof this test. You have: (y_1;y_2)=(X_1 0; X_2 I_(n2))(B;y)+(e_1;e_2) and H_0 is equal to y=0. The restricted is equal to: (y_1;y_2)=(X_1 0; X_2 I_(n2))(B;0)+(e_1;e_2) and the unrestricted is equal to: (y_1;y_2)=(X_1 0; X_2 I_(n2))(B;y)+(e_1;e_2). But now you test with e_r=(e_1;e_2) and e_1=(y_2-Bx_1). I don't understand why y_2 drops.


title: Hypothesis testing with Neyman–Pearson (finding cutoff quantile of Poisson dist)

text: I am currently trying to solve a problem that should be very easy, yet I am stumped, and help would be GREATLY appreciated!! _**The background of the question is this:_** Q: The number of typing errors in a newspaper can be modeled by a Poisson distribution with some mean $μ$. Each of the $31$ students in a journalism course are allocated at random a past edition of the The Guard newspaper and asked to find the number of typing errors it contains. The total number of errors found by the students is $120$, i.e. $\sum_{i=1}^{n}y_i=120$. It has been hypothesized by the instructor of the course that there are an average of $5$ typing errors in each edition of The Guard. The editor of the newspaper has objected, claiming that the mean is less than $5$ per edition. You want to test this with at $5\%$ significance level. _**Here is where I begin my work_** So I want to test: $H_0:\mu_o=5$ and $H_a:\mu_a<5$. According to the Neyman Pearson lemma, we can create the most powerful test by looking at the likelihood ratio (LR). I simplified the LR and got the following which I believe is correct, hopefully! $$LR= e^{-n\mu_a+n\mu_0}(\frac{\mu_a}{\mu_0})^{\sum_{i=1}^{n}y_i}$$ Now looking at the ratio, and the fact that $\mu_a<\mu_0$, we can see that _effectively_ our test statistic is $\sum_{i=1}^{n}y_i$, and as it increases, we get less comfortable rejecting $H_0$, since it makes our LR smaller. _Now this is where I am stuck:_ Now I know that to reject $H-0$, I need to find the cut off value where $\sum_{i=1}^{n}y_i \le k$ where $k$ is the cutoff value. So since $\sum_{i=1}^{n}y_i$ is our test statistic, we can say that we have $n=31$ independent Poisson distributions which leads our test statistic to have $n\mu_0$ Poisson distribution. I am unsure if what I just did is correct or not, but assuming it is, we next want to find $k$. Using $R$, we have our test statistic model $n\mu_0$ $=$ $31\cdot4=124$ which we can use in $R$. `qpois(0.05` ( = significance level) `,124)` ( = our new parameter for the mean under H_0) = 106 which is incorrect. Hopefully this clears up the confusion! Thanks so much everyone!!


title: Simultaneous imputation of multiple binary variables in R

text: I have a dataset with multiple correlated binary variables (0/1). Can anyone point me toward a solution for imputing completely random missing values based on information in the other variables? Below, I provide some code to create a simplified dataset with just 3 correlated binary variables.               set.seed(123)     # create correlated random binary (0/1) variables     x1 <- runif(100,0,1)          # N(0,1))     x2 <- x1 * runif(100,0,1)     # N(0,1))     x3 <- x2 * runif(100,0,1)+0.2 # N(0,1))     x1 <- round(x1)     x2 <- round(x2)     x3 <- round(x3)          #introduce random missing (MCAR)     x1[seq(1,100,7)]<-NA     x2[seq(2,100,7)]<-NA     x3[seq(3,100,7)]<-NA     # how can I impute missing values in this dataframe?     df <- as.data.frame(cbind(x1,x2,x3))     cor(df,use="pairwise.complete.obs")


title: Plot triple optimization boosting procedure

text: I am looking for a way to graph a triple optimization used in boosting. I did individual cross-validation on shrinkage, number of trees, and interaction depth and found the boundaries of interesting values. Now I am trying to figure out "the best boost". My dataset contains 1515 rows, 1010 are training and 505 are test data.               intd <- c(5,6,7,8,9) # interaction.depth parameters     nt <- seq(1:10)*500  # number of trees     shrink <- seq(from = 0.001, to = 0.010, by = 0.001) # shrinkage factor        opt <- cbind(rep(nt, each=50),rep(shrink,each=5,times=10),rep(intd,times=100))          m <- matrix(0,nrow=500,ncol=505) # build empty matrix to fill with predictions     for (i in 1:500)   {     m[i,] <- predict((gbm(branch2, data=df[train,],distribution="gaussian",n.trees=opt[i,1],interaction.depth=opt[i,3], shrinkage = opt[i,2])), newdata=df[-train,],n.trees=opt[i,1])      cat(i," ")             }      So the question is, what is a sensible way to graphically represent the results? I'm guessing calculating the MSE will be necessary but then given there are three parameters I'm not sure how to best visualise it               boost.err=with(df[-train,],apply( (t(opt)-v10)^2,2,mean))      `df` contains all my data and `v10` is my response variable


title: Picture of Robert Wedderburn

text: I am going to present generalized linear models (GLM) to a small group of people. I like to talk about the contributors of such remarkable work; however, it is hard to find a picture of Robert Wedderburn, perhaps because he passed away too early. Can anyone let me know where to find a picture of him, please? Thank you!


title: How to obtain means of factor from a 'svyglm' object in R?

text: I am using `svyglm` and `svyttest` to weight my outcome with propensity scores as per: Ridgeway & colleagues (2013) "Toolkit for Weighting and Analysis of Nonequivalent Groups: A tutorial for the twang package" So after:               glm1 <- svyglm(X ~ Y, design=design.ps)     summary(glm1)      or               svyttest(X ~ Y, design=design.ps)      (where X is my dependent variable and Y, my factor), I obtain my means using:               svyby(~X, ~ Y, svymean, design=design.ps, na.rm=TRUE)       I also need to run "doubly robust" regression to control for covariates that remain unbalanced after propensity score weighting using:               glm1 <- svyglm(X ~ Y + a + b + c, design=design.ps)     summary(glm1)        or               svyttest(X ~ Y + a + b + c, design=design.ps)      (where `a`, `b` and `c` are my covariates). My question is how to obtain adjusted means and standard errors relating to my new model. The use of:               svyby(X ~ Y + a + b + c, svymean, design=design.ps)      is extremely inefficient (and I am yet to obtain an estimate).


title: Time Series Analysis and Forecasting

text: I am looking at ways to forecast monthly time series data over a larger geographic region. I have time series weather data (e.g., temperature, precipitation) from multiple stations, and the stations are at certain proximity to each other. In the data set, temperature has slightly increasing trend while precipitation does not. I want to forecast monthly temperature and precipitation for next 10 years using prior year observations. I want to combine prior temperature relationship with precipitation to forecast precipitation, and vice versa, which I think basically involves multivariate time-series analysis. Using the mean forecasted/predicted values and the range (i.e., 95th lower and upper percentiles or Standard deviation), I want to introduce inter-annual variability in the weather patterns. While introducing inter-annual variability in the weather pattern in each stations, I want to preserve the existing spatial relation observed among stations. For example, if a station A has a temperature of 11 C, Station B should have similar temperature although not the same. Has anyone have any idea about the ways to handle this problem? I have been exploring approach such as ARIMA and ETS, but have not figure out quite well how to come up with better estimates. Any suggestions would be highly appreciated. Also, I have been exploring the use of AMELIA in this case. Does anyone have opinion about the use of AMELIA as a forecasting approach? I know AMELIA is a imputation model but I don't know quite well if it has forecasting capability.


title: Discrete time hazard models (cloglog) in R

text: The `survival` package in `R` appears to focus on continuous time survival models. I am interested in estimating a discrete time version of a proportional hazard model, the complementary log-log model. I have a fairly straightforward survival model, with simple right censoring. I know that one way to estimate this model is to create a data set that has a separate row for each observation for each period in which it isn't "dead." Then, a `glm` model with the `cloglog` link can be used. This approach seems very memory inefficient; indeed, it would likely produce a data set that is too large for the memory on my machine. A second approach would be to code up the MLE myself. That would be simple enough, but I am hoping that there is a package that has this survival model canned. It would just be easier for collaboration and to avoid coding errors to use a package. Does anyone know of such a package?


title: Exploring a scatter-plot matrix for many variables

text: I am analyzing a dataset with many parameters (say, 50-200) and I am interested in looking at relations between variables (e.g. in terms of 2-variable scatter plots or 2d histograms). However, for this number of parameters it seems unfeasible to draw a 200x200 array of plots (unless I print it and hang on a wall). On the other hand, doing only a correlation matrix does not give all information about 2-variable relations. Is there a way (library or workflow) to explore 2-variable relations for many variables? I am especially interested in showing results to others (perhaps after some data preprocessing). E.g. something interactive in JavaScript, were I can see scatter-plot matrix for selected fields from a correlation matrix. By scatter-plot matrix I mean something like that: ![enter image description here](http://i.stack.imgur.com/aYMev.png) (taken from pandasplotting blog; aviable in Python / Pandas, R, D3.js, etc).


title: How to test the significance of increase in sample interval range(s)?

text: Suppose we have two samples of a variable taken under different conditions: e.g. A1 without medical treatment and A2 after medical treatment. These are not necessarily normally distributed. Suppose A1 ranges from 0 to 10 and A2 ranges -2 to 15. Now how do I test for the significance of these range increases?   Any thoughts? Thanks in advance.


title: Find out which factors are more important when creating index

text: I am very new to statistics and econometrics, I am programmer, but still I need to do one task and get familiar with statistics / econometrics more in details. 1) I have some index that consists of some factors. In total I have more than 30 factors. 2) I have the result index value and values of each seperate factor. This information is displayed for each country (about 80 indexes together with factor values). 3) Each factor has some weight attached to it when computing the result index and there is linear formula of computing this index. For those factors I need to understand - which factors are more important a) for all countries together b) for one particular country Can you, please, point me:   * taking into account econometrics/statistics terminology, what could be meant by "factor is more important" when forming index? "The more important" is something that shows, that if we remove this factor, then index would change a lot. "Not so important" could be factors that can be removed or the can be unknown, but index would not change much. Does this make sense?   * which method(s) should I start with to find this out?  P.S. I will use SPSS (or maybe R) as I am familiar with those a bit. P.S.2. Pleaes, note that I am not familiar with statistics in details yet and I do not have any experience, but I know some basics and I am wishing to understand more in details if you could help me.


title: regression with indepdent variables that I know the relative weight of

text: I want to do a regression of a continuous dependent variable (y) (say number of crystals formed) using several independent variables (also continuous), each represent the temperature of each of the plates that are encountered in the crystal-making process (this is a made-up scenario to illustrate the problem I'm having) (say 4, lets call them x1, x2, x3, x4), the thing is, I know the relative "weight" of each IV (time spent on each plate). For example, it takes 100 minutes to form crystals, and this 100 minutes are spent staying on 4 different plates with different temperature. the amount of time spent on each plate is different, say (20, 30, 10, 20, 20). so in this case, the raw materials for the crystal stay on plate 1 (x1) for 20 %of the total time, and stay on plate 2 (x2) for 30 % of the total time. and the amount of time on each plate stays the same WITHIN each batch. but the amount of time on each plate is different BETWEEN batches. I have 10 batches. my question is, how do I do the regression of y (number of crystal formed) against the IV (x1-x4, the temperatures of each of the 4 plates used to make the crystal), taking into account of the time the raw materials stay on each plate, when I do separate regression for each batch individually (the amount of time spent on each plates stays constant)? and when I do one regression for all of the batches together (so there are different amounts of time for each plate for each batch)? many thanks!


title: Using the results of a non-metric multidimensional scaling to derive a measure of deviance

text: I'm working with a dataset containing the taxonomic compositions of floras all over the world and performed a non-metric multidimensional scaling (NMDS) based on the relative proportions of plant families within each flora. **My question:** Is it mathematically correct to use the (Euclidean) distance to the coordinate origin in order to quantify the deviance of a flora from a hypothetical average composition?


title: What do p-values for levels of a categorical variable represent in Poisson regression?

text: I have a Poisson model with varying densities:               set.seed(1)     df = data.frame(density = 1:5, events = rpois(2000, 1:5))      If I regress on this, I get that the intercept is approximately `log(3)`, which makes sense because 3 is the mean of 1:5.               glm(events ~ 1, df, family = poisson)  # returns 1.089      But now suppose I want to read back the coefficients of density:               glm(events ~ as.factor(density), df, family = poisson)      (For simplicity I've used `density` as both the ID of the field and its density.) I would expect the coefficient of `density[i]` to be `log(3-i)` because the intercept would still be 3. However, it doesn't seem like the intercept remains 3 - in this case, the intercept is set to `log(1)`. In playing around with this, it seems like glm sets the intercept to be the coefficient of the first factor. Now I'm starting to wonder what the p values in a glm regression indicate. Is the null hypothesis that `density[i]` is the same as the intercept (aka `density[1]`)? Or is it that `density[i] = mean(density)`?


title: Meta analysis of summary data or model of raw data?

text: I have data of the following form:               Sample  MeasureZZ  Time  Event     S1      M1         t1    e1     ....      I have calculated the influence of ZZ on survival by               coxph(Surv(t1,e1) ~ MeasureZZ)      I want to compare this result with external data from multiple sources. The data have the following format:               Sample  Dataset MeasureType  Measure  Time  Event     S1      D1      TypeA        M1A      T1    E1     S1      D1      TypeB        M1B      T1    E1     ....     S2      D2      TypeK        M2K      T2    E2     ....      Where MeasureTypes A..X are usually (but not always) well correlated between them (verified by Kendall's W) and should also be well correlated with ZZ above in the literature. Note that each dataset may use a subset of measure types A..X, but this subset is the same for all samples of that dataset. I expect different measure types across different datasets to give same survival information as ZZ (the underlying mechanism is similar). The most obvious solution, since I have raw data, would be to do:                coxph(Surv(time,event) ~ dataset + measuretype + measure)      The alternative would be to do survival models within each dataset for each measuretype, then do a meta-analysis of different measuretypes within each dataset and a new meta-analysis of all datasets (I can provide example code, but this is long). I have done both. I have also done a simpler (one-step) meta-analysis of measure vs the composite measuretype:dataset factor. All provide similar results, with few differences. Is there a reason to choose one approach over another? I am looking for the most robust method, although statistical power is obviously desirable. For what it's worth, differences between measuretypes are more predictable than differences between datasets. In fact, it could even be argued that pooling samples with the same measuretype from different datasets is not completely wrong. Note that I am using metafor (but can also use lme4 if necessary). Also, data have been converted to Z-scores (mean/sd) via scale() and trimmed at the 0.1 percentile for outliers. Thanks a lot!


title: Bayesian approach on games of chance with physical devices

text: Suppose we alter side 6 of a die to appear more than 1/6th of the time. We do not know the actual proportion of the time each side of the die will appear because all or some of the other 5 sides do not have the 1/6 proportion too. What do we need to be 99% confidence of each of the 6 sides of this altered die? How does the _a priori_ knowledge help us to determine how biased it is? What are the benefits of using Bayes instead of classic study?


title: Hyperplanes in Support vector machine

text: Why do we care about the Hyper-plane wT xi + b = -1 and wT xi+ b =1 if we are searching for an optimal hyper-plane wTxi + b =0 Also what is significance of the alpha's in the solution. could i get more clear ideas of Vapanik constant C used in soft margin SVM. Is there any applet so that i can view the effect of all the parameters of support vector machine.


title: Monte Carlo by time or by interval

text: Say I compute monte carlo output from input scenarios. Input are discrete time series. I choose time series as an example to make the problem more obvious - this could be also any curve. Computation for output is performed at each time: $f(input(t)) = output(t)$. On the time line, we can distinguish severale periods with their own characteristics. I want to compare MC outputs from two experiences. Both experiences are different in the way to compute input scenarios. For the first one, an input scenario is obtained from original time series adding to the value at each time an independant gaussian noise. For the second one, an input scenario is obtained from original time series adding to all time values in one zone the same randomly picked gaussian noise. For both experiences, we compute output for all scenarios. Then, we average each output curve zone by zone. If there are 4 zones for instances, we have, for each of the N outputs from MC, 4 points. We compute MC distributions for these outputs. I want to compare the distributions of outputs for both experiences, and I will use for that purpose a Kolmogorov Smirnoff test for example. In a theoretical viewpoint though, what should I expect ? That distributions are similar ? Please ask for precision if problem is not clear.


title: Jaccard Similarity - From Data Mining book - Home work problem.

text: Exercise 3.1.3 : Suppose we have a universal set U of n elements, and we choose two subsets S and T at random, each with m of the n elements. What is the expected value of the Jaccard similarity of S and T ? I am reading the book http://infolab.stanford.edu/~ullman/mmds/ch3.pdf


title: Converting odds ratios to Cohen's d for meta analysis

text: I'm conducting a meta-analysis largely on mean differences data, however I have several articles that report only an odds ratio, N, p value, and confidence interval, and I need to convert this information to Cohen's d. I am having particular trouble calculating the variance for d and the confidence intervals for d. Here's an example from the data: the odds ratio for depression among partner violent men is 3.37, p<.0001, CI=2.08-5.47, and N=128. I calculated a d of .2909, converted the p-value to a z-value and obtained a z value of 3.291, and then calculated the SE at .1603. From here, I'm stuck on calculating the variance and confidence intervals of d. Any help would be very much appreciated! I am using Wilson's effect size calculator (not any additional meta-analysis programs).


title: Neyman-Pearson Theorem Question

text: I have the following question for my homework: Suppose X~exp($\theta$). We want to test $H_0: \theta=1 vs. H_a:\theta=2$, based on a sample of size 2 - ${X_1,X_2}.$ a. Obtain the most powerful test (MPT) at a significance level $\alpha=.01$ What I am confused about is how to go about next steps. Since I have $H_0: \theta=1 vs. H_a:\theta=2$, should I just plug those in for $\lambda$ to get my likelihood function to be: $\frac{L_1}{L_0}=\frac{\Pi2e^{-2x}}{\Pi e^{-x}}$ Then I end up with: $\frac{L_1}{L_0}=\frac{2e^{-n2\Sigma x_i}}{e^{-n\Sigma x_i}}$ $\frac{L_1}{L_0}=\frac{2e^{-2\Sigma x_i}}{e^{-\Sigma x_i}}$ $\frac{L_1}{L_0}=2e^{-\Sigma x_i}$ I am not sure where to go from here or if I am doing this right. Could someone please let me know whenever they have the chance? Thanks!


title: Logistic Regression Vs Neural Networks

text: This is what I have understood about neural networks. There are $L$ layers from input to output (including both) and the equations as follows. $a(2) = \theta(1) x$ $a(3) = \theta(2)a(2)$ . . . $y = a(L) = \theta(n-1) a(n-1)$ So, we can write $a(L)$ as below $y = a(L) = \theta(1)\theta(2)\cdots\theta(n-1)x$ Let's suppose $\theta(1)\theta(2)\cdots\theta(n-1) = \Omega(L)$ Now we have $y = \Omega(L)x$ To solve the above equation, we can even have logistic regression techniques, which will be much easier. So, what is the point in doing it so complicated, when we are not including any higher order terms? Someone please explain


title: How to determine the sample size of a Latin Hypercube sampling?

text: I am designing an experiment with 5 variables. Some of the variables have 2 and others have 3 levels. The method that I am going to use is the Latin Hypercube, but I do not not what the sample size should be for a valid experiment?


title: Should I expect it to be a chicken or a penguin?

text: An alien is trying to classify a group of only chickens and penguins into, well, chickens and penguins by analyzing 3 independent boolean features A, B, C. If the animal (in reality) is a chicken, A holds for the animal 60% of the time, B holds 10% of the time, C holds 80% of the time. If the animal is actually a penguin, A holds for it 55% of the time, B holds for it 20% of the time, and C holds for it 90% of the time. The alien observes an animal (that must be either a chicken or a penguin) for which A and B hold, but C is false. Is the animal most likely a chicken? Is the answer just the greater of $0.6 * 1 + 0.1 * 1 + 0.8 * 0$ and $0.55 * 1 + 0.2 * 1 + 0.9 * 0$? Since the features are all independent?


title: How do I combine two predictors?

text: I am trying to classify a data set with 2 boolean values. I have two classifiers that may/may not be independent. The first one is 65% accurate, and the second one is 60% accurate. Can I combine the two somehow to get a better classifier? Or will the second one just weaken my prediction?


title: How to calculate confidence intervals for power function ($y=ax^b$) in R using NLS

text: I hope you can help me get some confidence in my confidence interval... I am trying to get the confidence interval for a particular (threshold) point on a predicted curve. I find the confidence interval (C.I.) looks funny, so I would like to know whether what I am doing is correct... Thank you for any advice you can give me!               #Here is the data to which I fit a nonlinear least squares regression (power function)     y <-c(0.5745373, 0.8255836, 0.7139635, 0.6318004, 0.8688738, 0.8341626, 0.9278573, 0.6548638,     0.6995722, 0.8410211, 1.0000000, 0.8512973, 0.7917790, 0.5722589, 0.6918069, 0.7117084,     0.5279689, 0.6121315, 0.5869292, 0.7332605, 0.5991816, 0.5470566, 0.5987166, 0.4440854,     0.3719892, 0.4394820, 0.4410862, 0.4966288, 0.4291826, 0.4221613, 0.3951223, 0.3595973,     0.4617549, 0.5106482, 0.5939970, 0.5340835, 0.6036920, 0.5181577, 0.3892170, 0.3667581)          x <- c(0.95149254, 0.57954545, 0.56763699, 0.57089552, 1.00000000, 0.66870629, 0.70833333,     0.53125000, 0.58776596, 0.78061224, 0.61551724, 0.63750000, 0.85000000, 0.52397260,     0.66870629, 0.50328947, 0.52192982, 0.40691489, 0.36016949, 0.37500000, 0.35915493,      0.53968254, 0.36334197, 0.23486842, 0.19615385, 0.35961538, 0.30039267, 0.26424870,     0.54838710, 0.23363874, 0.26936620, 0.09514925, 0.15324519, 0.32465278, 0.33909574,      0.31587838, 0.24401914, 0.28813559, 0.23181818, 0.29272959)          plot(x,y, xlim=c(0,1), ylim=c(0,1))     m1 <- nls(y~a*x^b, start=list(a=2,b=0.2))     summary(m1)     curve((0.86888*x^0.43042), add=TRUE)     points(0.2276312,0.4595148, col='red', pch=17) # this is the threshold value I                                                      # previously calculated      now I would like to know the Confidence Interval of the threshold               #lower CI     0.86888-(2*0.04578) # intercept minus 2* S.E.     0.43042-(2*0.06333) # exponent minus 2*S.E.     curve((0.77732*x^0.30376), add=TRUE) # plot the lower C.I.     0.77732*0.2276312^0.30376 # Now I enter the x coordinate of threshold point in the      # Confidence interval function, to get the y coordinate     points(0.2276312,0.4958526, col='red', pch="-", cex=1.5) # plotting the CI of threshold value          #Upper CI     0.86888+(2*0.04578) # intercept plus 2* S.E.     0.43042+(2*0.06333) # exponent plus 2*S.E.     curve((0.96044*x^0.55708), add=TRUE) # plot upper C.I.     0.96044*0.2276312^0.55708 # Now I enter the x coordinate of threshold point in the      # Confidence interval function, to get the y coordinate     points(0.2276312,0.4211113, col='red', pch="-", cex=1.5) # plot the CI of threshold value     lines(c(0.2276312,0.2276312),c(0.4211113,0.4958526), col="red")          # then check whether correct     confint(m1) # more or less, but slightly different, don't understand why?      ![nls confints](http://i.stack.imgur.com/PVhBA.png) To me, a C.I. that looks like this would make more sense...               curve((0.96044*x^0.30376), col="green",add=TRUE)     curve((0.77732*x^0.55708), col="green",add=TRUE)      ![NLS_confint2](http://i.stack.imgur.com/vTl6d.png) but that would be mixing the intercept of upper and exponent of lower C.I. and vice versa My question is then, are the black Confidence Intervals correct (and thus the thresholds red confidence interval?


title: How to analyse data from experiment?

text: I have a problem to choose a method for analysing data from my experiment. I have one within-subjects (all four stimuli presented to every subject four times) and four between-subjects factors and a categorical dependent variable (four possible responses). Also, I want to ask is it ok to use Chi square to test for gender differences in dependent variable since I have more males and females in groups because data are in long format and I have more than one row for every subject? I am using R and long format of data for the first time, so I am not sure how to test between-subjects factor effect with multiple rows for every subject. Thank you in advance! Alex ![data_example](http://i.stack.imgur.com/HvmYb.jpg)


title: Machine Learning on Percent/Continous Dependent Variable

text: I have a large dataset of 30,000 cases with 150 variables. I am looking for a few possible machine learning solutions/methods that I could try and use for cross validation. My dependent variable is a percent/continuous variable while all my independent variables are continuous or discrete categorical variables. I am only looking for one output which would be the prediction variable which provides a percent (continuous between 0 and 1). Currently I have run linear, logit and probit models, with probit showing the most promise. I believe Naïve Bayes is a simple one yet to try, but wondering what other types of machine learning models are out there that would give the desired output. **Update:** I am modeling Election Turnout based on past turnout within an individual precinct. One precinct may have 33.33% turnout while another might have 57.65% turnout. The possible outputs could be any continuous percent between 0 and 1. Any thoughts would be much appreciated! Thanks in advance!


title: How to reduce question set of a poorly designed questionnaire?

text: I have a very lengthy but poorly designed questionnaire, which has a very large number of questions. The questionnaire really needs to be completely rewritten, but prior to that being done, there is an organisational need to use it again in the short term. Given the length of the questionnaire, it would be helpful to reduce the number of questions. Some of the 'factors' are composed entirely of questions in a Likert score format - the average of these questions is taken to obtain the 'factor' score. It is likely that some of the questions only contribute noise to the score. Would it seem reasonable to undertake a factor analysis of only the questions lying under each 'factor' and look for a one factor solution, eliminating those questions which don't load to that factor?


title: How to verify if values are from Gamma distribution?

text: I have a set of numbers generated from Gamma distribution. How can I verify that these values were all randomly generated from Gamma or satisfy Gamma distribution?


title: AICc and K for categorical factors and interactions

text: I am new to multimodel inference. I am trying to create a model that has multiple categorical factors and possible interactions. For example say that my model is... > Y ~ X1 + factor(X2) + factor(X3) say each factor has two possible categories. R is only giving me AIC in my summary output. When calculating AICc would K be 3 or 5? also, if I have a potential interaction as in... > Y ~ X1 * X2 + X3 would K be 3 or 4? Also what do I do if there is an interaction? Do I just stop there or can I continue with analysis leaving the interaction term in the model.


title: GB2 model building

text: It's not clear how to build models when using the Stata `gb2lfit` module. Output includes two tables, parameter estimates and their logs. The significance tests can sometimes differ between the tables, with the paramter estimates being significant and the logs not. On which table ahould decisions about parameter significance be made? Since it's possible that models with fewer than 4 parameters will fit the data, one also needs some way of comparing models. If only one log or parameter is not significant, is it then safe to assume that it is equal to 0 (log) or 1 (parameter)? If I constrain parameters as shown in the documetation, then how are models compared? Should the log pseudo-likelihood be utilized with the lowest value representing the best model?


title: Can the bias introduced by lasso change the sign of a coefficient?

text: L1 penalized regression introduces a bias on your regression model but decreases the variance. When this bias is introduced, is it possible that the coefficient of $B$ changes sign? This would greatly affect the interpretation, where a regular regression might tell you an increase in $x$ increases $y$ by $B$, but maybe using lasso an increase in $x$ decreases $y$ by $B$


title: Weibull distributed data for correlation analysis

text: On visual inspection of P-P plots, my time data appears to fit a Weibull distribution better than a normal distribution. I will use the variable for correlations and as a predictor in regression. In SPSS, I could fit the data to a Weibull regression model. I have also read that you can Vincentize the data, but I am unfamiliar with this technique. However, for the correlation analysis: 1) Would I need to transform the data? 2) If so, how do I transform data from a Weibull to an approximate normal distribution?


title: Hypothesis testing for equivalence of two arrangements

text: I have two arrangements(i.e. permutations) of numbers. First one is the target/real arrangement. Second, is the observed arrangement. > e.g. > > Target := 1,2,3,4,5,6,7 > > Observed := 4,1,7,3,2,5,6 Any two elements in an arrangement is not equal. What kind of test should I use? p.s. I am not good in statistics. I am trying to evaluate a simulation model with real world data. Target arrangement is a sequence of real world events while Observed arrangement is the sequence of events which occurred in a simulation. My hypothesis is that these two are similar. --EDIT-- Can this be done using Sequence Alignment methods used in Bioinformatics? --EDIT-- Actually i have 30 samples (30 subjects participated). All target values and observed values for a particular sample are in the same range where range is [1,n] and n ~= 15.


title: Gibbs sampling for LDA -- does a small Dirichlet concentration parameter make a difference?

text: I'm using a Gibbs sampler for Latent Dirichlet allocation as described by Griffiths and Steyvers (http://www.ncbi.nlm.nih.gov/pmc/articles/PMC387300/). The sampling of a new topic $j$ for word $i$ is described by $P(z_i = j | \mathbf{z}_{-i}, \mathbf{w}) \propto \frac{n^{(w_i)}_{-i,j} + \beta}{n^{(\bullet)}_{-i,j} + W \beta} \frac{n^{(d_i)}_{-i,j} + \alpha}{n^{(d_i)}_{-i} + T \alpha}$ For the sake of the question, let's focus on $\frac{n^{(d_i)}_{-i,j} + \alpha}{n^{(d_i)}_{-i} + T \alpha}$, which describes how a word's topic is sampled in relation to the current topic counts. The term $n^{(d_i)}_{-i,j}$ is the number of words in the document (not including the word we're updating) assigned to topic $j$, while $n^{(d_i)}_{-i}$ is the number of total words in the document not including the current word, $T$ is the number of topics, and $\alpha$ is the concentration parameter for the Dirichlet prior. In many applications, it's reasonable for documents to only contain one or two topics, which leads to a very small concentration parameter $\alpha$, ($<0.01$). It seems that when $\alpha$ gets very small, the prior doesn't "exert much influence" in the sampling. The topics are mostly sampled in proportion to existing counts with only a very small $\alpha$ term added in. In my experience, this doesn't tend to steer the sampler toward distributions that favor only one or two topics (or at least not quickly at all), as the sampler just samples in proportion to the initial topic distribution. In the extreme case of $\alpha =0$, there should be exactly one topic per document, but even after thousands and thousands of iterations, my sampler does not seem to approach such a distribution. Is there any way for a small Dirichlet concentration parameter to "exert more influence" in Gibbs sampling?