title: Positive & negative shocks in a VAR model and impulse response function in R

text: I have two general questions about impulse response functions in R using the package **vars**. Take a look at this code:               library(vars)     data(Canada)          var.2c <- VAR(Canada, p = 2, type = "const")     plot(irf(var.2c, impulse = "e", response = c("prod", "rw", "U"), boot = T))      It produces the impulse responses to a shock in the variable "e". My Questions:   1. Is this a positive or a negative shock in e?   2. Assuming this a a positive shock, how can I change it into a negative shock? I think, the shocks in a VAR-model are symmetric so in order to obtain the negative shock I should change the sign of the impulse response function? Am I right so far? Thanks for any help!


title: Comparing the accuracy of nested distributions: is the larger model always better?

text: I have a random variable $X$ with density function $f(x|\theta)$. I could approximate $f(x)$ using one of two density functions: $g(x|\phi)$ or $s(x|\psi)$. Suppose that $g(x|\phi)$ is nested inside $s(x|\psi)$: for every value of the parameters $\phi$ I can find the corresponding set of parameters $\psi$ such that $g(x|\phi) = s(x|\psi)$, $\forall x$. Further, suppose that I have a infinite sample $x_1, x_2, \dots$ from $f(x|\theta)$, so I could estimate the two approximations exactly. My intuition is that $s(x|\psi)$ is more flexible than $g(x|\phi)$, so it should give me a better approximation. Hence my question is: is my intuition right? And if it is right, can we prove that "larger" distribution $s(x|\psi)$ is always "better" (for example it has lower Mean Integrated Squared Error MISE)?


title: How to draw a conclusion about an hypothesis if one study obtains p < .05 and a second study obtains p > .05?

text: I have a general question about comparing two studies: > Suppose study $A$ reports a $p$-value less than $0.05$. Also suppose study > $B$ reports a $p$-value greater than $0.05$. Assume that both studies have > the same $H_0$ and $H_a$. Both studies use different statistical techniques. Which study is "correct"? Would you have to compare the statistical techniques used? If study $A$ was published before study $B$ would there be a "publication bias" in that most people would accept study $A$ as correct?


title: Forward and backward algorithm for HMM

text: I am looking for example based explanation for HMM's. How is that the HMM has P(O|HMM model) will have 2T* N^T (T is the length of observation, N is number of hidden states.) And how this is solved by forward and backward algorithm. I need some example. \- Thanks in Advance


title: Is a Biplot suitable ? Domestic Violence VS Non-Violence data

text: I have a small data set that consists solely of counts. I have several variables. For example State, for each state (Mexican States), for example, Tamauplipas I have the counts of domestic violence and the counts of non- violence. For Age I have the same but for each age group... I would eventually have a table, for example, like this: ![Variables](http://i.stack.imgur.com/oFvm7.jpg) I naturally have many more columns. What these people want to know (plot) is what levels of what variable are associated with Non-Violent and the same for Violent. I was told in work to do a biplot but have been trying and dont get nice results. I'm not sure this is the best idea. I tried with a regular PCA biplot... What do you guys think? any comments would be appreciated. Greetings! J.


title: Holt-Winters and Abnormal termination in LNSRCH

text: I try to fit data with Holt-Winters function in R. Nevertheless, i am getting the following message:               ts1<-ts(data$nb_decl,frequency=53)     hw1<-HoltWinters(ts1)          Warning message:     In HoltWinters(ts1) :       optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH      'data' counts 313 lines and if I just delete or change the last value (or add 314th value), the code works... Am I doing something wrong? or what is the problem with my data? see bellow the data used:               data         week year nb_decl     1     00 2006       0     2     01 2006       0     3     02 2006       0     4     03 2006       1     5     04 2006       0     6     05 2006       1     7     06 2006       0     8     07 2006       0     9     08 2006       0     10    09 2006       0     11    10 2006       1     12    11 2006       0     13    12 2006       2     14    13 2006       1     15    14 2006       1     16    15 2006       0     17    16 2006       2     18    17 2006       0     19    18 2006       1     20    19 2006       0     21    20 2006       0     22    21 2006       0     23    22 2006       0     24    23 2006       1     25    24 2006       0     26    25 2006       1     27    26 2006       1     28    27 2006       0     29    28 2006       1     30    29 2006       0     31    30 2006       0     32    31 2006       0     33    32 2006       0     34    33 2006       0     35    34 2006       1     36    35 2006       0     37    36 2006       1     38    37 2006       0     39    38 2006       0     40    39 2006       0     41    40 2006       1     42    41 2006       1     43    42 2006       0     44    43 2006       0     45    44 2006       0     46    45 2006       1     47    46 2006       3     48    47 2006       2     49    48 2006       4     50    49 2006       2     51    50 2006       1     52    51 2006       1     53    52 2006       0     54    01 2007       0     55    02 2007       1     56    03 2007       1     57    04 2007       1     58    05 2007       0     59    06 2007       2     60    07 2007       0     61    08 2007       1     62    09 2007       1     63    10 2007       1     64    11 2007       1     65    12 2007       1     66    13 2007       1     67    14 2007       1     68    15 2007       1     69    16 2007       1     70    17 2007       0     71    18 2007       0     72    19 2007       1     73    20 2007       0     74    21 2007       0     75    22 2007       3     76    23 2007       0     77    24 2007       0     78    25 2007       1     79    26 2007       0     80    27 2007       2     81    28 2007       0     82    29 2007       1     83    30 2007       0     84    31 2007       0     85    32 2007       1     86    33 2007       0     87    34 2007       2     88    35 2007       1     89    36 2007       1     90    37 2007       1     91    38 2007       1     92    39 2007       2     93    40 2007       0     94    41 2007       3     95    42 2007       0     96    43 2007       0     97    44 2007       0     98    45 2007       3     99    46 2007       0     100   47 2007       0     101   48 2007       0     102   49 2007       0     103   50 2007       1     104   51 2007       1     105   52 2007       0     106   53 2007       0     107   00 2008       1     108   01 2008       9     109   02 2008       0     110   03 2008       0     111   04 2008       1     112   05 2008       0     113   06 2008       0     114   07 2008       2     115   08 2008       0     116   09 2008       2     117   10 2008       2     118   11 2008       1     119   12 2008       0     120   13 2008       0     121   14 2008       3     122   15 2008       1     123   16 2008       0     124   17 2008       1     125   18 2008       2     126   19 2008       1     127   20 2008       1     128   21 2008       1     129   22 2008       1     130   23 2008       2     131   24 2008       1     132   25 2008       1     133   26 2008       0     134   27 2008       0     135   28 2008       3     136   29 2008       0     137   30 2008       4     138   31 2008       1     139   32 2008       0     140   33 2008       1     141   34 2008       1     142   35 2008       2     143   36 2008       0     144   37 2008       1     145   38 2008       0     146   39 2008       0     147   40 2008       1     148   41 2008       0     149   42 2008       2     150   43 2008       1     151   44 2008       0     152   45 2008       1     153   46 2008       2     154   47 2008       0     155   48 2008       3     156   49 2008       3     157   50 2008       2     158   51 2008       0     159   52 2008       0     160   00 2009       2     161   01 2009       3     162   02 2009       3     163   03 2009       2     164   04 2009       3     165   05 2009       1     166   06 2009       1     167   07 2009       1     168   08 2009       1     169   09 2009       2     170   10 2009       3     171   11 2009       1     172   12 2009       1     173   13 2009       2     174   14 2009       1     175   15 2009       2     176   16 2009       4     177   17 2009       0     178   18 2009       0     179   19 2009       3     180   20 2009       0     181   21 2009       2     182   22 2009       0     183   23 2009       0     184   24 2009       1     185   25 2009       1     186   26 2009       1     187   27 2009       1     188   28 2009       3     189   29 2009       0     190   30 2009       3     191   31 2009       4     192   32 2009       0     193   33 2009       1     194   34 2009       3     195   35 2009       0     196   36 2009       2     197   37 2009       1     198   38 2009       1     199   39 2009       1     200   40 2009       2     201   41 2009       0     202   42 2009       0     203   43 2009       1     204   44 2009       3     205   45 2009       2     206   46 2009       2     207   47 2009       1     208   48 2009       0     209   49 2009       2     210   50 2009       0     211   51 2009       0     212   52 2009       0     213   00 2010       1     214   01 2010       2     215   02 2010       3     216   03 2010       2     217   04 2010       1     218   05 2010       2     219   06 2010       0     220   07 2010       2     221   08 2010       0     222   09 2010       1     223   10 2010       1     224   11 2010       1     225   12 2010       1     226   13 2010       1     227   14 2010       0     228   15 2010       1     229   16 2010       0     230   17 2010       1     231   18 2010       1     232   19 2010       1     233   20 2010       2     234   21 2010       0     235   22 2010       2     236   23 2010       4     237   24 2010       3     238   25 2010       4     239   26 2010       1     240   27 2010       1     241   28 2010       1     242   29 2010       2     243   30 2010       2     244   31 2010       2     245   32 2010       0     246   33 2010       0     247   34 2010       2     248   35 2010       0     249   36 2010       1     250   37 2010       1     251   38 2010       1     252   39 2010       1     253   40 2010       2     254   41 2010       4     255   42 2010       2     256   43 2010       1     257   44 2010       3     258   45 2010       1     259   46 2010       0     260   47 2010       6     261   48 2010       1     262   49 2010       1     263   50 2010       1     264   51 2010       2     265   52 2010      11     266    0 2011       0     267    1 2011       3     268    2 2011       2     269    3 2011       0     270    4 2011       3     271    5 2011       0     272    6 2011       4     273    7 2011       2     274    8 2011       1     275    9 2011       1     276   10 2011       1     277   11 2011       2     278   12 2011       2     279   13 2011       3     280   14 2011       1     281   15 2011       5     282   16 2011       3     283   17 2011       3     284   18 2011       1     285   19 2011       2     286   20 2011       1     287   21 2011       3     288   22 2011       0     289   23 2011       2     290   24 2011       1     291   25 2011       2     292   26 2011       1     293   27 2011       2     294   28 2011       2     295   29 2011       2     296   30 2011       1     297   31 2011       0     298   32 2011       0     299   33 2011       2     300   34 2011       0     301   35 2011       2     302   36 2011       1     303   37 2011       1     304   38 2011       1     305   39 2011       3     306   40 2011       3     307   41 2011       2     308   42 2011       2     309   43 2011       1     310   44 2011       2     311   45 2011       1     312   46 2011       1     313   47 2011       0      thank you for your help!


title: Maximum number of classes for RandomForest multiclass estimation

text: I have researched the internet|literature a lot on multiclass prediction to find out what is a realistic limit for the number of classes that can successfully be used for estimation when using a RandomForest method. The literature body on text mining sometimes comes up with really large numbers of classes (>1000), while most other "classical" cases described have a class count less than 6-8. Most of them describe handmade algorithms specifically designed for the particular problem, though, while I am interested in the performance of standard RF implementations (in R, for example). I have even started to analyse simulated data to find out more about it, but the problem is to generate data that simulates a lot of multiple classes yet has meaningful and realistic predictors. I know that the results depend largely on the number of observations in every class and the balance between class outcomes. For my data, I can safely assume that there will be enough observations per class, so that I can balance the data accordingly. So I am curious whether people have applied standard RandomForest implementations to multiclass problems with a class count >>10. Note that I am not talking about separating the estimation into multiple one-vs-all problems. Does anybody here have some real-life experience with that kind of data? Thanks in advance and best regards!


title: Validating a Varimax implementation?

text: I'm writing an implementation of factor analysis and I'm having trouble convincing myself (or even better, proving) that my varimax implementation is correct. What's the best way to prove that a Varimax implementation is working properly? I can imagine two different ways; does anyone know where I can find the following:   1. _Known good data sets._ A few "known good" data sets would be a great way to prove that my implementation is working right.   2. _Reference implementation._ A reference implementation (in any language is fine, but ideally Python or Java) would be useful because then I could just generate "known good" data sets myself. I have access to SPSS, but I can't figure out how to get SPSS to just give me a Varimax rotation of a certain data set, and I'd like to isolate the Varimax transformation from other computations. All of this is more complicated if the Varimax mapping is not one-to-one. (That is, if each unique input to Varimax has only one correct output.) If Varimax is one-to-one, then if the output of my implementation differs from the output specified by one of the methods above, then I know it's wrong. If Varimax is _not_ one-to-one, then I can conclude that my implementation is not the same, but I can't necessarily conclude that it's wrong. Does anyone know if Varimax has only one correct output for each unique, valid input? Thanks for your help, guys!


title: Multi-level beta regression package for R?

text: I want to run multi-level beta regression with two levels (years nested within geographic pixels). I know about betareg, but as far as I know it doesn't have a multi-level option. Thanks.


title: Is there a binomial regression model that captures data with fat tails?

text: Specifically, are there any binomial regression models that use a kernel with heavier tails and higher kurtosis than the standard kernels (logistic/probit/cloglog)? As a function of the linear predictor $\textbf{x}'\mathbf{\hat{\beta}}$, the logistic distribution   * Underestimates the probability of my data being in the tails of the distribution   * Underestimates the kurtosis, or clustering of data, in the middle of the distribution: This can be seen from a diagnostic plot of my fit: ![enter image description here](http://i.stack.imgur.com/ar6OW.png)   * The red line is the logistic CDF, representing a perfect fit   * The black line represents the fitted probabilities from my dataset (calculated by binning observations into 0.1 intervals of $\textbf{x}'\mathbf{\hat{\beta}}$, where $\mathbf{\hat{\beta}}$ is obtained from my fit)   * The grey bars in the background represent number of observations on which the true probabilities are based upon   * The grey areas are where the tail 10% of the data lie (5% each side). Ideally, any solution would use R. ## Edit Why am I talking about CDFs? Our GLM equation is: $$\mathbb{P}(Y = 1) = \mathbb{E}[Y] = g^{-1}(\textbf{x}'\mathbf{\beta})$$ Where $g$ is the link function. Further, if $g^{-1}$ is a valid probability distribution (i.e. monotonically increasing from 0 to 1, indeed the case with probit, logit, cloglog), then consider a latent (not directly observed) continuous random variable $Y^{*}$ whose distribution (CDF) is given by $g^{-1}$. Then by definition $$\mathbb{P}(Y^{*} \leq \textbf{x}'\mathbf{\beta}) = g^{-1}(\textbf{x}'\mathbf{\beta})$$ Equating the two equations above, we see the probability of $Y=1$ is exactly equal to the CDF of $Y^{*}$ $$\mathbb{P}(Y = 1) = \mathbb{P}(Y^{*} \leq \textbf{x}'\mathbf{\beta})$$ Hence I talk interchangeably about the expected response $\mathbb{E}[Y]$ and CDF of $Y^{*}$ over linear-predictor ($\textbf{x}'\mathbf{\hat{\beta}}$) space.


title: How to compute the steady state gain from the transfer function of a discrete time system?

text: Having the Transfer Function of a discrete system as such: $$H(z) = \frac{0.8}{z(z-0.8)}$$ I am asked to find the Steady State Gain of the system. I have the solution and it simply states:   Steady State Gain:   $$H(1) = \frac{0.8}{1-0.8} = 4$$ I have no idea why this is the case, and the book, the slides, the solutions, my (limited) notes and the internet are seemingly empty on the subject. Is the Steady State Gain of a system always the outcome of the Transfer Function applied to 1? That just sounds ridiculous, especially since I'm not finding any references to it online.   I was chased out of mathoverflow with this question, those guys really hate homework... Then again, who doesn't.


title: How to perform coefficient test in linear model with robust covariance matrix?

text: I estimate a linear model and then calculate a (White) heteroscedasticity robust variance covariance matrix. This allows me to do the following (where `hr.cov` is the robust covariance matrix):               > coeftest(model, hr.cov)          t test of coefficients:                        Estimate Std. Error t value  Pr(>|t|)         (Intercept) 157.146959  14.760426 10.6465  8.99e-15 ***     x             0.134499   0.033050  4.0696 0.0001575 ***     d            20.921722  16.043356  1.3041 0.1978438         x:d           0.068120   0.047416  1.4366 0.1566992         ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1       I would like to do a one sided test whether d<0 (as opposed to different from 0 in both directions which is what is shown above). I would like to be able to type something like `test(model,d<0,hr.cov)`. How should I go about this?


title: How to calculate the correlation between a variable and a matrix of variables

text: I'm trying to improve a factory quality control. I have some variables from the melting process (something like ten control variables) that changes trough time (a matrix of the values of those control per minute), and in the end I have a quality score for the final product (one single variable). I have one of this for each production batch. I want to know if you guys can help me saying how can I look for correlations with the quality score inside that matrix. I know that I can look to each control variable alone, but those variables interfere with each other. So it is necessary to look at the sistem as a one. Thank you.


title: Optimal Portfolio Selection using a single index model

text: 1. Stock | 2. (Return of stock i - Rf) / Beta(i) | 3. SUM [(Return of stock i - Rf) * Beta(i) / Var(error of i)] | 4. Beta(i)^2 / Var(error i)^2 | 5. SUM [Beta(i)^2 / Var(error i)^2]     * y 0.0395 13.7170 13.7170 342.2115 347.2115 C1     * x 0.0080 4.3307 18.0477 538.4952 885.7067 C2     * z 0.0067 1.9733 20.0210 294.4627 1180.1694 C3 Assume Rf = 0.02 and that the variance of the returns of the market is 0.0023. **If short sales are not allowed, what is the composition of the optimal portfolio?** I calculated C1 = 0.017541, C2 = 0.013667, C3 = 0.012397 So I compute that only the first stock has a excess return to beta ratio greater than C, which means the Z/Sum(Z) = 1? What is the return of the portfolio and standard deviation then? This doesn't make any sense. How could the optimal portfolio only contain stock Y? No other information at all is given.


title: What are common statistical sins?

text: I'm a grad student in psychology, and as I pursue more and more independent studies in statistics, I am increasingly amazed by the inadequacy of my formal training. Both personal and second hand experience suggests that the paucity of statistical rigor in undergraduate and graduate training is rather ubiquitous within psychology. As such, I thought it would be useful for independent learners like myself to create a list of "Statistical Sins", tabulating statistical practices taught to grad students as standard practice that are in fact either superseded by superior (more powerful, or flexible, or robust, etc) modern methods or shown to be frankly invalid. Anticipating that other fields might also experience a similar state of affairs, I propose a community wiki where we can collect a list of Statistical Sins across disciplines. Please submit one "sin" per answer.


title: Clustering microblogs

text: I have a microblog dataset with about 100k tweets and I would like to cluster them effectively using some less advanced algorithm. Is there anywhere I can find like a single-pass K-Means or similar algorithms. Thanks. I have tried using K-Means and it takes about 5 minutes to cluster 100k tweets, using standard pre-processing (tokenization, stemming). I have also used Online LDA, but I can't retrieve tweets from the topics.


title: Mean and Covariance of Office Hours

text: I have a homework assignment with the following question. I've attempted to solve the problems and my answers are below, but I feel like they are incorrect, and have no idea how to check. I would be very grateful for any help you can provide:               Posted office hours notwithstanding, students make calls at the offices of their      professors and have formed the impression that professors are more likely to be away     from their offices on Friday than any other working day. A review of calls,      1/5 of which are on Fridays, indicates that for 16% of Friday calls, the professor      is away from his or her office, while this occurs for only 12% of calls on every      other working day.  Define two random variables as follows:          X = 1 if call is made on a Friday; otherwise, X = 0          Y = 1 if professor is away from his or her office; otherwise, Y = 0               a.  Find the mean value of X, namely, E(X).          b.  Find the mean value of Y, namely, E(Y).          c.  Find (to four decimal places) and interpret the covariance between X and Y.        * * * My Answers: `A` : 0.2(1) = **0.2** `B` : (0.12)0.2 + (0.12)0.2 +(0.12)0.2 +(0.12)0.2 +(0.16)0.2 = **0.128** `C` : E(xy) – E(x)*E(y) :


title: If X and Y are correlated, but Y and Z are independent, is X and Z always independent?

text: If random variables X and Y are correlated, but Y and Z are iid, is X and Z always independent? I am able to prove they are independent for specific types of structure on X and Y (for example if $X=Y+U$ where $U$ is iid), but is this true in general. Not sure why I am having such a difficult time proving the general case?


title: Gradient and Hessian of a likelihood function where y is defined implicitly

text: $\newcommand{\implicit}{\mathrm{implicit}}$ I'm attempting to optimize model parameters $\theta$ by maximizing the likelihood function $$ f(y) = \ln \Bigl(\frac {n!}{k!(n-k)!}y^k(1-y)^{n-k}\Bigl) $$ where $y$ must be calculated iteratively, as it is defined implicitly by $g(\theta)=0$, as follows: $$ g(\theta) = -1 + \frac{a\cdot \theta_2}{\ln \bigl(\frac{-y}{y - 1} \bigl) - \theta_1} + \frac{b\cdot \theta_4}{\ln \bigl(\frac{-y}{y - 1} \bigl) - \theta_3} $$ ($n$, $k$, $a$, and $b$ are constants). As far as I can tell, using partial implicit differentiation, the gradient (first partial derivative with respect to $\theta_i$) would be defined by $$ \frac{\partial \operatorname{f}}{\partial \operatorname{\theta_i}} = \frac{\frac{\partial \operatorname{f}}{\partial \operatorname{y}}\cdot \frac{-\partial \operatorname{g}}{\partial \operatorname{\theta_i}}}{\frac{\partial \operatorname{g}}{\partial \operatorname{y}}} $$ This seems to match numerical approximations computed by R, but when I try to calculate the Hessian by taking the partial derivative of the gradient with respect to $\theta_i$, I get strange results that are nowhere close to numerical approximations. This is the formula I came up as my attempt at calculating the Hessian (second-order mixed partial derivative of $f$ with respect to $\theta_i$): $$ \frac{\partial}{\partial \operatorname{\theta_i}} \Bigl(\frac{\partial \operatorname{f}}{\partial \operatorname{\theta_j}}\Bigl) = \frac{\frac{\partial}{\partial \operatorname{\theta_i}} \Bigl(\frac{\partial \operatorname{f}}{\partial \operatorname{y}}\Bigl) \cdot \frac{-\partial \operatorname{g}}{\partial \operatorname{\theta_i}}}{\frac{\partial \operatorname{g}}{\partial \operatorname{y}}}+ \frac{\frac{\partial \operatorname{f}}{\partial \operatorname{y}} \cdot \frac{\partial}{\partial \operatorname{\theta_i}} \Bigl(\frac{-\partial \operatorname{g}}{\partial \operatorname{\theta_i}}\Bigl)}{\frac{\partial \operatorname{g}}{\partial \operatorname{y}}}+ \frac{\frac{\partial \operatorname{f}}{\partial \operatorname{y}} \cdot \frac{-\partial \operatorname{g}}{\partial \operatorname{\theta_i}}\cdot\frac{\partial}{\partial \operatorname{\theta_i}}\Bigl(\frac{\partial \operatorname{g}}{\partial \operatorname{y}}\Bigl)}{\Bigl(\frac{-\partial \operatorname{g}}{\partial \operatorname{y}}\Bigl)^2} $$ Considering that the values I get from these calculations don't match (not even close) the automatically generated numerical approximations returned by R, I suspect that I've done something wrong here. Can anyone spot an error with the Hessian formula? Is there a property of higher-order mixed partial implicit differentiation that requires a different approach? Thanks! UPDATE: Alecos Papadopoulos posted a solution that eliminates the need for iteration by solving directly for $g(\theta)$, thus providing an exact value for $y$. This works perfectly for this problem, as the calculation of the gradient and Hessian does not require implicit differentiation in this case! For proof of concept (and in case I come across any $g(\theta)$ functions that can't be solved directly), I'm still interested in figuring out what went wrong with my attempt at the Hessian. If anyone has any insight into a general solution for the Hessian using implicit differentiation, it is certainly welcome.


title: What is the difference in Bayesian estimate and maximum likelihood estimate?

text: Please explain to me the difference in Bayesian estimate and Maximum likelihood estimate?


title: What are the general strategies in creating a Probabilistic Graphical Model?

text: While there is lot of theory and probability in the background to understand, I wanted to know if there are any resources/quick pointers on what to consider while modeling a problem using Bayesian Networks. Specifically, I had the following questions: \- what kind of priors to use? \- what do we learn from the data? and how to do prediction? Is there a blog/document that explains these from the ground-up with simple use-cases? (eg: modeling user-click behavior, weather forecasting, diagnosis etc). I am looking for examples which are a little detailed beyond toy- examples (a bit close to real life situations) used in courses. Any help is appreciated.


title: How can you define below average, average and above average according to percentile rank?

text: A particular school evaluation lists abilities by percentile rank against the abilities of a population, and each ranking/ability is then given a group ranking of "Below Average", "Average", or "Above Average". Oddly, the "Average" group ranges from 16% to 84%, so a 16% in reading comprehension and an 84% in reading speed both receive an "average" for each ability, which strips the labels of almost all meaning. Ignoring the obvious observation that _all_ points above/below the average are "above/below average", how might one define, say, "poor, average, and excellent" or even "poor, fair, average, good and excellent" according to their percentile ranks, as measured against the population of students considered? This is not merely an academic question, but impacts how we advocate for our child's future schooling. Thank you.


title: How to built a software recognition algorithm

text: I'm trying to design an algorithm to determine which softwares are installed on a computer. I have several millions of files each corresponding to the list of softwares which are on a computer. I've built several regexps which allow me to recognize few key softwares. I've computed joint probabilities in order to know better which software are almost always installed on a same computer and conditional probabilities. With these I'd like to be able to determine which softwares are installed on a computer having just to detect few key softwares. How should I do ? What rule can I derive from joint et conditional probabilities ?


title: What to do with coefficients of log-transformed predictors in meta regression of effect sizes?

text: I am preparing data for a meta-analysis regression of effect sizes, in which I will study the determinants of effect size for a specific variable across several studies. I have some studies in my dataset that have used the logarithm of my predictor variable of interest in their regression models. How should I deal with these studies where this independent variable is log- transformed, so that they are comparable to all the other effect sizes (in the other studies) for that same variable that are not log-transformed?


title: Repeated Measures of Dependent Randomized Conditions

text: I'm designing an experiment to evaluate how learnable a software tool is. The basic setup is to have $n$ subjects perform $m$ tasks in a random order, and then test whether the subjects performed significantly better on the second set of the tasks than the first:               | Subject 1 | Subject 2 | ...      -------------------------------     |   Task 1  |   Task 3  |     |   Task 4  |   Task 6  |     |   Task 2  |   Task 4  |     -------------------------------     |   Task 3  |   Task 2  |     |   Task 5  |   Task 1  |     |   Task 6  |   Task 5  |      The tasks are not expected to take the same amount of time, but the subject can "time-out" on a task and be assigned the maximum allowed time. The maximum time is the same for each task. The null hypothesis is "subjects perform equally well on the first set of tests as the second set of tasks." Or more precisely: subjects complete the first and second set of tasks in the same amount of time. Based on this, I have a few questions:   * Is this this a reasonable study design for testing the hypothesis?   * Which statistical test(s) should be used to test the null hypothesis?    * How does the randomized assignment of the tasks affect the statistical test?   * If we didn't group the tasks into two sets, is the an appropriate test for determining if task order (# of tasks the user has already seen) significantly affects performance? In some ways, a paired t-test seems appropriate, however I'm not sure whether the random assignment of tasks to groups is permissible.


title: Is the book wrong or me?

text: it's my first question on this forum. I just purchased Probability and Statistics - Basic book from amazon to refresh and learn what I've never learned. I like the book it explains the subject well. However, I got a different result on one of the examples. **Example:** 2 fair dice are rolled. What is the probability of getting a sum less than 7 or a sum equal to 10? The books claims that P(A) = 15/36 and P(B) = 1/36. And base on that P(A or B) = 4/9. I agree on P(A) but I believe that P(B) should be equal to 3/36, including {(4,6), (5,5), (6,4)} rolls. Therefore, my final result is P(A or B) = 3/6. Am I missing something or the book got it wrong?


title: Cluster quality measures

text: Does Matlab provide any facility for evaluating clustering methods? (cluster compactness and cluster separation. ....) Or is there any toolbox for it?


title: Bernoulli or Binary RBM

text: For the classical RBM: $P(\tilde{h}|\tilde{x})=\sigma\left(\tilde{b}+W\tilde{x}\right) a$ and $P(\tilde{x}|\tilde{h})=\sigma\left(\tilde{c}+W^{T}\tilde{h}\right)$ for for hidden layer $\tilde{h}$ and visible layer $\tilde{x}$ , and weight matrix $W$ and bias vectors $\tilde{b}$ and $\tilde{c}$ Is the correct (/preferred) name for this:   * a Binary-Binary RBM,   * a Bernoulli-Bernoulli RBM,   * a Binary RBM,   * a Bernoulli RBM I see Binary RBM most often, in papers -- if it is specified at all, which it usually isn't, this type of RBM is the default. Bernoulli-Bernoulli RBM makes the most sense to me, as the elements in the visible and in the hidden layers are assumed to be Bernoulli distributed. Which means the take Binary values. Bernoulli-Bernoulli also works better if we have Gaussian-Bernoulli RBMs also being talked about, as this speaks of the distrobutions of each layer. However I read papers which say Gaussian-Binary RBM, which to me is not right as Binary is not a distribution, but Gaussian is. What is preferred?


title: What machine learning algorithms are used for internet advertising?

text: I.e. based on stats of site visiting, clicks, classifying visitors based on what they are interested in, time, day of week etc. Any frameworks/resources you can recommend? For example: we have a lot of advertisements with different prices. Task is to build a set of advertisements in a block to achieve the best CTR (in other words we need to show what people click)


title: Variance and multiple levels in clmm

text: I'm trying to fit a multi-level model where people are nested in groups. V1 is my outcome variable. G2 is my highest level group and G1 is a group nested within G2               m1 <- clmm(V1 ~ 1 + (1 | G2), data = df, Hess = T)     m2 <- clmm(V1 ~ 1 + (1 | G1) + (1 | G2), data = df, Hess = T)      Based on AIC values, m1 is the best fit, so I'll go from there.               m3 <- clmm(V1 ~ 1 + (1 | G1:G2), data = df, Hess = T)     m4 <- clmm(V1 ~ 1 + (1 | G2/G1), data = df, Hess = T)      Is the syntax for m3 nesting G2 in G1 or G1 in G2? What about for m4? Lastly, is the variance produced by summary(m1), summary(m2) etc. the variance that remains in the model? If I do a cross level interaction between G1:G2, how should I be interpreting the variance in the RE effects portion of the results?                            Var   Std.Dev     G1:G2    9.017e-11 9.496e-06     G1       1.683e-01 4.103e-01      I'm sorry for the fundamental-ness of this question and the lack of data.


title: Why does a kernel density estimate of a normally distributed vector have a non-smooth second derivative?

text: I had some normally distributed data:               mu <- 3     sigma <- 5     x <- rnorm(1e5, mu, sigma)      I took a kernel density estimate with a fairly high bandwidth:               kernel_density_of_x <- density(x, bw = 5)      Then I differentiated it:               differentiate <- function(x, y)     {       diffOfX <- diff(x)       data.frame(         x      = x[-length(x)] + (diffOfX / 2),          dyByDx = diff(y) / diffOfX       )     }          first_derivative <- with(kernel_density_of_x, differentiate(x, y))      This looked just as expected:               library(ggplot2)     (p1 <- ggplot(first_derivative, aes(x, dyByDx)) + geom_line())      ![the first derivative looks smooth, as expected](http://i.stack.imgur.com/CMr4R.png) When I differentiated again, I expected another smooth curve, but I saw an odd cyclical effect.               second_derivative <- with(first_derivative, differentiate(x, dyByDx))     (p2 <- p1 %+% second_derivative + ylab("d2yByDx2"))      ![the second derivative is unexpectedly noisy](http://i.stack.imgur.com/0BzBa.png) I tried a few different options for the `kernel` argument, but the noisiness persisted.   Dropping the bandwidth down to, for example, `0.5` gave a lower frequency noise that dominated the plot (making it nonsense). Dropping the number of sampling points down from `n = 512` to `n = 32` stopped the issue, but that causes other problems. Why does this effect occur? Is it an artifact of the `density` function, or have I done something silly? * * * We can redraw the plot using the probability density function of the normal distribution that `x` was generated from to see the shape I expected:               xx <- seq.int(-20, 20, 0.1)     pdf_of_xx <- dnorm(xx, mu, sigma)     first_derivative_of_xx <- differentiate(xx, pdf_of_xx)     second_derivative_of_xx<- with(first_derivative_of_xx, differentiate(x, dyByDx))     ggplot(second_derivative_of_xx, aes(x, dyByDx)) + geom_line()      ![the second derivative created directly from the probability density function is smooth](http://i.stack.imgur.com/egoGg.png)


title: Power for experimental design

text: It's easy to calculate power for paired t-tests and ANOVA, but how do you calculate power for more complicated designs such as split plot designs or factorial designs?


title: What is the advantage of clustering (multilevel sample)?

text: (This is a question about the case study on page 636 of the book Applied Longitudinal Analysis by Fitzmaurice, Laird and Ware) The data consists of 1028 weights of mice fetuses, they are clustered in 94 mother mice. The mother mice were assigned to a dose of poison. I am puzzled as I don't see the advantage of clustering the data. Can you help me comparing the following two approaches:   * Standard linear regression model of the weight against the treatment. (blue dashed line on the graph).   * Linear mixed model with a random intercept for mother mice. (red line on the graph, with the random intercept for one mother mice displayed next to the weights of the fetuses of that mother mice).  I thought the advantage would be that the standard errors of treatment effect of the clustered model would be a lot smaller than the one of a linear model, but the opposite is true (the standard errors are 10 times as large as the ones of the linear model)! So what is the true advantage of clustering the data and adding a random intercept? ![enter image description here](http://i.stack.imgur.com/vBc8o.jpg)


title: ANCOVA or repeated measures ANOVA?

text: I have 3 groups randomised to intervention A or B or control. Participants were measured pre and post intervention - should I use repeated measures ANOVA or ANCOVA with the variable measure at time 1 (pre) as the covariate? Variables are continuous. thanks


title: Understanding forecasting in R

text: I am presently trying to learn R. I would like to be able to apply it more in my work environment as I am an analyst in the Health Care industry. I am presently trying to use R to forecast. What is the best forecasting package in R? I am presently using the forecast package. I have tried to fit the `ets` models to my data but I feel that it is giving me some fairly unreasonable solutions. The data is flat, meaning that it does not linearly increase and there are some fluctuations, but I have not been able to assess whether or not those fluctuations are seasonal. I am assuming they are not. How can I calculate the out of sample error when I am comparing forecasting models? Also, is there a way to plot my forecasted data against the actual values? Lastly, how can I determine the model that is generated from the forecast? Thanks for all of your help in advance.


title: How to calculate fitted values when using AdaBoost via gbm package in R? What does fitted values mean?

text: I poseted this question about two weeks ago. But it was closed because I cross-posted here and on SO. So I deleted the post on SO and re-post here. My question is as below. * * * First of all I'm not good at English. So please be patient to interpret my words. I really need your help. I used the R package GBM as probably my predictive modeling. I made a model with adaboost distribution. In these cases you would use `predict.gbm` code to predict some new data. But I have to predict them without R (maybe with SAS IML etc). So I want to know how to calcultate fitted values based on the model that I made. Data is as below (tab delimited). Variable 'sex' and 'houstype'(APT, DANDOK, ETC, VILLA) are categorical variables and y is a target variable that has 0,1 categories.               sex age income  area    housprice   houstype    y     M   30  72  31  99  APT 1     M   40  73  1   97  DANDOK  0     M   40  34  1   90  DANDOK  0     M   29  45  4   86  DANDOK  0     F   29  37  39  86  DANDOK  0     F   29  53  4   98  DANDOK  0     M   35  63  4   61  DANDOK  0     F   27  47  4   74  DANDOK  0     M   28  38  3   79  DANDOK  0     F   30  45  31  99  APT 0     M   31  32  4   82  DANDOK  0     M   30  47  64  98  APT 1     M   35  34  4   92  DANDOK  0     M   26  56  3   83  DANDOK  0     M   32  43  15  85  VILLA   1     M   24  32  4   69  DANDOK  0     M   31  35  3   68  DANDOK  0     M   40  47  8   52  DANDOK  0     F   26  40  56  14  DANDOK  0     F   29  45  61  63  VILLA   0     M   32  43  25  24  DANDOK  0     M   37  36  58  47  APT 1     M   36  49  49  53  APT 1     M   27  42  31  71  APT 1     M   33  53  7   48  DANDOK  0     M   29  60  8   77  VILLA   0     M   31  43  4   88  VILLA   0     M   28  35  3   83  DANDOK  0     M   38  53  12  58  VILLA   1     M   39  30  3   61  DANDOK  0     F   35  32  24  63  APT 0     F   26  36  31  97  APT 0     M   31  30  31  10  DANDOK  1     F   24  39  31  18  VILLA   0     M   38  96  4   23  APT 1     M   37  69  39  56  APT 1     M   39  37  58  17  APT 1     M   31  42  26  69  APT 1     M   32  64  64  47  APT 1     M   31  54  61  8   VILLA   1     F   29  35  18  88  APT 0     M   38  79  7   40  DANDOK  0     M   26  39  7   48  DANDOK  0     M   32  57  76  24  DANDOK  1     M   32  89  18  74  DANDOK  0     F   25  31  3   38  DANDOK  0     M   38  74  93  38  APT 1     F   34  37  76  80  APT 0     M   32  49  39  67  APT 1     M   33  60  18  69  APT 1     F   29  61  18  69  APT 0     M   37  30  83  73  APT 1     M   36  38  4   43  DANDOK  0     M   28  45  27  63  APT 1     M   39  89  39  76  APT 1     M   39  66  17  83  APT 1     M   29  67  1   61  DANDOK  0     M   30  41  1   69  DANDOK  0     M   32  44  1   83  DANDOK  0     M   32  57  49  41  VILLA   1     M   32  49  25  73  DANDOK  1     M   30  39  12  38  DANDOK  0     M   36  37  17  40  DANDOK  0     M   28  48  1   97  DANDOK  0     M   29  51  8   69  VILLA   0     M   29  50  8   69  VILLA   0     M   30  48  31  93  APT 1     M   40  86  39  98  APT 1     M   29  56  18  93  APT 1     F   29  79  8   92  APT 0     F   30  88  8   92  APT 0     M   29  34  2   69  DANDOK  0     M   34  41  25  54  DANDOK  0     M   27  34  2   92  DANDOK  0     F   34  33  49  39  VILLA   0     F   30  56  49  6   DANDOK  0     F   34  54  8   62  APT 1     F   23  37  31  38  VILLA   0     M   35  56  31  65  VILLA   1     F   26  42  39  60  APT 0     F   27  60  12  23  APT 0     F   29  52  8   37  APT 0     M   26  73  61  16  VILLA   1     F   33  45  49  39  APT 0     M   35  34  39  69  APT 1     F   25  30  8   55  APT 0     M   27  70  46  43  APT 1     M   30  43  18  81  APT 1     F   31  33  15  48  APT 0     M   32  80  24  84  APT 1     M   33  59  76  80  APT 1     M   23  32  27  32  VILLA   0     M   34  87  76  66  APT 1     F   30  42  31  46  VILLA   0     M   26  64  12  11  DANDOK  0     F   34  36  56  19  APT 0     M   35  50  4   7   DANDOK  0     M   32  75  4   7   DANDOK  0     F   29  58  69  72  APT 0     M   32  33  8   7   DANDOK  0     F   31  71  89  69  DANDOK  0     M   33  45  12  23  APT 1     M   26  56  4   76  ETC 0     M   38  48  3   44  DANDOK  0     M   37  62  39  69  APT 1     M   32  34  48  41  APT 1     M   32  54  89  20  APT 1     M   36  42  4   3   DANDOK  0     M   34  38  30  30  APT 1     M   37  51  18  38  DANDOK  0     F   27  32  24  23  APT 0     F   34  38  18  27  DANDOK  0     M   29  41  4   4   ETC 0     M   30  56  39  17  APT 1     M   32  68  8   77  APT 1     M   37  89  64  1   APT 1     M   31  39  3   57  DANDOK  0     M   38  36  83  73  VILLA   1     M   35  41  8   17  APT 1     M   33  30  49  1   DANDOK  1     M   32  32  4   4   DANDOK  0      Let the dataname as 'dataflc1'. And code is as below. Because I just want to know the method, I simplified the code.               dataflc1<-read.table("data.txt",header=TRUE)          as.character(dataflc1$sex)     as.character(dataflc1$houstype)          library(gbm)          dataflc1.output<-gbm(y~sex+age+income+area+housprice+houstype,     data=dataflc1, distribution="adaboost", n.trees=2,     keep.data=TRUE,     interaction.depth = 3)      Decision trees are as below.               pretty.gbm.tree(dataflc1.output,i.tree = 1)     pretty.gbm.tree(dataflc1.output,i.tree = 2)      which yields                 SplitVar SplitCodePred LeftNode RightNode MissingNode ErrorReduction Weight    Prediction     0        3 10.0000000000        1         2           9     17.7855736     60 -0.0002277533     1       -1 -0.0008718615       -1        -1          -1      0.0000000     24 -0.0008718615     2        0  0.0000000000        3         7           8     29.1469266     36  0.0002016521     3        2 48.5000000000        4         5           6      0.9517364     23  0.0008808467     4       -1  0.0007259475       -1        -1          -1      0.0000000     10  0.0007259475     5       -1  0.0010000000       -1        -1          -1      0.0000000     13  0.0010000000     6       -1  0.0008808467       -1        -1          -1      0.0000000     23  0.0008808467     7       -1 -0.0010000000       -1        -1          -1      0.0000000     13 -0.0010000000     8       -1  0.0002016521       -1        -1          -1      0.0000000     36  0.0002016521     9       -1 -0.0002277533       -1        -1          -1      0.0000000     60 -0.0002277533            SplitVar SplitCodePred LeftNode RightNode MissingNode ErrorReduction Weight    Prediction     0        5  1.0000000000        1         5           9      20.248252     60 -0.0001990321     1        0  2.0000000000        2         3           4      28.015752     30  0.0003333333     2       -1  0.0010000000       -1        -1          -1       0.000000     20  0.0010000000     3       -1 -0.0010000000       -1        -1          -1       0.000000     10 -0.0010000000     4       -1  0.0003333333       -1        -1          -1       0.000000     30  0.0003333333     5        3 28.0000000000        6         7           8       2.523166     30 -0.0007313975     6       -1 -0.0010000000       -1        -1          -1       0.000000     20 -0.0010000000     7       -1 -0.0001941924       -1        -1          -1       0.000000     10 -0.0001941924     8       -1 -0.0007313975       -1        -1          -1       0.000000     30 -0.0007313975     9       -1 -0.0001990321       -1        -1          -1       0.000000     60 -0.0001990321      Fitted values are as below.               dataflc1.output$fit      which yields                 [1] -0.2249587 -0.2288306 -0.2288306 -0.2288306 -0.2281529 -0.2288306 -0.2288306 -0.2288306 -0.2288306      [10] -0.2289587 -0.2288306 -0.2252328 -0.2288306 -0.2288306 -0.2272328 -0.2288306 -0.2288306 -0.2288306      [19] -0.2281529 -0.2281529 -0.2272328 -0.2252328 -0.2249587 -0.2252328 -0.2288306 -0.2288306 -0.2288306      [28] -0.2288306 -0.2269587 -0.2288306 -0.2289587 -0.2289587 -0.2264270 -0.2281529 -0.2268306 -0.2249587      [37] -0.2252328 -0.2252328 -0.2249587 -0.2261529 -0.2289587 -0.2288306 -0.2288306 -0.2261529 -0.2269587      [46] -0.2288306 -0.2249587 -0.2289587 -0.2249587 -0.2249587 -0.2289587 -0.2252328 -0.2288306 -0.2252328      [55] -0.2249587 -0.2249587 -0.2288306 -0.2288306 -0.2288306 -0.2261529 -0.2269587 -0.2272328 -0.2272328      [64] -0.2288306 -0.2288306 -0.2288306 -0.2252328 -0.2249587 -0.2249587 -0.2288306 -0.2288306 -0.2288306      [73] -0.2272328 -0.2288306 -0.2281529 -0.2281529 -0.2288306 -0.2281529 -0.2261529 -0.2289587 -0.2289587      [82] -0.2288306 -0.2261529 -0.2289587 -0.2252328 -0.2288306 -0.2249587 -0.2252328 -0.2289587 -0.2249587      [91] -0.2249587 -0.2272328 -0.2249587 -0.2281529 -0.2269587 -0.2289587 -0.2288306 -0.2288306 -0.2289587     [100] -0.2288306 -0.2281529 -0.2252328 -0.2288306 -0.2288306 -0.2249587 -0.2252328 -0.2249587 -0.2288306     [109] -0.2252328 -0.2269587 -0.2289587 -0.2289587 -0.2288306 -0.2249587 -0.2268306 -0.2249587 -0.2288306     [118] -0.2264270 -0.2268306 -0.2264270 -0.2288306      I have some questions. First, the SplitCodePred value of variable 'houstype' that has index 5 is 1.0000000000 in the second decision tree. Which category does 1.0000000000 belong to? APT? DANDOK? ETC? VILLA? How should I interpret the result? Second I want to know how to calculate 'ErrorReduction', 'Weight', 'Prediction' in each of the trees and fitted values based on the two trees. Third, What is 'MissingNode' in the results about decision trees? By the way, all variables don't have missing values. I sincerely hope your advices. Thank you.


title: how to calculate information value

text: i have the following tree               Outlook(        [Sunny] => (Yes, Yes, No, No, No),        [Overcast] => (Yes, Yes, Yes, Yes),        [Rainy] => (Yes, Yes, Yes, No, No)     )      i want to know how to calculate the information value of each set               info([2,3])     info([4,0])     info([3,2])      i already know that               info([2,3]) = 0.971 bits     info([4,0]) = 0.0 bits     info([3,2]) = 0.971 bits      but i don't know the formula * * * EDIT: i just found it now info([2,3]) = entropy(2/5,3/5) = $- (0.4 * log2(0.4) + 0.6 * log2(0.6)) = 0.971$ in general $$ entropy(p_{1},..,p_{n}) = - \sum^{n}_{i=1} (p_{i} * log2(p_{i})) $$


title: How can I prove that the cumulative distribution function is right continuous?

text: I've learned in my probability courses that the cumulative distribution function $F$ of a random variable $X$ is right continuous. Is it possible to prove that?


title: Question about prediction bands for non-linear regression computation?

text: I am interested in computing prediction bands for a non-linear regression (log-logistic function with 3 parameters). I have read the Prism help page: > The calculation of the confidence and prediction bands are fairly standard. > Read on for the details of how Prism computes prediction and confidence > bands of nonlinear regression. > > First, let's define G|x, which is the gradient of the parameters at a > particular value of X and using all the best-fit values of the parameters. > The result is a vector, with one element per parameter. For each parameter, > it is defined as dY/dP, where Y is the Y value of the curve given the > particular value of X and all the best-fit parameter values, and P is one of > the parameters.) > > G'|x is that gradient vector transposed, so it is a column rather than a row > of values. > > Cov is the covariance matrix (inversed Hessian from last iteration). It is a > square matrix with the number of rows and columns equal to the number of > parameters. Each item in the matrix is the covariance between two > parameters. > > Now compute c = G'|x * Cov * G|x. The result is a single number for any > value of X. > > The confidence and prediction bands are centered on the best fit curve, and > extend above and below the curve an equal amount. > > The confidence bands extend above and below the curve by: = > sqrt(c)*sqrt(SS/DF)*CriticalT(Confidence%, DF) > > The prediction bands extend a further distance above and below the curve, > equal to: = sqrt(c+1)*sqrt(SS/DF)*CriticalT(Confidence%, DF) I also read the topic How to compute prediction bands for non-linear regression?. But I still do not understand something : how to compute the Hessian matrix. I do not think it is (I do not get values between -1 and 1), but, is it the second-order partial derivatives of my log-logistic function ? In Prism, it refers to the "Hessian from last iteration" while in the topic, it is written "the second derivatives of the likelihood function at your estimates". The likelihood function I use is ( - _N_ / 2 ) * log( sum( ( _dat_ \- _model_ ) .^2 ) ) where _N_ is the number of observations, _dat_ is a vector with data being fitted and _model_ is a vector of model values. I do not really understand how it would be possible to do a second order derivatives of this function. Could anyone give me a hand with this ?


title: What part of an ARMA model requires a stationary time series - the AR or the MA?

text: Could I use a non-stationary time series with simply an Autoregressive model?


title: How to find the support of the posterior distribution to apply Metropolis-Hastings MCMC algorithm?

text: I am trying to sample from a posterior distribution using a MCMC algorithm using the Metropolis-Hastings sampler. How should I deal with the situations where I'm stuck in regions of the posterior with zero probability? These regions are present because the posterior distribution is truncated and also due to numerical limitations on the computer the likelihood can become zero if you are very far from the mean. That is, say the likelihood is distributed normally; if you are 100 standard deviations away from the mean you get what appears as zero probability to the computer. What I want to know is how to chose the initial value of the chain in order to be sure that it is contained in the support of the posterior.


title: Test when and how often 2 time series have similar values

text: I am monitoring two moving cars and collecting data on their instantaneous speeds three times per second for approximately two hours. I would like to create some kind of hypothesis that will tell me something about how often they have similar speeds (and if possible at what time) and how often not. Is there a statistical method that will allow me to do this?


title: How to compute conditional expectations with respect to a sigma field?

text: **Example:** Toss a coin twice. Letting $\mathbb P$ be a probability measure, suppose $\mathbb P(HH)=p^2,\mathbb P(HT)=\mathbb P(TH)=p(1-p), \mathbb P(TT)=(1-p)^2.$ I would like to answer the following questions:   1. Define $Y$ to be the number of heads in the example. Derive the $\sigma$-field generated by $Y$.    2. Find the expectation of the random variable $Y$ from the previous exercise, and also the conditional expectation of $Y$ given $\mathcal F_{1}$ s.t $\mathcal F_{1}=\\{\varnothing,\\{HH,HT\\},\\{TH,TT\\},\Omega\\}$. Check that in this case $E[E[Y|\mathcal F_{1}]]=E[Y]$.  **My answers:** the sample space is $\Omega=\\{HH,HT,TH,TT\\}$. I think the answer to 1 is $Y=\\{\\{TT\\},\\{HT,TH\\},\\{HH\\}\\}$, the $\sigma$-field generated by $Y$ is $\mathcal F(Y)=\\{\varnothing,\\{TT\\},\\{HH,HT,TH\\},\\{HT,TH\\},\\{HH,TT\\},\\{HH\\},\\{HT,TH,TT\\},\Omega\\}.$ But for question 2 I could only get $E(Y)=2p$. How is the conditional expectation computed?


title: Advice for interpolating a model

text: I'm new in Stack Exchange, so I hope no to be off topic. I'm also new in bioinformatics and I was asked to perform an analysis. Briefly, I have a dataset of 29 cell lines and the IC50 values of a test drug. I want to find a relation between the gene expression profiles of each cell line (nearly 31000 genes) and the IC50 values. My problem is the huge number of response variables (the genes) and the low number of samples. I'm trying to perform a linear regression using lasso to reduce the number of genes, dividing the samples in a train set of 14 cell lines and a test set of 15 cell lines. The division is performed by randomly sampling among the 29 samples. The problem is that lasso is not stable and every time I train the model I get different results. So I tried to reduce the dimensionality using the PCA, but as far as I have read, the PCA doesn't work well when the number of covariates is greater then the number of samples. Is this true? Can you suggest me some kind of regression which is robust when the number of samples is low? Thank you very much


title: How to calculate model error (like MSE) for a multivariate proportional response?

text: I have data where the response is multivariate and proportional (rows [observations] sum to 1). I am modelling this response using a Dirichlet regression via the **DirichletReg** R package where the predictors are the first $m$ principal components of a larger set of variables. I could use AIC/BIC to select the model size (the value for $m$), but as this is a "prediction" exercise, I'd like to be able to compute some form of error for each of the $m$ models. If this were a simple linear regression I'd use the mean squared error or its square root. In my data I have 4 proportions measured for each response and I can derive the predicted proportions for each training sample given the fitted model. What would a suitable measure of "error" be for a multivariate proportional response? ### Further details The response variable is the proportion of each of 4 vegetation types around a lake, for a set of $n$ lakes. The predictors are the first $1, ..., m$ principal components of a data set of counts of pollen grains on $p$ species. The aim is to predict the proportions of the 4 vegetation types for samples where only the pollen counts are know (like the past 10,000 years). I'd like to compute an "error" statistic based on the lack of fit to the response given the model. Then I would use a $k$-fold CV to select the "best" model for predicting the proportions of the vegetation types.


title: Frequency Distribution Test Vs. Chi Square Goodness of FIt

text: I need to prove that the higher frequency of an occurrence within a particular range is statistically significant. Let me explain this with an example: X-axis has ages and Y-axis has the frequencies of humans sleeping more than 10 hours a day.   I need to prove that people below 10 years and above 80 years of age sleep more than 10 hours a day compared to those in other age groups (such as 10-20 yrs, 20-30,.... 50-60, etc.). A simple histogram can show you higher frequency but then each group also has different universe size. So something like hitrate might not be statistically correct. I have looked into `Chi-Square Goodness of Fit` test but I'm not sure what a good choice of `expected frequency` would be in this case. Can you suggest the right approach to do this? Thanks.


title: Increasing or decreasing sample size with higher correlation among repeated measures in G*Power

text: When I conduct an a priori power analysis to calculate the required sample size with G* power for a repeated measures design with between factors, the sample size increases when I insert a higher value for correlations among repeated measures. However, when I choose the 'repeated measures, within-between interactions' variant, the opposite is true: the sample size decreases as a result of a higher correlation among repeated measures value. What am I not understandig here? For me, it would be logical that with a higher test-retest reliability, you would need fewer participants than with a lower reliability, in general. We will do a study with one within subjects variable (time: pretest vs. posttest) and one between subjects variable (condition: training vs. control). So, we are interested in the interaction between time and condition, therefore I assume we want to choose the within-between interaction option in G*power, but I am getting slightly confused here. I hope you can help. With kind regards, Johanna


title: Corrective random generation to fit a preset distribution when combined with biased data

text: I have a source of data that will have an unknown distribution. I want to generate "random" data such that, when I combine the two, the overall distribution completely hides what the original data's distribution is. It doesn't matter what I get at the end, so long as it is completely unrelated to the original. I can't simply drown the signal in noise because:   1. I can only add a finite amount of random data.   2. This needs to be repeatable a large number of times with the same starting distribution, without someone being able to infer it. The only idea I have so far is to use a probability density function that accounts for the original/target distributions to eliminate the original bias, but I'd prefer to use known methods, if they exist.


title: Regression coefficients by group in R?

text: I have a dataframe with a group variable GRP (ranging from 1-100) and an X and Y for each one. I'd like to get a list of the regression intercepts and slopes for lm(Y~X) within each group. The intercepts and slopes don't need to be in the same dataframe. Any suggestions? R beginner here, so simplicity would be great!


title: Why is it claimed that a sample is often more accurate than a census?

text: When learning the course of sampling, I meet the following two statements: 1) Sampling error leads to mostly variability, nonsampling errors lead to bias. 2) Because of nonsampling error, a sample is often more accurate than a CENSUS. I do not know how to understand these two statements. What is the underlying logic for getting these two statements?


title: Relationships between health conditions - Spearman's correlation?

text: I am examining the relationships between health issues in pigs. I want to know which conditions tend to co-occur. All animals have at least one condition and 5 at most. I have coded each condition in SPSS as 0 = absent, 1 = present. How should I examine the relationships between these? I have tried a Spearman's correlation but something is not right, maybe I have organised the data in the wrong way. This is an example of how I have organised the data:               Pig ID  Tail lesions    Lameness    Bursitis     1            1              0           0     2            0              1           0     3            0              0           1      When I am looking at the output in SPSS I can't tell if/how these conditions are related to each other when I carry out a Spearman's correlation. ![enter image description here](http://i.stack.imgur.com/qhUFY.jpg)


title: Calculating the oob strength and correlation estimates for a regression random forest

text: I'm currently working on a regression application of random forests. I'd like to produce strength and correlation graphs of the tree ensemble with varying number of random features selected (mtry in the R implementation) similar to the ones in Breiman's original paper (Figure 1 page 30): http://oz.berkeley.edu/~breiman/randomforest2001.pdf Breiman describes how to calculate the out-of-bag estimates for strength and correlation in appendix II, but only for the classification case. I was not able to deduct the calculation for the regression case. So I tried to find other resources, but couldn't find anything besides a more recent paper treating only the classification case too. http://hal.archives- ouvertes.fr/docs/00/59/84/66/PDF/icic10.pdf So I'd be really grateful if someone knows a different source. The next question would be about the implementation in R, so far i checked out the randomForest and party packages. They don't seem to provide ready made functionalities to retrieve the oob estimates for strength and correlation. I pretty sure all needed information can be subtracted from the randomForest implementation, as one basically just needs the out-of-bag indices and generalization errors per tree. I'm happy to do the coding, but I'd need the proper calculation as a starting point. cheers, yannick


title: How to generate new Topic for new documents?

text: what approach would help me generate new topics for new documents? I read this page in order to learn more about the effect of specifying keywords for the topics that we care about detecting in new documents, but the problem of generating new topics having a given set of topics is still unresolved.


title: Power of a Mann Whitney test compared to a t test

text: So a Mann Whitney U test is supposedly about 95% as powerful as a t-test when t-test assumptions of normality and homogeneous variance are satisfied. I also know that a Mann Whitney U test is more powerful than a t-test when these assumptions are not satisfied. My question is, is a Mann Whitney test on data where assumptions aren't satisfied as or almost powerful as a t-test on data where assumptions are satisfied? I'm asking because I often see people doing power calculations based on the assumption that they will perform a t test. After they collect the data, they explore data and decide to use a Mann Whitney test instead and don't really revisit how changing the test affected power. Thanks!


title: Calculate the Burning Time for a Lamp

text: If you have a lamp with burning time 4000 hours. If the time goes forward until the lamp will be destroyed the exponential distribution is 3675 hours, what is the probability of a lamp to be working at least 4000 hours? I don't know how to calculate in order to achieve the result? I have read the document about calculation but still I don't know how to do it.


title: How to check if factor structures in two groups are similar?

text: After doing semantic differential experiment with two different group of people, there is need to do Factor Analysis, and after that to compare group's coordinates in factor space. I did factor analysis separately for two groups, and unfortunately the result is their factor structure are little (or not little) different. What than should I do for compare? Need I to "merge" two groups at one file and make Factor Analysis? Or it would be wrong? I use SPSS.


title: Adjusting for experimentally-caused panel attrition when evaluating treatment effects

text: This question involves a questionable hypothetical scenario, but please bear with me. Suppose I ran an experiment in a coffee stand where the treatment was playing country music instead of the usual soft jazz while people wait. I have panel data on expenditures from the loyalty program, and I am interested in whether playing country music has a negative effect on the _subsequent_ transaction. Suppose there's no long-term effect or any fade-out of the effect to simplify things. The data for each person is at transaction level, so I see how much she spends the first time she came in, the second time, and so on. Let $i$ index people, and $t$ index transactions: \begin{equation} \overbrace{\ln y_{it}}^\text{Logged \$} =\overbrace{\alpha_i}^\text{Customer Fixed Effect} + \underbrace{\gamma \cdot C_{it-1}}_\text{Country Music at Previous Trans.} +\overbrace{x_{it}^\prime \beta}^\text{Covariates} + \underbrace{\varepsilon_{it}}_\text{Idiosyncratic Error} \end{equation} The null hypothesis is that $\gamma<0$. Consistent estimation of $\gamma$ requires eliminating $\alpha_i$ (otherwise $\hat \gamma >0$) because being a regular is correlated with treatment. I first-difference the transactions data to get rid of the fixed effect, and I am left with: \begin{equation} \ln y_{it} - \ln y_{it-1} = \gamma \cdot \left( C_{it-1} - C_{it-2} \right) + \left(x_{it} - x_{it-1}\right)^{\prime}\beta + \left(\varepsilon_{it}-\varepsilon_{it-1}\right) \end{equation} I run this regression and $\hat \gamma_{FD}<0$. I go to get an espresso. One the way to the coffee stand I realize that the people who don't return to the coffee shop are not in the data. There's no row because that transaction does not exist. If you spend less, you're in the data. If you miss a day because the memories of Garth Brooks haunt you (so the effect is -100%), you're not in the data. So the FD approach estimates an intensive margin effect, ignoring the extensive margin response. So $\hat \gamma_{FD}$ is an upper bound on the true effect. I could aggregate the data to weekly level and get some zeros that way. I can also model time between transactions. But I want the bi-marginal effect on the next transaction only. Are there any solutions to this? Hacks are perfectly acceptable if they are well-motivated.


title: Probability of continuous random variable

text: It is claimed that probability of continuous R.V. is 0 for ALL x in R. That is, for every point the probability is zero. But, somehow when we sum all these zero probabilities over the entire domain $\mathbb{R}$, the total becomes equal to 1. Can somebody explain how 0s can add up to 1? This defies kindergarden math.


title: Using AIC, for model selection when both models are equally weighted, and one model has fewer parameters

text: I am using AIC (Akaike information criterion) for model selection. There are 2 models. The first model has 2 parameters with log likelihood of -10182.0284 and the second model has 3 parameters with the same likelihood when tried on a specific dataset that displays the need for only two parameters. The weighting I get with AIC is equal for both models. The equality seems to come from the fact that not all significant figures are taken into account and with such low log likelihood the number of parameter penalization is insignificant. The results:               AICmodelSelect(-10182.0284,-10182.0284)     AIC_min     null model min AIC     relprob_null          1     relprob_alt          1     weight_null         0.5000     weight_alt         0.5000      AIC equally favours both models. I am also doing likelihood ratio test cause the models are nested and the p-value is below 0.01 for the null model (simpler constrained model). But how do I justify choosing the simpler model with AIC when there is equal weighting given here?


title: Multi-factor model for a quantitative response variable

text: I am analyzing data in order to develop a model for power generation of a solar system based on three predictor variables. I have done t-tests showing significant improvement in the generation while varying one variable at a time. Each variable has only two possible values in my situation, and for each variable I have a 95% confidence interval on the change in generation. While it is great to understand the effect of each variable alone, I would like to somehow develop a model that can predict the change in generation given the affects of multiple variables. I'll write this out symbolically if that helps explain my situation. Three predictor variables A B C Each has two possible values, 1 or 0. My response variable is G I know how G changes when I vary A whilst B and C are held constant. I know the same for when I vary B and C (holding the other two constant). I would like to know how G changes when I vary more than one of A, B, or C. What sort of statistical process would I use to create such a model? **Comment to rolando2:** After reading your article, I am not quite sure that i'm working with interactions. Let me put numbers to what I have right now to make this more concrete. One of the predictor variables is angle. I know that when changing the angle from 15 to 30 degrees, everything else being constant, I generate 80 more watts on average. This I found using a paired 2-sample t-test with daily data over a year for both angles side-by-side. Another variable is inverter type. I know that when I use inverter B instead of inverter A, I produce 40 more watts on average. I am guessing that if we put together a system with 30 degree tilt and inverter B, and compare that to our "standard" system with 15 degree tilt and inverter A, we would not generate 120 more watts on average. How would I know how much more I would produce without having to experimentally test each combination? P.S. My experience in multivariable statistics is very limited, as you probably might have guessed. I have a semester course of statistics based in R, following a year of AP Stats in HS.


title: How to check mediating effect of categorical data

text: Kindly suggest a statistical test in SPSS for the following type of model:   _X_ = dependent variable (categorical with 3 categories)   _Y_ = dependent variable (scale)   _M_ 1 = mediating (categorical with 3 categories)   _M_ 2 = moderating variable (moderator scale variable)


title: GLM analogue of weighted least squares

text: **The short version:** I can fit a model using Weighted Least Squares, given a diagonal matrix of weights $W$, by solving $(X^TWX)\hat{\beta}=X^TWy$ for $\hat{\beta}$. ~~Is there a GLM analogue? if so, what is it?~~ There seems to be a GLM analogue, e.g. with the `weights` argument in R's `glm` function. How is R using these weights? * * * **The long version:** ### the situation As a follow-up to my IPTW question, I just want to double check that I understand how to fit a parametric model using inverse probability(-of- treatment) weights (IPTW). The idea with IPTW is to simulate a dataset in which the relationship between my independent variables $(a^1,a^2,a^3)$ and dependent variable $y$ is unconfounded and therefore causal. For argument's sake let's say I already estimated an IPT weight $\hat{w}_i$ for each observation. These weights are hypothetical probability weights from the simulated dataset. ### the question I now want to fit a GLM. I'd just use WLS, but I'm working with a binary outcome and an outcome truncated at zero. So I have a linear model $\eta_i=a^T\beta$, a link $\mu_i=g(\eta_i)$, and a variance $V(y_i)$ derived from my likelihood for $y$. Then the likelihood equations are $$ \sum_{i=1}^N \frac{y_i-\mu_i}{V(y_i)}\frac{\partial\mu_i}{\partial\beta_j}=\sum_{i=1}^N \frac{y_i-\mu_i}{V(y_i)}\left(\frac{\partial\mu_i}{\partial\eta_i}x_{ij}\right)=0,~\forall j $$ as per _Categorical Data Analysis_ , Agresti, 2013, section 4.4.5. So all I have to do is multiply $var(\mu_i)$ by the weight $\hat{w}_i$, right? The same way I might if I wanted to incorporate an overdispersion parameter? If so, is this because the variance of, say, 5 independent observations is 5 times the variance of one independent observation? **Follow-up idea:** since the likelihood is the product of the likelihood for each observation, is there some weighting procedure I can use to just weight the likelihoods?


title: Probability limit formula for coefficient in errors in variables regression

text: I found an online resource which lists the plim formula for a simple regression under the classic errors-in-variables assumption as: $$ \text{plim }\beta_1=\frac{{\rm Cov}(y, x_1)}{{\rm Var}(x_1)} $$   * Why is this the case?    * Where does this formula come from?


title: Details on mutual information equation

text: I have been trying to implement mutual information in a java program that tries to identify terms that occur often together in documents. Basically, the program does the following: There is a text collection of just over 4.000.000 summaries from Wikipedia. This is searched through based on a query, resulting in a collection of say 100 texts that are deemed as a match to the query by Lucene index. I then create bigrams (structures of two words) so that all words are paired with all the other words in the result collection. Each bigram is then calculated a mutual information value for. The idea is that the highest scoring bigrams should be (hopefully) directly relevant to the original query. On to how I understand the Mutual Information equation, which for the record looks like this: P(x,y) * log(P(x,y) / (P(x) * P(y)), where P(x,y) is the probability of term x and y occurring in a document, P(x) is the probability of just term x occurs in a document, and likewise for P(y) for term y. I am unsure if I understand correctly what this means in terms of calculating the probability values. It is calculated by finding the frequency of a term divided by the number of articles in the result set. But is it the term frequency or the document frequency for a given term? Ie. is P(x) in a result set of 100 articles the number of times term x occurs in the collection divided by 100? Or is it the number of articles in which term x occurs at least once divided by 100? And if it is the former, how then is P(x,y) calculated?


title: How to compare the results of two benchmarks?

text: I am writing a program and would be interested in knowing how much faster does a version run compared to another. For that, I wrote a benchmark script that times the run time of my program, repeats the process multiple times and computes the average of the results. Let's say I run it on the version 1 of my program, I get the value `X1`. Then I run it on the version 2 and get the value `X2`. How can I tell if the difference between the two results is significant or not? In addition, how can I objectively know how many times I need to run the benchmarks to get useful results?


title: Selecting link function in GEE with binary dependent variable

text: In my experiment participants had to make a binary (yes-no) decision about various stimuli. I have two categorical (stimulus characteristics coded as -1 0 and 1 and treatment group coded as 0 1) and three continuous (questionnaire scores) variables. After useful feedback on a previous question of mine, Difference between generalized linear models & generalized linear mixed models in SPSS, I decided to analyse this dataset with Generalized Estimating Equations (in SPSS). There are several options for selecting link functions that are suitable for binary data. I ran the model with binary logistic binomial identity, binomial Logit, binomial probit and some other ones. If I compare the goodness of fit values (Quasi Likelihood under Independence Model Criterion and Corrected Quasi Likelihood under Independence Model Criterion), these are smallest when I select the link function binomial log complement. Does this mean that this is my best option? The results actually change a little bit in terms of significance. In another post I read that I should plot my data. But I don't really know what condition I should put on the x-axis, I don't have a time variable or anything like that. Is looking at the goodness of fit a good idea to do instead?


title: Hybrid (K-means + Hierarchical ) clustering

text: I have a huge dataset (50,000 2000-dimensional sparse feature vectors). I want to cluster them in to k (unknown)clusters. As hierarchical clustering is very expensive in terms of time complexity (though it provides better result), I have designed my clustering framework as follows:   1. do K-means clustering to partition the data into several bins (k is unknown so I make it reasonably large. eg. k=500)   2. get centroids of all 500 partitions   3. do hierarchical clustering on those 500 centroids(kind of merging based on some threshold value t)   4. assign the datapoints to the nearest centroid(centroids emerged from hierarchical clustering) I would like to know, whether my approach is efficient and if possible any other good solution to this problem. Thank you.


title: Bonferroni correction with Pearson's correlation and linear regression

text: I am running stats on 5 IVs (5 personality traits, extroversion, agreeableness, conscientiousness, neuroticism, openness) against 3 DVs Attitude to PCT, Attitude to CBT, Attitude to PCT vs CBT. I also added in age and gender to see what other effects there are. I am testing to see whether personality traits can predict attitudes of the DVs. I initially used Pearson's correlation for all variables (45 tests). The main finding was that extroversion was correlated to attitude of PCT at p=0.05. But as I was running 45 tests I did a Bonferroni correction of alpha = 0.05/45 = 0.001, therefore making this finding insignificant. I then ran a simple linear regression on all variables, again extroversion was significant with attitude to PCT. If I do the Bonferroni correction this it comes out insignificant again. **Questions:**   1. Do I need to Bonferroni correct at Pearson's correlation?   2. If I do, and therefore making extroversion with attitude to PCT insignificant, is there still a point in doing linear regression?   3. If I do a linear regression, do I need to do the Bonferroni correction for this also?   4. Do I only report corrected valued or both uncorrected and corrected values?


title: Understanding the intercept value in a multiple linear regression with categorical values

text: I'm failing to understand the value of the intercept value in a multiple linear regression with categorical values. Taking the "warpbreaks" data set as an example, when I do:               > lm(breaks ~ wool, data=warpbreaks)          Call:     lm(formula = breaks ~ wool, data = warpbreaks)          Coefficients:     (Intercept)        woolB          31.037       -5.778      I'm able to understand that the value of intercept is the mean value of breaks when wool equals "A", and that adding up the "woolB" coefficient to the intercept value I get the mean value of breaks when wool equals "B". However, if I also consider the tension variable in the model, I'm unable to figure out the meaning of the intercept value:               > lm(breaks ~ wool + tension, data=warpbreaks)          Call:     lm(formula = breaks ~ wool + tension, data = warpbreaks)          Coefficients:     (Intercept)        woolB     tensionM     tensionH          39.278       -5.778      -10.000      -14.722      I thought it would be the mean value of breaks when either wool equals "A" or tension equals "L", but that isn't true for this dataset. Any clues on interpreting the value of intercept?


title: Probability of selecting each item at least once when sampling with replacement

text: Sampling with replacement $N$ items from $K$, where $N \geq K$, what is the probability that each item in $K$ is selected at least once? For example, from $10,000$ individuals, what is the probability that every birthday is represented? (assume each birthday $p = 1/365$) I've been using simulations to work this out, but would like to know the general solution. Thanks!


title: How to deal with ceiling effect due to measurement tool?

text: I have collected psychophysiological data measuring the subjects' (two groups) ability to perceive vibration. A vibrating probe moves against the skin at smaller and smaller displacements, and the subject indicates when they feel the vibration. Unfortunately, at high frequencies, the probe can only move a short distance, and sometimes the largest distance that the probe can move is still not large enough for the subjects to perceive. Thus, I have accurate threshold values for some subjects, but for some who never felt the vibration, I simply have a value that I know their threshold is greater than. Is there any way for me to still include this data? And what is the best way to analyze it?


title: When does LASSO select correlated predictors?

text: I'm using the package 'lars' in R with the following code:               > library(lars)     > set.seed(3)     > n <- 1000     > x1 <- rnorm(n)     > x2 <- x1+rnorm(n)*0.5     > x3 <- rnorm(n)     > x4 <- rnorm(n)     > x5 <- rexp(n)     > y <- 5*x1 + 4*x2 + 2*x3 + 7*x4 + rnorm(n)     > x <- cbind(x1,x2,x3,x4,x5)     > cor(cbind(y,x))                 y          x1           x2           x3          x4          x5     y  1.00000000  0.74678534  0.743536093  0.210757777  0.59218321  0.03943133     x1 0.74678534  1.00000000  0.892113559  0.015302566 -0.03040464  0.04952222     x2 0.74353609  0.89211356  1.000000000 -0.003146131 -0.02172854  0.05703270     x3 0.21075778  0.01530257 -0.003146131  1.000000000  0.05437726  0.01449142     x4 0.59218321 -0.03040464 -0.021728535  0.054377256  1.00000000 -0.02166716     x5 0.03943133  0.04952222  0.057032700  0.014491422 -0.02166716  1.00000000     > m <- lars(x,y,"step",trace=T)     Forward Stepwise sequence     Computing X'X .....     LARS Step 1 :    Variable 1     added     LARS Step 2 :    Variable 4     added     LARS Step 3 :    Variable 3     added     LARS Step 4 :    Variable 2     added     LARS Step 5 :    Variable 5     added     Computing residuals, RSS etc .....      I've got a dataset with 5 continuous variables and I'm trying to fit a model to a single (dependent) variable y. Two of my predictors are highly correlated with each other (x1, x2). As you can see in the above example the lars function with 'stepwise' option first chooses the variable that is most correlated with y. The next variable to enter the model is the one that is most correlated with the residuals. Indeed, it is x4:               > round((cor(cbind(resid(lm(y~x1)),x))[1,3:6]),4)         x2     x3     x4     x5      0.1163 0.2997 0.9246 0.0037        Now, if I do the 'lasso' option:               > m <- lars(x,y,"lasso",trace=T)     LASSO sequence     Computing X'X ....     LARS Step 1 :    Variable 1     added     LARS Step 2 :    Variable 2     added     LARS Step 3 :    Variable 4     added     LARS Step 4 :    Variable 3     added     LARS Step 5 :    Variable 5     added      It adds both of the correlated variables to the model in the first two steps. This is the opposite from what I read in several papers. Most of then say that if there is a group of variables among which the correlations are very high, then the 'lasso' tends to select only one variable from the group at random. Can someone provide an example of this behavior? Or explain, why my variables x1, x2 are added to the model one after another (together) ?


title: What are acceptable validation or cross validation error rates?

text: Is there a commonly acceptable error rate for validation? As in, if the error rate is less than X %, then my machine learning method would be considered "successful". I'm looking for something analogous to a p-value of 0.05 used for many experiments, but for cross-validation. I would use 5% as an error rate but thats really hard to achieve especially if you have small training and validation sets. (I only have 6 subjects total).


title: Seeing if estimators are unbiased

text: I have the pdf $$f(y ; \theta) = \frac{1}{\theta} \exp( \frac{-y}{\theta}), \ y > 0$$ and I'm supposed to determine if the following two estimators are unbiased or not: $ \hat \theta = nY_{min} $ and $ \hat \theta = \frac{1}{n}\sum_{i=1}^n Y_i $. I'm running into some problems because when I try to find the expected value of $ Y_{min} $ and $ Y_i $, the integral is undefined.


title: Dummy coding from a bayesian stand of view

text: If we want to use categorical variables in regression context we are allowed to use dummy codings such as: http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm Is this also required in bayesian (MCMC) context, such as with WinBUGS, OpenBUGS, that we have to model k factors with k-1 dummy variables or are we allowed to use k dummy variables and linear dependency of the variables is not an issue?


title: What do I do with zero columns in a tensor.prod.model.matrix?

text: Suppose my two marginal bases are given by (for variable $p$ and $t$ the degree is equal to $1$ and $5$ inner knots are used).               > bp          p.1  p.2  p.3  p.4  p.5 p.6     [1,]   1 0.00 0.00 0.00 0.00   0     [2,]   0 0.75 0.25 0.00 0.00   0     [3,]   0 0.00 0.50 0.50 0.00   0     [4,]   0 0.00 0.00 0.25 0.75   0     [5,]   0 0.00 0.00 0.00 0.00   1      and               > bt          t.1  t.2  t.3  t.4  t.5 t.6     [1,]   1 0.00 0.00 0.00 0.00   0     [2,]   0 0.75 0.25 0.00 0.00   0     [3,]   0 0.00 0.50 0.50 0.00   0     [4,]   0 0.00 0.00 0.25 0.75   0     [5,]   0 0.00 0.00 0.00 0.00   1      Obviously defining the tensor-product-matrix as $B_i = bp_i\otimes bt_i$ I will get zero columns if, for instance, the elements of the first column of $bp$, i.e., $p.1$ are multiplied with anything which is not the first column of $bt$, i.e., $t.1$. The problem here is, that I use some quasi-Newton approach for optimization which uses a predefined gradient which involves multiplications with the columns of the tensor matrix $B$. Hence if there is a zero column, the gradient will be zero at this point and any abitrary parameter value can be used since the procedure will not change this parameter. In the optimization step, will those columns just be eliminated?


title: Confidence of a variance estimate

text: > **Possible Duplicate:**   >  Calculating required sample size, precision of variance estimate? I would like to present some measure of how variable a particular phenomenon is. This phenomenon appears to be normally distributed, so a reasonable measure of variability might be variance. I've seen people report variances before, but I've always been bothered by how no margin of error is ever reported. > **First half-question** : Why do I never see people report margins of errors > of variances? Or can you give me examples of where someone does? Or maybe > there's another measure of variability that I should use? Although I've never seen someone report a confidence measure of a variance, I feel like it would be incorrect not to report the standard error of the variability measure. So I considered how one might do that. Cochran's theorem tells us that sample variance follows a $X^2$ distribution with $n-1$ degrees of freedom. Call the sample variance $s^2$. $\frac{ns^2}{\sigma^2}\sim\chi^2_{n-1}$ Estimate the variance. Call that $\hat\sigma$. $s^2\sim\frac{\hat\sigma^2}{n}\chi^2_{n-1}$ > **Second half-question** : How do I get a confidence interval of a variance? > Do I just compute the cdf of $\frac{\hat\sigma^2}{n}\chi^2_{n-1}$ within > some bounds? And since $X^2$ is not symmetrical, is there a standard way of > getting a two-sided confidence interval? In case it's helpful, I post a formula for the standard error of a variance that I found in my code, although I don't think it's at all useful. $Var(s^2)=\frac{2\sigma^4}{n-1}$ Approximate $\sigma$ as $\hat\sigma$, and we get $\widehat{Var}(s^2)=\frac{2\hat\sigma^4}{n-1}$ $\widehat{SE}(s^2)=\sqrt{\frac{2\hat\sigma^4}{n-1}}$


title: Statistical Analysis for Perception Survey

text: ## Context: A sample of people were asked to provide ratings of their impression of jeepneys using a Likert Scale. Age, gender, household income, and so on was also obtained. **Dangerous:** strongly agree; agree; neutral; disagree; strongly disagree **Convenient:** strongly agree; agree; neutral; disagree; strongly disagree **Comfortable:** strongly agree; agree; neutral; disagree; strongly disagree **Age:** $\quad<18\quad\quad$ $19-25\quad\quad$ $25\geq\quad$ **Gender:** male; female **Household Income:** $\quad<3000\quad\quad$ $3001-5000\quad\quad$ $5000\geq\quad$ ## Questions:   * What statistical analyses would you recommend?    * How could I implement this in stata?


title: Looking for a reference to a famous quotation about precision of numbers

text: I dimly remember from my youth a quotation from a famous engineer or mathematician about the need for expressing values with a suitable number of digits, but can't remember it at the moment. It was something like > "Nothing shows ignorance faster than using more digits than is appropriate." However my google-fu has failed me. Anyone out there aware of it?


title: Standard measure of clumpiness?

text: I have a lot of data and I want to do something which seems very simple. In this large set of data, I am interested in how much a specific element clumps together. Let's say my data is an ordered set like this: {A,C,B,D,A,Z,T,C...}. Let's say I want to know whether the A's tend to be found right next to each other, as opposed to being randomly (or more evenly) distributed throughout the set. This is the property I am calling "clumpiness". Now, is there some simple measurement of data "clumpiness"? That is, some statistic that will tell me how far from randomly distributed the As are? And if there isn't a simple way to do this, what would the hard way be, roughly? Any pointers greatly appreciated!


title: Custom power analysis in R

text: Distance sampling is a common technique used to estimate the abundance of wildlife species (such as deer or turkey). This technique employees line- transect based sampling where perpendicular distance to individual animals is measured. Those distances are then used to model a detection function. That detection function is then used to correct counts of animals for incomplete detection (Buckland et al. 2001). Wildlife surveys that employee distance sampling used Program Distance (Thomas et al. 2010) to estimate abundance. Program Distance results in an abundance estimate $\hat N$ and it's coefficient of variation $cv(\hat N)$ where $cv(\hat N)$ is $se(\hat N)/\hat N$ and $var(\hat N)$ is $se(\hat N)^2$. To compare abundances from two different areas, one can use a simple z-test: $$ Z = \frac{(\hat N_1-\hat N_2)}{\sqrt{(var(\hat N_1-\hat N_2))}} $$ where $$ var(\hat N_1-\hat N_2) = var(\hat N_1)+var(\hat N_2) $$ provided $\hat N_1$ and $\hat N_2$ were independently estimated (Buckland et al. 2001). After conducting a pilot survey, it is easy to estimate the amount of effort (i.e., length of transects $\hat L$) needed to obtain a target coefficient of variation desired from future survey efforts. $$ \hat L = \frac{(L_o*cv(\hat N)^2)}{(cv_t(\hat N)^2)} $$ where $L_o$ is the total line length surveyed during a pilot survey, $cv(\hat N)$ is the coefficient of variation of the abundance estimate from the pilot survey, and $cv_t(\hat N)$ is the target coefficient of variation desired from future survey efforts (Buckland et al. 2001). To determine a target coefficient of variation, I have developed a `R` function to estimate the power of the above test to detect a given difference between $\hat N_1$ and $\hat N_2$.               # A function to estimate power to detect a change between two abundance estimates      # using a 2-tailed z-test. Where N1 and N2 are the two abundances, change is the      # desired percent change or difference, and CV is the coefficient of variation.     ZPOWER.2tail = function(change,N1,CV,alpha){       N2=N1+(N1*change)       SD1=N1*CV       SD2=N2*CV       VAR=(SD1^2)+(SD2^2)       qL=qnorm(alpha/2)       qU=qnorm(1-(alpha/2))       Z=(N1-N2)/sqrt(VAR)       power=(pnorm(qL-Z))+(1-pnorm(qU-Z))       cbind(N1,N2,change,CV,SD1,SD2,VAR,qL,qU,Z,alpha,power)}          # Example use of function.     ZPOWER.2tail(change=-0.2,N1=300,CV=0.15,alpha=0.05)      A one-tailed verison would be:               # A function to estimate power to detect a change between two abundance estimates      # using a 1-tailed z-test. Where N1 and N2 are the two abundances, change is the      # desired percent change or difference, and CV is the coefficient of variation.     ZPOWER.1tail = function(change,N1,CV,alpha){       N2=N1+(N1*change)       SD1=N1*CV       SD2=N2*CV       VAR=(SD1^2)+(SD2^2)       if(change>0) q=qnorm(alpha)       if(change<0) q=qnorm(1-alpha)       z=(N1-N2)/sqrt(VAR)       if(change>0) power=pnorm(q-z)       if(change<0) power=1-pnorm(q-z)       cbind(N1,N2,change,CV,SD1,SD2,VAR,q,z,alpha,power)}      Of course these functions could easily be combined in to one function after all the bugs are worked out. After writing this code, I found the following code posted at CREEM.               #--------------------------------------------     # Power to detect a change between two time     # points using a z-test.     #--------------------------------------------     #alpha level (type I error)     alpha<-0.05      #vector of coefficients of var     cvs <- seq(0.05,0.5,by=0.05)      #vector of true differences     fs <- seq(0,1,by=0.1)     lcvs<-length(cvs)     lfs<-length(fs)     res<-matrix(0,lcvs,lfs)     q1<-qnorm(1-(alpha/2))     q2<-qnorm(alpha/2)     s2<-sqrt(2)     #work out power at each level of CV and true diff     for(i1 in 1:lcvs) {       cv<-cvs[i1]       for (i2 in 1:lfs) {         f<-fs[i2]         z<-f/(s2*cv)         power = (1-pnorm(q1-z)) + pnorm(q2 - z)         res[i1,i2]<-power       }     }     #output as a perspective plot     persp(cvs,fs,res)     #--------------------------------------------      I know it was a long time getting here, but here is my question. The major difference between my code and the CREEM code is: $$ Z=(\hat N_1-\hat N_2)/\sqrt{var(\hat N_1-\hat N_2} $$ vs. $$ Z=d/(\sqrt{2}*cv(\hat N)) $$ Where $d$ is the "true difference." Can someone help me understand the use of $\sqrt{2}*cv(\hat N)$ in the denominator instead of $\sqrt{var(\hat N_1-\hat N_2)}$ ? Thanks. Buckland, S. T., D. R. Anderson, K. P. Burnham, J. L. Laake, D. L. Borchers, and L. Thomas. 2001. Introduction to distance sampling: estimating abundance of biological populations. Oxford University Press, New York, New York, USA. Thomas, L., S. T. Buckland, E. A. Rexstad, J. L. Laake, S. Strindberg, S. L. Hedley, J. R. B. Bishop, T. A. Marques, and K. P. Burnham. 2010. Distance software: design and analysis of distance sampling surveys for estimating population size. Journal of Applied Ecology 47:5–14.


title: Lasso modification for LARS

text: I am trying to understand how Lars algorithm can be modified to generate Lasso. While I do understand LARS, I am not able to see the Lasso modification from the paper by Tibshirani et al. In particular I don't see why the sign condition in that the sign of the non-zero coordinate must agree with the sign of the current correlation. Can someone please help me with this. I guess I am looking for a mathematical proof using KKT condition on the original L-1 norm problem i.e the Lasso. Thanks much!


title: Logistic regression: fixed effects for firms, countries & years

text: I am trying to use logistic regression on a sample of 20,000+ firms across 50+ countries, from 2000-2010. Do I need to use logistic regression with fixed effects for year and firm + dummy variables for countries? Firm Failure (fail = 1, otherwise 0) = $f(X_1, X_2, X_3, X_4, X_5, X_6)$   $X_1$ and $X_2$ are country-level independent variables   $X_3$ and $X_4$ are firm-level independent variables   $X_5$ is a country-level independent variables variable (control variable)   $X_6$ is a firm-level independent variable (control variable) Logistic regression model for sample firms $i$ and country $c$: $$\text{Fit} = β_0 + β_1X_1c + β_2X_2c + β_3X_{3i} + β_4X_{4i} + β_5X_5c + β_6X_{6i} + ε$$