title: Pesaran Cross sectional independence test in panel data

text: I was testing the cross sectional independence test using `Pesaran CD` test in `stata 10` for the panel data using the `N=50 and T=18`. After running the fixed effect with time dummies (effects), I got a negative value for the test with the `p` value >1 (which is not possible). I checked the paper on this: http://www.stata-journal.com/article.html?article=st0113 ,but they didn't talk about the negative value at all. I would really appreciate if you can help me to interpret the results.               xtreg Y  X1 X2 i.year,fe          Fixed-effects (within) regression               Number of obs      =       900     Group variable: state                           Number of groups   =        50          R-sq:  within  = 0.4336                         Obs per group: min =        18            between = 0.0417                                        avg =        18            overall = 0.3772                                        max =        18                                                          F(19,825)          =     33.24     corr(u_i, Xb)  = -0.0215                        Prob > F           =    0.0000          ------------------------------------------------------------------------------                 Y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------                X1 |  -.0008466   .0005767    -1.47   0.142    -.0019787    .0002854                X2 |    .185078   .0304043     6.09   0.000     .1253991     .244757                  |             year |            1993  |  -.0101846   .0029659    -3.43   0.001    -.0160061   -.0043631            1994  |  -.0173621   .0028749    -6.04   0.000     -.023005   -.0117192            1995  |  -.0046886   .0029062    -1.61   0.107     -.010393    .0010159            1996  |  -.0175546   .0028879    -6.08   0.000     -.023223   -.0118862            1997  |   -.009572    .002962    -3.23   0.001     -.015386    -.003758            1998  |   -.012109   .0029123    -4.16   0.000    -.0178254   -.0063926            1999  |  -.0097348   .0029043    -3.35   0.001    -.0154356   -.0040341            2000  |  -.0137918   .0028784    -4.79   0.000    -.0194416   -.0081419            2001  |    .004144   .0029323     1.41   0.158    -.0016116    .0098996            2002  |   .0188509   .0028925     6.52   0.000     .0131733    .0245285            2003  |   .0058601   .0028772     2.04   0.042     .0002127    .0115076            2004  |  -.0005801   .0028731    -0.20   0.840    -.0062195    .0050594            2005  |  -.0085907   .0029534    -2.91   0.004    -.0143877   -.0027937            2006  |   -.018118   .0028702    -6.31   0.000    -.0237516   -.0124843            2007  |  -.0164648   .0028771    -5.72   0.000    -.0221121   -.0108174            2008  |  -.0352191   .0028805   -12.23   0.000    -.0408731   -.0295652            2009  |   .0094679   .0032904     2.88   0.004     .0030094    .0159265                  |            _cons |   .0467515   .0024088    19.41   0.000     .0420234    .0514797     -------------+----------------------------------------------------------------          sigma_u |  .00686051          sigma_e |  .01433916              rho |   .1862708   (fraction of variance due to u_i)     ------------------------------------------------------------------------------     F test that all u_i=0:     F(49, 825) =     4.02             Prob > F = 0.0000               xtcsd, pesaran abs      Pesaran's test of cross sectional independence =    -2.673, Pr = 1.9925       Average absolute value of the off-diagonal elements =     0.215


title: Should the group decisions be independent in Wilcoxon rank sum test?

text: In case you want to compare the average income of a group of male employees against the average income of a group of female employees, the observations are clearly independent. Now, I have a network of a certain number of nodes. These nodes are linked by edges and I can characterize each node by the number of links it has to other nodes. (this is $k$: called the degree) I can also characterize the nodes by their average nearest neighbour degree. That is the sum of the degree of all nodes, to which one node is linked to. (this is $k_{nn}$; $k_{nn}$ of node i = $\sum k_j$ for any node $j$ that is linked to $i$). When I create a scatter plot of these nodes ($k$ vs $k_{nn}$) I can clearly distinguish two groups of nodes by a certain threshold value for $k$ and $k_{nn}$. My nodes also have a color. Now I want to test if a certain color in these two groups is overrepresented. I can do that using the Wilcoxon rank test, because the color is an independent observation. Fine. But is the color really an independent observation? Implicitly the association to a group is not only based on the node's own property, but also on the properties of the other nodes (because of $k_{nn}$). So can I really use the Wilcoxon rank test here? Actually, my question is: Does the Wilcoxon rank test require only an independent observation. Or does it also require an association to a group that is based on independent observations?


title: What is the difference between Binary Clasification and Multiclass classification?

text: Apology for posting almost one question daily. I am trying to learn some aspects of Statistical Machine learning, so every day many questions coming and if I am not finding answer in my offline peer community, I am trying to ask to you. You answer so nicely it encourages people like me to ask. I was trying to go around Naive Bayes model. It seems it supports Multi-class classification. Now there are some binary classifiers also like logistic regression. I was thinking how would be the difference of estimation between binary classification and multiclass classification.( not multinomial logistic or multilabel). If any one of the esteemed members of the group may kindly absolve my query. Thanking you in Advance, Regards, Subhabrata Banerjee.


title: Comparison of Box-Cox transformed data with different lambdas.

text: I'm trying to compare some data and see if there is a significant Pinteraction value between them. The data is highly skewed and thus I would like to use a transformation; a log transformation results in highly non-normal residuals, thus I am looking for a more appropriate transformation, if it exists. I came upon the Box-Cox transformation, and I'm trying to see if it will work. However, for every dataset I have a unique lambda, and thus a different equation of the form Box-Cox using the former, because my lambda value was found to be not zero on all occasions. My question, therefore, is if I can statistically compare two data sets transformed with different lambda values, or if there is a way to find a lambda value which is the maximum likelyhood for both data sets. Or if I've made a horrible mistake. Thank you! ~~~~~ This is how I found my lambda value, just to make sure I did not make a mistake. Assume data sets Data1 and Data2, where Data1 is the response variable.               library('MASS')              #Initial regression to get regression object     LM <- lm(Data1 ~ Data2)     LM.b <- boxcox(LM)         #x = lambda values, y = likelihood values     lam <- LM.b$x         lik <- LM.b$y     lam.lik <- cbind(lam,lik)         #Sort by likelihood to get maximum likelihood lambda      lam.lik.sort <- lam.lik[order(-lik),]     LAM <- lam.lik.sort[1,1]              #Perform regression on transformed values     Data1.trans <- ((Data1^LAM) - 1)/LAM     LM.trans <- lm(Data1.trans ~ Data2)     shapiro.test(LM.trans$residuals)      Thank you for your time!


title: Linear Matrix Inequality

text: Lyapunov's inequality states that A^T*P + P*A < 0\. Given a 2X2 A how do I go about representing this as a canonical linear matrix inequality?


title: Error message with robust regression in R

text: Whenever I run a robust regression model in R with `rlm` and with either M or MM methods, I get the following error message: > Error: 'lqs' failed: all the samples were singular I am not sure what it means and what I can do to bypass this error.


title: On using an orthogonal series to estimate a regression function

text: Suppose I have a function $g\in L_2(\mathbb{R})$, and we observe variables two -vectors $(Y_i,X_i)$ such that $Y_i = g(X_i) + U_i$ for some IID error terms $U_i$. If I want to estimate $g$, I want to use an orthonormal basis of $\mathbb{R}$, which we will call the set $\\{e_i\\}_{i\geq 0}$. Then we could write $g(x)$ as $$ \begin{align*} g(x) \sim \sum_{k=0}^{\infty} \alpha_k e_k(x), \end{align*} $$ and in the spirit of having a finite sample amount of data, we have truncate the series to a sum from 0 to $p_n$ such that $p_n \rightarrow \infty$, and $p_n = o(n)$. Then using a least squares algorithm, we obtain estimates $$ \begin{align*} \boldsymbol{\alpha}_p = (\boldsymbol{A}^\text{T}\boldsymbol{A})^{-1}\boldsymbol{A}^\text{T}\boldsymbol{Y} \end{align*} $$ where $\boldsymbol{\alpha}_p = (\alpha_0,...,\alpha_{p_n})^{-1}.$ My question has to do with the invertibility of the matrix $\boldsymbol A^\text{T}\boldsymbol A$. $\boldsymbol{QUESTION}$: How do I guarantee that $\boldsymbol A^\text{T}\boldsymbol A$ is invertible? Or more importantly, is it stochastically invertible? I know that the nonzero eigenvalues of $\boldsymbol A^\text{T}\boldsymbol A$ are the same the as the nonzero eigenvalues of $\boldsymbol{A}\boldsymbol{A}^\text{T}$, but the size of the latter matrix is smaller than the latter of the former matrix, which implies that some of the eigenvalues of $\boldsymbol{A}^\text{T}\boldsymbol{A}$ must be zero, however when I run simulations, I seem to get decent estimators using the above formula for $\boldsymbol{\alpha}_p$. So I am not sure how to justify the invertibility of $\boldsymbol{A}^\text{T}\boldsymbol{A}$. Any help would be appreciated.


title: Practically speaking, how do people handle ANOVA when the data doesn't quite meet assumptions?

text: This isn't a strictly stats question--I can read all the textbooks about ANOVA assumptions--I'm trying to figure out how actual working analysts handle data that doesn't quite meet the assumptions. I've gone through a lot of questions on this site looking for answers and I keep finding posts about when not to use ANOVA (in an abstract, idealized mathematical context) or how to do some of the things I describe below in R. I'm really trying to figure out what decisions people actually make and why. I'm running analysis on grouped data from trees (actual trees, not statistical trees) in four groups. I've got data for about 35 attributes for each tree and I'm going through each attribute to determine if the groups differ significantly on that attribute. However, in a couple of cases, the ANOVA assumptions are slightly violated because the variances aren't equal (according to a Levene's test, using alpha=.05). As I see it, my options are to: 1. Power transform the data and see if it changes the Levene p-val. 2. Use a non-parametric test like a Wilcoxon (if so, which one?). 3. Do some kind of correction to the ANOVA result, like a Bonferroni (I'm not actually sure if something like this exists?). I've tried the first two options and gotten slightly different results--in some cases one approach is significant and the other is not. I'm afraid of falling into the p-value fishing trap, and I'm looking for advice that will help me justify which approach to use. I've also read some things that suggest that heteroscedasticity isn't really that big of a problem for ANOVA unless the means and variances are correlated (i.e. they both increase together), so perhaps I can just ignore the Levene's result unless I see a pattern like that? If so, is there a test for this? Finally, I should add that I'm doing this analysis for publication in a peer- reviewed journal, so whatever approach I settle on has to pass muster with reviewers. So, if anyone can provide links to similar, published examples that would be fantastic.


title: Dealing with different time series data in Machine Learning

text: I am trying to create a stock market model based on fundamental variables for the US economy. I am using R. Some of the variables I am looking to include are: GDP, Unemployment Rate, Initial Claims, etc... Given that some of these indicators have different time frames from weekly to quarterly will I get error messages because of NA values if I try different models? Below is the output of the df I created showing all the NA values.                           Var1 Var2 Var3     Var4     2013-12-08     NA  NA 358000     NA     2013-12-15     NA  NA 368000     NA     2013-12-22     NA  NA 339000     NA     2013-12-29     NA  NA 344000     NA     2014-01-01    6.6 144     NA 317602     2014-01-05     NA  NA 333000     NA     2014-01-12     NA  NA 329000     NA     2014-01-19     NA  NA 334000     NA     2014-01-26     NA  NA 345000     NA     2014-02-01    6.7 197     NA 317760     2014-02-02     NA  NA 328000     NA     2014-02-09     NA  NA 343000     NA     2014-02-16     NA  NA 330000     NA     2014-02-23     NA  NA 351000     NA     2014-03-01    6.7 192     NA 317920     2014-03-02     NA  NA 325000     NA     2014-03-09     NA  NA 319000     NA     2014-03-16     NA  NA 323000     NA     2014-03-23     NA  NA 310000     NA     2014-03-30     NA  NA 326000     NA      Is there a model that can deal with all the different time frames (daily, weekly, monthly)? Or should I try to consolidate the time frames into the highest order (i.e. monthly)? Any insight or assistance would be greatly appreciated.   * This is my first question on this exchange.


title: Fitting data sample to a distribution

text: I'm trying to fit a data sample to a distribution. So far I have created a histogram and fitted the data with a lognormal distribution in R and made a Q-Q plot in excel (of log(benefits paid) against theoretical normal quantiles). Here are my histogram and Q-Q plot: ![enter image description here](http://i.stack.imgur.com/ZL94H.png) ![enter image description here](http://i.stack.imgur.com/sNcQu.jpg) However, I don't think this distribution is fitting the data closely enough. I was hoping for an opinion on whether there are any distributions which should fit the data better. Or is this as closer fit as I'm going to be able to get?


title: Model selection and parameter estimation in forecasting with a Dynamic Linear Model

text: I am implementing a general purpose prediction tool for time series. I want to tolerate missing values, so I decided to settle for DLMs. To make it as relevant as possible on a large number of datasets, I want it to try several different models and select the best parameters. It would then operate the prediction with the one that fits the best. This should allow me to extract as many relevant patterns as possible to make the forecasting as relevant as possible. Here is my interrogation: in most papers, meaning all the sources I've read until now, they want to use likelihood and and other similar criteria like AIC. This doesn't seem optimal in the case of forecasting. You evaluate if your model fits statistically, but that doesn't tell you if it will predict well future values. I think it would make way more sense to assess models on their predictive power within the data set. For instance you could compute the mean square error between predictions at all intermediate time points and actual realizations of your time series. This is possible thanks to the recursive nature of a DLM. With this technique you wouldn't take the risk of over-fitting, as you are assessing forecasting capacity and forecasting capacity is what you are looking for. Do you see any reasons why I would be wrong ? Why is everybody using maximum likelihood ? Have you seen any references that use something close to what I'm suggesting ?


title: Choosing a predictive model after k-fold cross-validation

text: I am wondering how to choose a predictive model after doing K-fold cross- validation. This may be awkwardly phrased, so let me explain in more detail: I understand how K-fold cross-validation works. However, whenever I run K-fold cross- validation, I use K different sets of training data, and end up with K different models. I am wondering how to go about choosing among these K models, so that I can present it to someone and say "this is the best classifier that we can come up with." Is it OK to pick any one of the K models? Or is there some kind of best practice that is involved, such as picking the model that achieves the median testing error?


title: Correction for multiple comparisons with multiple response variables

text: I'm trying to determine what constitutes a family of hypotheses when measuring multiple response variables. In particular I have an A/B type test where I'm comparing several marketing campaigns against a baseline campaign. For each campaign I want to measure increase in revenue, clicks, and purchases. What I'm not sure of is how many families of hypotheses exist in this situation? Is it more correct to group tests together based on the response variable of interest, or should all the tests be considered one big family? Concretely, assume I have 4 campaigns, and I want to compare each campaign against the baseline. Does this mean I have three families of hypothesis, one each for revenue, clicks, and purchases, and each family contains 4 hypothesis? Or is it better to say I have one family of hypothesis, containing 12 individual hypotheses (4 comparisons over three variables)?


title: Multivariate laplacian distribution

text: I am finding it strange but what does multivariate laplacian distribution looks like. What is its pdf? I googled for a while but couldn't find a good description. I wasn't paying attention to Laplacian. But all of a sudden when I need it today, I couldn't find the multivariate case.


title: Relations and differences between time-series analysis and statistical signal processing?

text: I was wondering what relations and differences are between time-series analysis and statistical signal processing? I found some recommendations of books in time series including some books in statistical signal processing. But I am not sure how these two areas are related and differ? Thanks and regards!


title: ANOVA:How to detect non-normality with a QQPlot in the presence of non-homogeneous variance

text: This is a pretty general question but, I often find statistical textbooks claiming that, in order to justify the within groups normality assumption of a one way ANOVA, you can look at a QQ plots of the residuals. However, the qq plots can only detect non-normality when the variance (or standard deviation) across all groups is homogeneous. I was wondering how I could use this type of plot when I use Welch's ANOVA. For example, in the case of Welch's ANOVA, you don't need to adhere to any assumption of homogeneous variance, but if your variance is not homogeneous, how can you use a QQ plot to test for normality? I was thinking of standardizing each of the residuals by their within group standard deviations first so that all variances would be equal. I feel like this approach is reasonable but I have not seen it discussed in any textbook or anywhere online. Lets just assume the sample sizes of each group are unequal and too small to invoke the central limit theorem.


title: Why do correlations matter in simulating data to compare classifiers when p >> N?

text: In genomics and computational biology, expression data sets contain a much larger number of features (p) than the number of observations (N). I wanted to simulate data where p>>N to compare the performance of classifier methods (e.g., regularized methods, automatic feature selection methods, and filtering/feature selection before applying any methods). In hopes of mimicking gene expression data, the informative genes were sampled from the Asymmetric Laplace Distribution (ALD) (Purdom, Holmes, 2005, suggested that gene expression data follows the ALD better than the normal distribution) and the uninformative/noisy genes were sampled from the standard normal distribution. For case 1, with n = 100 and p = 1000, 300 informative genes were sampled from X ~ ALD(-3, 0.5, 1) for class 1 and from Y ~ ALD(3, 2, 1) for class 2, and 700 uninformative/noisy genes were sampled from e ~ N(0,1). Then, 5% of the informative gene expression were transformed into outliers by increasing the scale parameter from 1 to 10. Finally, the samples (rows) were shuffled. We have E(X) = -1.93934, E(Y) = 1.93934, and Var(X) = Var(Y) = 4.761039. Here are the density curves: ![https://www.dropbox.com/s/y3n7vrybauvste0/stackexchange.jpeg](http://i.stack.imgur.com/J4iqS.jpg) The red line (left pointy hat) is X, the blue line (right pointy hat) is Y, the black line is e, the orange line is X transformed into outliers, and the green line is Y transformed into outliers. After talking to a professor about this, she suggested that I add in correlations within the features because since I am sampling independently each time, my covariance matrix for the features is the identity matrix. She suggested transforming my covariance matrix to be rho^|i-j|, where 0 My questions are:   1. **Why is it better to have the covariance matrix to be rho^|i-j| for simulation? Is that more realistic in gene expression data? Is there a specific reason why it's not all just rho in the off diagonal and 1's in the diagonal?**   2. **If I have the first 300 features to be informative and the 301th-1000th features to be noninformative (~N(0,1)), would I only add in correlations to the first 300 features? If so, would I have the covariance matrix to be {rho^|i-j| for 1 < i,j <300 if i does not equal to j, and 0 for i, j > 300 and i does not equal to j}?** This suggestion might have been related to when I told her that my CART result revealed only one different feature in R...


title: How to approximately solve system of linear equations that has no solutions

text: When I prepare data for Bayesian network meta analysis using model published in Woods B.S, Hawkins N et al - Network meta-analysis on the log-hazard scale, combining count and hazard ratio statistics accounting for multi-arm trials; A tutorial(2010), I need to solve a system of (usually trivial) simultaneous equations that might be over-specified, like the one below: $ \left\\{ \begin{array}{ll} \text{hr}_1 = 0 \\\ \text{hr}_2 - \text{hr}_1 & = \text{HR}_{2,\,1}\\\ \text{hr}_3 - \text{hr}_1 & = \text{HR}_{3,\,1}\\\ \text{hr}_3 - \text{hr}_2 & = \text{HR}_{3,\,2} \end{array} \right. $ where $\text{HR}_{i,\,j}$ is a published hazard rate, and $\text{hr}_k$ is latent cumulative hazard that I want to estimate. What people usually do, is ignoring random equation(s) from the above set and solve the rest algebraically. But I think that a better solution would be to systematically (not randomly) produce all distinct and algebraically solvable truncated equation sets, and than average all the results. It is doable, of course, and I can program it in e.g. R. (I'd need to take all possible permutations of truncation of this equation, decide which is solvable, solve them, and average the results). In the abundance of the packages already existing in R, maybe someone has already published functions that I might use to this end? The problem is, that I don't know how to find them. Maybe at least you can help my creativity and help me with the search phrase I could put to Google?


title: Evaluating parametric vs non-parametric methods

text: I am having a hard time finding comparisons between non-parametric and parametric methods, specifically for the task of density estimation (e.g. GMM vs using Dirichlet Processes). More than tractability or running time, I am interested in their statistical performance on general datasets. Are there any known studies or papers on this topic?


title: Choice between Type-I, Type-II, or Type-III ANOVA

text: We have a dataset with three variables (dV: self-reported measure on scale 1-5, assumed to be metric; iV1: factor with 4 levels; iV2: factor with 8 levels). We are interested whether the dV differs in regard to both iVs and whether there is an interaction between the iVs. Idea: Calculating an ANOVA with both main effects and the interaction between both iVs using R. The question: What Type of Sum of Squares should be used for this research question? Using aov() in R calculates Type-I Sum of Squares as standard. SPSS and SAS, on the other hand, calculate Type-III Sum of Squares by default. However, using Anova() {car} in combination with options(contrasts=c("contr.sum", "contr.poly")) in R gives the same Type-III ANOVA tables as calculated in SPSS. I have already read the following discussions:   * http://afni.nimh.nih.gov/sscc/gangc/SS.html   * http://myowelt.blogspot.de/2008/05/obtaining-same-anova-results-in-r-as-in.html However, I am still confused which Type of Sums of Squares is the most adequate for our question. The results (F and p values) differ considerably.


title: Estimating Change in Mean (Linear Regression Example)

text: A question states: Estimate the mean change in the number of flycatchers killed when the nest box tit occupancy increases by 10%. The slope for this regression is .10766X. Hence, says the answer, take .10766(10) = the mean change. My question is: The reason they are inserting 10 is not because it is a 10% change, correct? I cannot find another explanation but I am confident that this can't be the reason. Once I know for certain I can assess other possibilities.


title: What statistical methods are archaic and should be omitted from textbooks?

text: In answering a question about a confidence interval for a binomial proportion I pointed out the fact that the normal approximation is an unreliable method that is archaic. It should not be taught as a method, although there might be an argument that it be included as a part of a lesson about what makes an adequate method. What are other 'standard' statistical approaches that have passed their use-by date and should be omitted from future editions of textbooks (thereby making space for useful ideas)?


title: Factor Analysis: calculate maximum likelihood factor loadings from only the correlation (R) matrix and/or covariance (S) matrix?

text: Does anybody know how to calculate the maximum likelihood factor loadings from _only_ the correlation (R) matrix and/or covariance (S) matrix in Factor Analysis "by hand" (i.e., by Excel)? Or, even better, point me to a clear explanation with a worked example? I don't have the underlying data, so I can't just use a software program to solve the problem for me. (Any answer involving typical FA software commands or menu sequences won't be helpful, I don't think.) To the extent it's helpful, I'm using Johnson and Wichern's "Applied Multivariate Statistical Analysis", 6th Edition. In the exercises in Chapter 9, we're given just the R or S matrix and then asked to perform the maximum likelihood calculation repeatedly (e.g., exercises 9.20 b; 9.24, 9.26, and 9.27), yet the book never shows an example (as far as I can tell, nor does the answer key, which my professor makes available to us). Thanks in advance and regards!


title: First stage of TSLS and the matrix of instruments W

text: If we assume we have 2 equations and each equation contains the other dependent variable. $y_1 = \beta_0 + \beta_1 y_2 + \beta_2 z_1 + u_1$ $y_2 = \alpha_0 + \alpha_1 y_1 + \alpha_2 z_2 + u_2$ For exact identification, we would need two instruments. Let's say these instruments are $z_3$ and $z_4$. The first stage of TSLS is equal to the regression of the assumed endogenous variables $y_1$ and $y_2$ on the matrix of instruments $W$ which contains all the instruments and exogenous variables. So the columns of the matrix of instruments would look like $W = [z_1, z_2, z_3, z_4]$ and the first stage of TSLS like $y_1 = W\pi_1 + v_1$ $y_2 = W\pi_2 + v_2$ For identification, I need the significance of exact one of those instruments, i.e., of $z_3$ and $z_4$ in each of those reduced form equations. **My question is: does matrix $W$ indeed contain both instruments or do we regress onto $W_1=[z_1, z_2, z_3]$ and $W_2=[z_1, z_2, z_4]$ if we take $z_3$ as a instrument for $y_1$ and the instrument variable $z_4$ for the endogenous variable $z_4$?** **Is that right that we would need at least 3 instruments if we observe something like?** $y_1 = \beta_0 + \beta_1 y_2 + \beta_2 y_3 + \beta_3 z_1 + u_1$ $y_2 = \alpha_0 + \alpha_1 y_1 + \alpha_2 y_3 + \alpha_3 z_2 + u_2$ $y_3 = \gamma_0 + \gamma_1 y_1 + \gamma_2 y_2 + \gamma_3 z_3 + u_3$ so that the matrix of instruments would look like $W = [z_1, z_2, z_3, z_4, z_5, z_6]$ if $z_4, z_5$ and $z_6$ are our instruments.


title: Ridge regression coefficients

text: When running ridge regression, how do you interpret coefficients that end up larger than their corresponding coefficients under least squares (for certain values of lambda)? Isn't ridge regression supposed to monotonically shrink coefficients? On a related note, how does one interpret a coefficient whose sign changes during ridge regression (i.e., the ridge trace crosses from negative to positive on a ridge trace plot)? Thanks!


title: DOE design matrices and orthogonality

text: Some time ago I read this great short paper on the relationship between 'correlation', ‘independence’ and ‘orthogonality’. In short, two vectors can be uncorrelated, they can be orthogonal, they can be both, or they can be neither. I have noticed in differing regression textbooks the authors state different requirements for ensuring that coefficient values are not affected by the inclusion of additional independent variables. For example, in the textbook _Introduction to Econometrics_ by Stock & Watson the authors state that if regressor variables are uncorrelated with one another then their coefficient estimates will be invariant with respect to the inclusion of additional independent variables. In the book _Introduction to Linear Regression Analysis_ by Montgomery the author says that regressor variables must be orthogonal to one another in order to ensure that their coefficient estimates are not affected by one another. The professor who taught the class that used the second book says that when the book described orthogonality as being necessary, it was referring to the centered versions of the regressor variables. If that is true then it is in fact _uncorrelated_ variables that have invariant coefficient estimates - at least with respect to the inclusion of additional variables, as the first book suggests. When studying Design of Experiments one learns that the design matrices are specifically created in such a way that the coefficient estimates associated with a regressor variable are invariant with respect to the other variables used in the analysis. Interestingly, though, these design matrices don't just ensure that the IV's are uncorrelated (which is all that is necessary to achieve the result), they _also_ ensure orthogonality. **My question is: is there a greater purpose for ensuring orthogonality in addition to 'uncorrelated-ness'? And/or is it just too difficult to ensure 'uncorrelated-ness' without also ensuring orthogonality?**


title: Accessing errors of 3D surface generation algorithm

text: > **Possible Duplicate:**   >  Assessing error of a spatial interpolation algorithm This is a similar question to this and this one. I have a set of 3D points that are sampled from a terrain. The points have three components $x$, $y$, $z$. You can think of these points as the surveyor points that one collects from measuring a terrain for GIS purposes. I have a few algorithms that take in this set of 3D points as the reference: the algorithms will then   1. decide on the mesh size   2. compute a polygonal mesh for the area of the convex hull of the input 3D points. Now, given that I have a few algorithms for this, is there anyway I can evaluate and compare the correctness of these algorithms? Is there any kind of formula that allows me to know which algorithms are better at generating the 3D surface? The two most intuitive criteria for the evaluation are:   1. How near the interpolated points are to the given input points   2. How smooth is the surface  I have no idea how to quantitatively define the above two, what do you think?


title: simulate dirichlet process in R

text: I am reading the paper of "Dirichlet Process Mixtures of Generalized Linear Models" authored by L. A. Hannah. If I would like to simulate the following model $$\mathcal{P}\sim \text{DP}(c\mathbb{G}_0)$$ $$\theta_i|\mathcal{P}\sim\mathcal{P}$$ $$X_{i,j}|\theta_{i,x}\sim\mathcal{N}(\mu_{ij},\sigma^2_{ij}), j=1,...d$$ $$Y_i|X_i,\theta_{i,y}\sim\mathcal{N}(\beta_{i0}+\sum^d_j=\beta_{ij}X_{ij},\sigma^2_{ij})$$ In R, how can I get $\mathcal{P}$ and $\theta_i|\mathcal{P}$ in $$\mathcal{P}\sim\text{DP}(c\mathbb{G}_0)$$ $$\theta_i|\mathcal{P}\sim\mathcal{P}$$


title: In bivariate linear regression is there a direct relationship between $n$, $r^2$ and coefficient error?

text: In bivariate linear regression is there a direct relationship between sample size $n$, coefficient of determination $r^2$ and $\sigma_\beta$ (the standard error of coefficient $\beta$)? Assume data have been normalized so both target and predictor variable have $\sigma=1$. Putting the question another way, does $\sigma_\beta$ tell me something different to $r^2$ or are they measures of the same thing? Or, is it possible to have a strong, certain but unreliable link between two variables (large $\beta$, small $\sigma_\beta$, but small $r^2$)? (In multiple regression this doesn't apply as even with high $r^2$, $\sigma_\beta$ can indicate uncertainty as to which of the multiple predictors is causing the response). **EDIT** Just got this out of my software (without standardized data):               regression coeff 0.023     stderr of coeff 0.0046     p=0.000002          n=131     multiple r2=0.17     adjusted r2=0.16          predictor std=22.5     target std=2.24      Standardized coefficient is presumably $0.023*22.5/2.24 = 0.23$. If standardized coefficient is the same as correlation, then $r^2 = 0.23^2 = 0.053$ ...not the same as the software gave. What am I doing wrong?


title: Comparing the performance of two models

text: Suppose that I have a set of data sets $(m_i)$. Each of this data set $m_i$ is related to a parameter $\hat{\theta}_i$ that I know. Then I have two models providing me HPDs intervals for each $\theta_i|m_i$. I want to assess/compare the performance of the two models. Obvisouly I can simply look how many of the $\hat{\theta}_i$ falls in the intervals for both method but this seems to be a very rough characterisation (the two models can give identical scores but HPDs of very different nature) and I wonder is there any more specific tools to do so ?


title: Goodness of fit using frequency tables with low counts

text: I have two frequency tables, one representing observed data and one representing modeled data. I am looking for a Goodness of Fit measure, checking whether the model data fits the observed data. Problem is, my counts are rather small (most are smaller than 5), and so Pearson Goodness of Fit fails to operate on these data sets (I am using R). R reports the following:               X-squared = Inf, df = 534, p-value < 2.2e-16      Should I use another Goodness of Fit measure? Any suggestions? **EDIT** (providing some more info): The _model_ contains frequencies of "topics" that are associated with "documents" read by users in the last month. The _observed data_ contains "topics" that are associated with "documents" read by users in the last day. (If "documents" are Music CDs, "topics" are musical Genres, for example) I am trying to understand if the model (data over time) fits the observed data (data of last day). If the model fits the observed data, it should mean that users stay around the same "topics" (keep listening to the same musical genres).


title: Linear mixed model for repeated measures analysis with missing values

text: we have ran an experiment where we compared three interaction techniques for a 3d docking task. So we had two factors: the aforementioned technique type and a factor representing the direction translation (i.e.: if participants had to move an object that appeared close to their viewpoint and move it in depth or vice versa). Each trial was repeated 5 times. Some of those trials were skipped because of the difficulty. If I run a regular repeated measures anova then each participant that even a single missing value will be dropped from the analysis. This means that I'd have to remove more than half of the participants. By reading around it seems I can use a linear mixed model instead. My doubt is, can I use a mixed model for this type of situation? I am confused as to whether mixed models are only relevant when you have a between-subjects factor such as the classic treatment/control groups. In my case every participant was subjected to the same conditions. There were no between- subjects factors. I ran the mixed model analysis by using Technique, Direction and Repetition as _Repeated_ , the ID of each participant as the _subject_ and technique and direction as fixed factors. Are my assumptions correct or did I do a terrible mistake? If so, what alternatives do I have when dealing with missing values? Thanks!


title: Uniform updating

text: I have a uniform prior with certain values $a$ and $b$ (not standard uniform). How do I update this distribution to take into account the results from my data? If it was $U(1,1)$, I could convert it to a beta distribution, but since it is not, I'm not sure how to approach this.


title: Validity of proportion estimate from finite sample

text: If I have a batch of 130 parts and I inspect a random sample of 14 parts, either accept or reject the part, and calculate the percent defective in the sample, how confident can I be that the percent defective in the sample represents the percent defective in the batch? Here's why I need to know: We have a process with very poor capability - generally around 20% defective. We want to pull a sample as soon as a batch of 130 parts is done and be able to predict the yield with some confidence level so we can react by either continuing to run if less than 20% defective would be expected from the batch within our confidence interval or shut down if greater than 20% defective is expected based on the sample. I will be asked how confident am I that the sample yield represents the final batch yield. P.S. I am not a statistician, but I have taken a lot of engineering mathematics.


title: PROC KDE convergence problem

text: I'm using PROC KDE to fit a density to a large set of continuous data. I can't figure out how to get around the convergence issue below. Any suggestions? I could plot percentiles, but I prefer the densities. Here's my log:               285  ods graphics off;     286  ods graphics on;     287  proc kde data=sasdata.have;     288    univar x(bwm=2) x(bwm=0.25)/ plots=(histdensity);     289  run;          NOTE: Processing UNIVAR statement number 1.     WARNING: Algorithm for Sheather-Jones Plug In did not converge; using     maximum bandwidth value.     WARNING: Algorithm for Sheather-Jones Plug In did not converge; using     maximum bandwidth value.     NOTE: PROCEDURE KDE used (Total process time):       real time           14.63 seconds       cpu time            5.67 seconds               290  ods graphics off;


title: Hypothesis testing on data with overlapping groups

text: I frequently run into the following example: Lot's of genes (order of hundreds of thousands), a few hundred samples stratified into multiple groups. And those groups overlap (i.e. healthy/sick + young/old or something similar). The task - I need to run a hypothesis test for differences in means for each of the genes. So I normally do t-test with multiple-comparison correction of p-values. However now if I would do that - I would mix some groups. Ideally I would have to run t-test for all different combinations of two groups (healthy youngs vs sick youngs; healthy olds vs sick olds). But this will reduce my already low sample size even more. So I am curious for suggestions - what is the typical way to approach this situation? ANOVA taking each combination of groups as separate group?


title: Using Pearson's correlation coefficient for probability

text: After a sample size of 400+ I was able to get a Pearson's coefficient of .25. How am I supposed to break this down into a probability or a percentage. Rather, how can I explain my findings in laymen terms? I should give some more information. We have two different tests. One of these tests has 1 question, How satisfied are you? It is scored between a 1-4, with 4 being the most satisfied and 1 being the lowest. The other test has 18 questions. These were an internal review. The review consisted of such questions (did the customer service rep use the customer's name, did the customer service rep give the correct technical answer, did the rep teach the customer how to fix the problem themselves if applicable). The first test is asked to the customer, the second test with its 18 different questions is filled out by the supervisor. They either pass (-1), Not Applicable or Neither Pass/Fail (0), or Fail (1). Our goal is to find which variables in the second test best predict scores in the first test. We want to do this to find the problem areas that we can fix internally in order to better serve our customer and give them the best experience possible (More 3's and 4's on the first test).


title: Comparing Incidence Rate Ratios

text: I am trying to compare incidence rate ratios (IRR). This is what I have: Group A (exposed vs. unexposed) IRR and Group B (exposed vs. unexposed) IRR So how would I correctly test for the difference between the IRR from Group A and the IRR from Group B? The question I would like to answer is whether the incremental differences observed in Group A are similar to those observed in Group B? Please help!


title: What do I need to know about Bernoulli distributions to build a Naive Bayesian classifier?

text: I'm building a naive bayesian classifier for a binary classification. Right now I have an estimator for Bernoulli distributions, and real distributions (using a kernel mixture distribution). I can build my distributions nicely, and take the PDF's from my distribution and samples. But when I combine the probabilities together, my bernoulli probabilities simply overwhelm all other evidence. Even for datasets I know have zero chance of being members of some class, they end up having probabilities in the 90-99% range. I think I'm being naive (hah!) here, am I missing something? EDIT: This is a little long for a comment, so I've attached it to my question. My data looks like a Needle in a haystack problem - Very few in my True class (about 3-5%), most in the false class. So my false probability is fire-walled, which is fine. But my True class is also overestimating, which is bad. I'm okay with NBC not being accurate wrt probabilities - just the overall classification. I get these extreme results when I start including my boolean values. Which is what's bothering me. Perhaps I'm missing some understanding of how to implement this. To give some idea on my training set, I've got 20,000 training samples, 10,000 in the true class, 10,000 in the false class. I'm predicted around 400K rows. If I use a simple probtrue > probfalse estimate, I get 126 out of 400K as being in the true class (true value is a little more than 10,000). If I simply us probtrue > 0.99, that shrinks down to about 84K. an improvement, But still not what I was expecting. If I can shrink this set down to 40K, I'll be ecstatic. FYI: I was previously using Prior Knowledge's prediction engine using which I was able to get down to a reasonable accuracy. Then they got bought up by Salesforce and made their CEO Director of predictive analytics.


title: What is $EY$, if $Y=max(X_{1},X_{2},...,X_{n})$ where $X_{i}$ are observations from uniform distribution over $(0,a)$

text: What is $EY$, if $Y=max(X_{1},X_{2},...,X_{n})$ where $X_{i}$ are observations from uniform distribution over set $(0,a)$, $EY$ goes to $a$ as $n$ goes to infinity ?


title: Moralization and triangulation on belief networks

text: Assume that I have a belief network with a set of nodes. In order to create a valid junction tree I have to moralize the graph. Assume now that I have nodes with more than 2 parents (e.g 3 parents) then I have to link each of the parents. After doing that, the graphical models contain some unnecessary cliques. Can I use only a part of the cliques for the Junction Tree Creation? Or I must use every clique in the junction tree (e.g., a case with 8 nodes and 8 cliques)?


title: Using sparse inverse covariance matrix in estimating least squares coefficients

text: I am reading the paper introducing the graphical lasso, which is a way to estimate a sparse inverse covariance matrix. http://www-stat.stanford.edu/~tibs/ftp/graph.pdf Finding a sparse inverse covariance matrix can be useful because 0's in the inverse of the covariance matrix imply that two variables are independent. In this paper, it is mentioned that it is related to the least squares solution because $$B = (X^T X)^{-1} X^T Y$$. In this solution $(X^TX)^{-1}$ can be thought of as the inverse covariance matrix. I am curious about two things: 1) Say I have some data where I suspect that variables of $X$ are independent of some other variables in $X$. Perhaps computing the empirical covariance matrix of $X$ doesn't reflect that idea, so instead I try to find a sparser inverse covariance matrix using the graphical lasso. Now, I introduce a response variable $Y$ and want to model a linear relationship that takes into consideration the idea that some variables in $X$ are completely independent of other variables in $X$. Does this mean I can solve a "graphical" least squares problem by taking a response matrix $Y$ and a design matrix $X$, and computing the sparse inverse covariance matrix of $X$ (call this $W$) then multiply by the crossproduct $X^T Y$, ie $$B = W X^T Y$$? When computing the crossproduct $X^T Y$, should I still use the original design matrix $X$, or should it be some function of $W$ now? 2) I am curious about the case where the response is a tuple, maybe the response $Y$ is a three column matrix. Suppose I want to measure a linear model for $Y$ with a design matrix $X$ again, but this time I just want sparsity in the regression coefficients, taking us back to the standard lasso problem. I could iterate over each column of $Y$ (call these $j_1, j_2, j_3$), and run a lasso with $j_i$ and $X$ which would yield 3 coefficient vectors, each vector having sparsity. Maybe instead I would like to estimate an entire coefficient matrix, where the matrix as a whole is sparse but each column individually need not be sparse. The paper has this to say: (does this mean solving for a sparse B matrix is more appropriately done using a graphical lasso?) > Blockquote There is a simple, conceptually appealing way to view this > procedure. Given a data matrix X and outcome vector y, we can think of the > linear least squares regression estimates $(X^TX)^{-1}X^TY$ as functions not > of the raw data, but instead the inner products $X^TX$ and $X^TY$ . > Similarly, one can show that the lasso estimates are functions of these > inner products as well. Hence in the current problem, we can think of the > lasso estimates for the pth variable on the others as having the functional > form lasso(S11; s12; ): (9) But application of the lasso to each variable > does not solve problem (1); to solve this via the graphical lasso we instead > use the inner products W11 and s12. That is, we replace (9) by lasso(W11; > s12; ): (10) The point is that problem (1) is not equivalent to p separate > regularized regression problems, but to p coupled lasso problems that share > the same W and $\theta = W^{-1}$ . The use of W11 in place of S11 shares the > information between the problems in an appropriate fashion. Blockquote


title: Interpreting a binned residual plot in logistic regression

text: I am carrying out a logistic regression with $24$ independent variables and $123,996$ observations. I am evaluating the model fit in order to determine if the data meet the model assumptions and have produced the following binned residual plot using the `arm` `R` package: ![enter image description here](http://i.stack.imgur.com/cS39X.png) Obviously there are some bad signs in this plot: many points fall outside the confidence bands and there is a distinctive pattern to the residuals. My question is - can I attach these issues to specific assumptions of the logistic regression model? For instance, can I say that there is evidence of non-linearity in the independent variables or of heteroscedasticity? If not, are there other diagnostics I can produce to help identify where the problem lies? * * * Based on Daniel's answer, it appears that the main issue is I was using residuals on the logit scale but expected values on the response scale. If I reproduce the plot with the residuals also on the response scale it looks like this: ![enter image description here](http://i.stack.imgur.com/rgyyD.png) which is much more believable.


title: Confidence interval for measured value with error

text: I have the (inexact) measurement of N points. I also have access to the true value for M < N of such points. Suppose the measurement error is normally distributed. Based on this knowledge I could probably compute the mean ("error_mean") and SD ("error_sd") of the measurement error based on the M points for which I have both the true value and the measurement. However I am not sure how to do this. See below what I need exactly. My objective is: considering one of the N-M points for which I do not have access to the true value, I want to give a 95% confidence interval for such value. I think that the confidence interval should be something around:               CI_lower_end = measurement -error_mean -2*error_sd     CI_upper_end = measurement -error_mean +2*error_sd      However I suspect that error_mean and error_sd are not simply the mean and SD of the M values that result from performing "measurement - true_value" of the M points above. In particular, if M is low such as 2, the error_sd would probably be much underestimated with such naive approach. Can you tell me how to compute the abovementioned confidence interval? Thank you.


title: Regression with repeated measures in minitab

text: I do not have a strong statistics background, and I am trying to figure out the best way to determine whether a correlation exists between two variables $x$ and $y$ using Minitab. I have repeated measures for both variables for each of 8 different subjects, but for any given subject the $n$ for variable $x$ is different from the $n$ for variable $y$. Similarly across subjects the $n$ for $x_i$ is different from $n$ for $x_j$, and $n$ for $y_i$ is different from n for $y_j$. I suppose I could just calculate the means for each $x_i$ and $y_i$ separately and perform a regression using those, but it feels like that would be overly simplistic and might lead me to an invalid conclusion. Is there a better direction I should be heading down?


title: Difference of Gamma random variables

text: Given two independent random variables $X\sim \mathrm{Gamma}(\alpha_X,\beta_X)$ and $Y\sim \mathrm{Gamma}(\alpha_Y,\beta_Y)$, what is the distribution of the difference, i.e. $D=X-Y$? If the result is not well-known, how would I go about deriving the result? Thank you!


title: How to best handle subscores in a meta-analysis?

text: I am conducting a meta-analysis of effect sizes **_d_** in R using the metafor package. **_d_** represents differences in memory scores between patients and healthy. However some studies report only subscores of the measure of interest **_d_** (e.g. several different memory scores or scores from three separate blocks of memory testing). Please see the following dummy data set with **_d_** representing the studies' effect sizes as well as their standard deviations sd.:               d <- round(rnorm(5,5,1),2)     sd <- round(rnorm(5,1,0.1),2)     study <- c(1,2,3,3,3)     subscore <- c(1,1,1,2,3)     my_data <- as.data.frame(cbind(study, subscore, d, sd))          library(metafor)     m1 <- rma(d,sd, data=my_data)     summary(m1)      I would like to ask for your opinion for the best way how to handle these subscores - e.g.:   1. Select one subscore from each study that reports more than one score.   2. Include all subscores (this would violate the assumption of independency of the rfx model as subscores of one study come from the same sample)    3. For each study that reports subscores: calculate an average score & average standard deviation and include this "merged effect size" in the rfx meta-analysis.    4. Include all subscores and add a dummy variable indicating from which study a certain score is derived.


title: Lasso regression with GAM models in R

text: I'm using a regression package that uses `gam()` from the `mgcv` package. Is it possible to include an L1 penalty with the `gam` function? The documentation indicates that an L2 penalty can be used using the _H_ parameter, but I do not understand how, or if it can be used for the lasso as well. Splines are not in the model; the only reason why I'm not using `glmnet` is because the function I use has GAM implemented. Example code:               formula.eq = Y ~ X     model.out = gam(formula.eq, binomial(link = "logit"), gamma = 1)      Where `Y` is a vector of zeros and ones, and `X` includes > 400 features, so I need to perform some shrinkage method. Thanks.


title: What statistical analysis is most appropriate for a cross-over treatment design? How many participants?

text: I am asking on behalf of a colleague. I was hoping to have an answer for him, but I rather seek some guidance and be a bit more confident. He is designing a treatment study, which I am inventing here for illustrative purposes. The treatment (singing) is supposed to increase mood, while the control (humming) may increase mood, but we are not sure. Mood will be measured by a continuous variable. N participants will be randomized to one of two groups: baseline, singing, humming control (ABC), or baseline, humming control, singing (ACB). The experimental phase (BC or CB) of each group will consist of three weekly sessions for six weeks. In other words, after baseline, participants will either complete three weeks of singing, followed by three weeks of a humming control, or vice-versa. He expects that the ABC group will show improvements early on (during the singing phase), which will either be maintained or slightly revert toward baseline during the humming control phase. The ACB group is expected to show little change from baseline until they begin to sing. Questions: **What statistical analysis is most appropriate for the design?** (I was thinking a simple paired samples t-tests design, wherein we would compare A to B, B to C, and A to C for each group. Please correct me.) **How many participants would be needed, assuming a "moderate" effect size?** Moderate can perhaps be conceptualized as a change of .5 to .7 standard deviation on the continuous dependent variable. I referred him to two books (DESIGN AND ANALYSIS OF CLINICAL TRIALS and STATISTICS APPLIED TO CLINICAL TRIALS) but they are pretty technical, and a discussion may be more helpful. He cannot have too many participants, and hopes that the power from the cross- over design will be helpful. Hopefully 10 or so is enough, but he is unsure. Any help or guidance would be appreciated!


title: Linear Regression and ANOVA

text: I found two very useful posts about the difference between linear regression analysis and ANOVA and how to visualise them: Why is ANOVA taught / used as if it is a different research methodology compared to linear regression? How to visualize what ANOVA does? As stated in the first post, to test whether the average height of male and females is the same you can use a regression model ($y = \alpha + \beta x + \epsilon$, where $y$ denotes height and $x$ denotes gender) and test whether $\beta = 0$. If $\beta = 0$, then there is no difference in the height between males and females. However, I am not quite sure how this is tested when you have three groups. Imagine the following example:               height (y) -  group (x)     5          -  A     6          -  A     7          -  A     6          -  A     30         -  B     32         -  B     34         -  B     33         -  B      20         -  C     19         -  C     21         -  C     22         -  C      The regression model would look like: $$y = a+ b x + \epsilon$$ I quickly visualized the data (see image below) They way I understood the regression model is that it would now test whether any of the three slopes (AB, AC or BC) has a slope $b$ which is significantly different from 0. If that's the case one can conclude like in an ANOVA that there is at least one group in which height is significantly different from one or more groups. Afterwards, one could use a post-hoc test of course to test which of the groups really differ. Is my understanding of how the regression models tests this hypothesis correct? ![enter image description here](http://i.stack.imgur.com/6LD5Q.png)


title: Summarization of text documents (legal domain) using deep learning techniques

text: I am referring to the site deeplearning.net on how to implement the deep learning architectures. I have read quite a few research papers on document summarization (both single document and multidocument) but I am unable to figure to how exactly the summary is generated for each document. Once the training is done, the network stabilizes during testing phase. So even if I know the set of features (which I have figured out) that are learnt during the training phase, it would be difficult to find out the importance of each feature (because the weight vector of the network is stabilized) during the testing phase where I will be trying to generate summary for each document. I tried to figure this out for a long time but it's in vain. If anybody has worked on it or have any idea regarding the same, please give me some pointers. I really appreciate your help. Thank you.


title: Hidden Markov Model: How to deal with multiple symbols emitted from a given state

text: I would like to build a HMM that would work similarly to a parts-of-speech tagger. But rather than tag words in a sentence with their part of speech, I would tag words in user searches on e-commerce sites with their category. For instance, both "red party dress" and "black formal tie" would be tagged with `[color] [occasion] [item]`. I have a basic understanding of my approach with one exception: How should I deal with multi-term elements in the search string? For instance "American Eagle" is a two token `[brand name]`; "bathing suit" is a two token `[item]`; and "periwinkle blue" is a two token `[color]`. So I would like for "periwinkle blue American Eagle bathing suit" to have only three tags, namely `[color] [brand name] [item]`. The approach that I'm following comes from chapter 6 of this really excellent book, but it just doesn't deal with multi-token symbols or equivalently multiple symbols emitted during a particular state.   * What are your ideas?   * Are there any papers that you can refer me to?


title: Sample size calculation for mixed models

text: I am wondering if there are any methods for calculating sample size in mixed models? I'm using `lmer` in R to fit the models (I have random slopes and intercepts).


title: Interpreting range bars in R's plot.stl?

text: I have trouble figuring out what the range bars in `plot.stl` exactly mean. I found Gavin's post on this question and read the documentation as well, I understand that they tell the relative magnitude of the decomposed components, but still I am not entirely sure how they work. E.g.: data: tiny bar, no scale seasonal: full bar, with scale ranging from -0.6 to 0.2 trend: another tiny bar (seems to be equal to data), no scale remainder: medium sized bar with scale from -1.5 to 0.5 I do not understand what's the basis to relation and why trend has no scale. I tried `stl` and `decompose` with identical results for multiplicative and additive methods.


title: Discussing binomial regression and modeling strategies

text: Today I have got a question about binomial/ logistic regression, its based on an analysis that a group in my department have done and were seeking comments upon. I made up the example below to protect their anonymity, but they were keen to see the responses. Firstly, the analysis began with a simple 1 or 0 binomial response (e.g. survival from one breeding season to the next) and the goal was to model this response as a function of some co-variates. However, multiple measurements of some co-variates were available for some individuals, but not for others. For example, imagine variable x is a measure of metabolic rate during labour and individuals vary in the number of offspring that they have (e.g. variable x was measured 3 times for individual A, but only once for individual B). This imbalance is not due to the sampling strategy of the researchers per se, but reflects the characteristics of the population they were sampling from; some individuals have more offspring than others. I should also point out that measuring the binomial 0\1 response between labour events was not possible because the interval between these events was quite short. Again, imagine the species in questions has a short breeding season, but can give birth to more than one offspring during the season. The researchers chose to run a model in which they used the mean of variable x as one covariate and the number of offspring an individual gave birth to as another covariate. Now, I was not keen on this approach for a number of reasons 1)Taking the average of x means losing information in the within-individual variability of x. 2) The mean is itself a statistic, so by putting it in the model we end up doing statistics on statistics. 3) The number of offspring an individual had is in the model, but it is also used to calculate the mean of variable x, which I think could cause trouble. So, my question is how would people go about modeling this type of data? At the moment, I would probably run separate models for individuals that had one offspring, then for individuals that had two offspring etc. Also, I would not use the mean of variable x and just use the raw data for each birth, but I am not convinced this is much better either. Thank you for your time (PS: I apologize that its quite a long question, and I hope that the example is clear)


title: Some doubts about conditional expectation

text: Suppose I have two random variable $X_t\sim NID(0,1)$ and $Y_t\sim NID(0,4)$ and $Cov(X_t,Y_t)=2$ . Consider the random variable $Z_t = X_t + Y_t + Y_tX_t$. $E(Z_t) = E(X_t)+E(Y_t X_t)+E(Y_t) = 0 + 2 + 0 = 2$ because, $E(X_t)=0$ $E(Y_t)=0$ $E(X_tY_t) = E(X_t)\times E(Y_t)+cov(X,Y) = 0 + 0 + 2$ But if I want to obtain $E(Z_t|X_t)$ what it is the reasoning behind? I am not able to figure out. The same with variance: $Var(Z_t) = Var(X_t)+Var(Y_t X_t)+Var(Y_t) = 1 + 2 + 4 = 7$ but what about $Var(Z_t|X_t)$ ?


title: kolmogorov smirnov test for two sample of different size

text: I have two samples A and B. A contains 100 elements while B has 150 elements. I want to know whether Kolmogorov-Smirnov test can be performed to compare these two samples A and B of different size? Kindly reply. With Regards Puglu


title: Notation conventions for random variables and their distributions

text: I get confused on the proper notations of meanings, as well as the meanings of some notations relating to random variables and their distributions. Below, I will list things that I think are true, as well as things that I don't understand, and I would love input/corrections. I have labeled each point/question with a number for ease of reference. If it is not appropriate to list items in a single question like this, please let me know. I thought it would be ok since they are all short.   1. A random variable is notated by a capital letter, e.g. $X$.   2. What does an operation on a random variable mean? (e.g., how do you interpret $X^2$ in words?).   3. A specific draw from a random variable is notated by either the lowercase letter (e.g. $x$) or the lowercase letter with a subscript (e.g. $x_1$) or an uppercase number with a number(e.g. $X_1$).   4. The random variable that is the $kth$ order statistic of $n$ draws from a random variable $X$ is notated as $X_{kn}$.   5. Is there a shorthand way to write "X is the random variable that is distributed by F(x) (or "cdf F(x)" or "B(a,b)" or any way to characterize a distribution)"?   6. Can I write $\mathbb{E}F(x)$ to mean the expectation of the variable distributed according to $F(x)$?   7. If I perform an operation on a variable X's cdf, for example, $F_{new}(x) = F_{old}(x)^2$ to get the cdf of the maximum of 2 draws from $X$, can I notate that in terms of $X$ somehow?   8. Is the appropriate way to write $(F(x))^2$ succinctly $F^2(x)$ or $F(x)^2$?   9. Is there any notational difference between a discrete and a continuous variable?


title: How can I test heteroskedasticity in a Tobit model with Stata 12?

text: 1. I want to test heteroskedasticity in a Tobit model with Stata 12. But I don't know how to do that.   2. When I used an OLS model, I tested heteroskedasticity and autocorrelation, and didn't find much, but the OLS model isn't significant (Prob > _F_ = 0.54) and _R_ ² is small = 0.01. How can I improve that problem? Thanks.


title: Consequences of non normal-distribution multiple regression

text: I am currently writing my thesis using multiple regression. However my data does not meet the regression assumption of normal distribution. I want to describe that, due to the non normal distribution, the interpretation of the data is limited. I don't want to transform the data, i just want to say what the impact of the non normal distribution on my regression results (N =110) is. If have been googling a lot, however only found the assumptions, but never really consequences when i still conduct my tests. Does any of you know, what the consequences of non normal distribution for my regression results are? And also very important, which paper could be cited? I would be very grateful for your hint!


title: Acceptance rates for Metropolis-Hastings with uniform candidate distribution

text: When running the Metropolis-Hastings algorithm with uniform candidate distributions, what is the rationale of having acceptance rates around 20%? My thinking is: once the true (or close to true) parameter values are discovered, then no new set of candidate parameter values from the same uniform interval would increase the value of the likelihood function. Therefore, the more iterations I run, the lower the acceptance rates I should get. Where am I wrong in this thinking? Many thanks! Here is the illustration of my calculations: $$Acceptance\\_rate = \exp \\{l(\theta_c|y) + \log(p(\theta_c)) - [l(\theta^*|y) + \log(p(\theta^*) ]\\},$$ where $l$ is the log-likelihood. As $\theta$ candidates are always taken from the same uniform interval, $$p(\theta_c) = p(\theta^*).$$ Therefore acceptance rate calculation shrinks down to: $$Acceptance\\_rate = \exp \\{l(\theta_c | y) - [l(\theta^* | y) ]\\}$$ The acceptance rule of $\theta_c$ is then as follows: If $U \le Acceptance\\_rate $, where $U$ is draw from uniform distribution in interval $[0,1]$, then $$\theta^* = \theta_c,$$ else draw $\theta_c$ from uniform distribution in interval $[\theta_{min}, \theta_{max}]$


title: Is my data fit "extreme value distribution" or "normal distribution"?

text: I have a large data.frame in R. I would like to double if its distribution fit normal distribution or extreme value distribution better Here is my simplified data.frame.               x <- data.frame(A=c(1,3,1,5,4,5,5,7,3,2,2,1,1,1,4,9,10))      Could you mind to let me know how to do so? Can I did this analysis with R?


title: Multivariate orthogonal polynomial regression?

text: As a means of motivating the question, consider a regresison problem where we seek to estimate $Y$ using observed variables $\\{ a, b \\}$ When doing multivariate polynomial regresison, I try to find the optimal paramitization of the function $$f(y)=c_{1}a+c_{2}b+c_{3}a^{2}+c_{4}ab+c_{5}b^{2}+\cdots$$ which best fit the data in a least squared sense. The problem with this, however, is that the parameters $c_i$ are not independent. Is there a way to do the regression on a different set of "basis" vectors which are orthogonal? Doing this has many obvious advantages 1) the coefficients are no longer correlated. 2) the values of the $c_i$'s themselves no longer depend on the degree of the coefficients. 3) This also has the computational advantage of being able to drop the higher order terms for a coarser but still accurate approximation to the data. This is easily achieved in the single variable case using orthogonal polynomials, using a well studied set such as the Chebyshev Polynomials. It's not obvious however (to me anyway) how to generalize this! It occured to me that I could mutiply chebyshev polynomials pairwise, but I'm not sure if that is the mathematically correct thing to do. Your help is appreciated


title: Expected required sample length to train a hidden Markov model

text: Say one wishes to train a hidden Markov model with $n$ hidden states, and (accidentally) the problem itself can be described with a hidden Markov model with $n$ (or less states). What is the expected number of samples one needs to learn such model correctly (with an error $<\epsilon$)?


title: Continue from my previous question - distribution for a set of data using R results

text: Follow the very useful answers from Peter Flom, Wayne and many others. I have now started using R and it gives me a feeling of python :) The results are below but I am not sure how should I go from here ? The density certain looks much better after log transformation. Can you please shed some light on how to do further analysis ? Thanks a lot. R - Results below: plot (density (messages$length)) ![enter image description here](http://i.stack.imgur.com/4YlXo.png) plot (density (log (messages$length))) ![enter image description here](http://i.stack.imgur.com/3jZDL.png) summary (messages)               > summary(message$mb)         Min.  1st Qu.   Median     Mean  3rd Qu.     Max.       0.00665  0.32610  0.88450  2.08500  2.35000 49.13000       qqnorm (messages$length) ![enter image description here](http://i.stack.imgur.com/y3R43.png) ===================================================================== EDIT: Thanks all all the answering ! I have tried the qqnorm with log(x) and it looks like a straight line ! Does this mean my data is pretty much following a Log-normal distribution ? qqnorm (log(messages$length)): ![enter image description here](http://i.stack.imgur.com/TK5CA.png) Also I have tried to fit my data with a log-normal and below is the result. > fitdistr(message$mb, densfun="log-normal") meanlog sdlog   >  -0.19019347 1.45795269 ( 0.02003787) ( 0.01416891) Does this mean anything ?


title: How to explain ANOVA output from a multiple regression?

text: If I do a multiple regression such as:               df<-data.frame(y1=rnorm(100,2,3),     y2=rnorm(100,3,2),     x1=rbinom(100,1,0.5),     x2=rnorm(100,100,10))          fit<-lm(cbind(y1,y2)~x1+x2,data=df)     > anova(fit)     Analysis of Variance Table                      Df  Pillai approx F num Df den Df Pr(>F)         (Intercept)  1 0.75423  147.306      2     96 <2e-16 ***     x1           1 0.00720    0.348      2     96 0.7069         x2           1 0.00928    0.450      2     96 0.6391         Residuals   97                                               ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      I am wondering how to explain this ANOVA object where two models have different responses and the same set of predictors.


title: parameter estimation in bivariate autoregressive sem with individual heterogeneity - troublesome likelihood surface

text: I am working with models for panel data of individuals, wherein:   my latent variables Yt= $\beta$Yt-1+Ai+Q   Individuals level of the process A~N(l,h)   Innovation of the latent Q~N(0,q)   observed variables Xt=$\lambda$Y+$\epsilon$   measurement residuals $\epsilon$~N(0,e) Basically, a multivariate vector autoregressive model with parameters fixed over time, where the individuals can differ in mean level of the process, and where there is spurious guasssian noise in the observed variable. The model is not critical, but I thought it best specify. I am using the sem package openmx within R, which relies on the NPSOL optimizer for FIML estimation. I am concerned because even though my initial estimation appears to converge and generate reasonable estimates, when I vary my starting parameters or add upper and lower parameter bounds to the model, I am able to arrive at a better solution, with lower -2ll. That is, my likelihood surface appears difficult to navigate. What are the best approaches in such a case - my current line of thinking is to iterate over all viable parameter estimates rather than try to optimize, or some combined strategy where I perhaps iterate over some of the parameters while optimizing the others. Any suggestions or reading material on the topic would be welcome.


title: PLA vs Regression

text: Just started getting into machine-learning, and I'm wondering if there is a relationship between the Perceptron learning algorithm and linear regression?


title: Given samples from multiple normal RVs, how do we recover the histogram of their means?

text: Let $X_1,...,X_N$ be independent normal random variables. $X_i$ is normal with mean $\mu_i$ and standard deviation $\sigma_i$. Let $x_i$ be a single random sample from $X_i$. **Input:** We get all $x_i$'s and all $\sigma_i$'s, but we don't get the $\mu_i$'s. **Question 1** Estimate the histogram of the $\mu_i$'s. **Question 2** Assume the means $\mu_i$ are independently drawn from a distribution $\mathcal{M}$. Generate an estimate of $\mathcal{M}$. Note: Question 1 was answered below by @soakley but the solution didn't help my application, so I added Question 2. Note that the goal is not to estimate the $\mu_i$'s individually, but rather to get a good estimate for the histogram, or the distribution, of all the $\mu_i$'s together. This estimate should hopefully be better than what we'd get by simply mixing the Gaussians around the $x_i$'s. A maximum-likelihood estimate should work, but I don't know how to produce it. **Warm up question:** A simpler question is the above when all the $\sigma_i$'s are the same. This is easy: see answer at the bottom. **More details:** I need a method that I can program and will run in reasonable time. So exponential-time algorithms will not be sufficient. In my input, $N$ is around 5000, the standard deviations are mostly between 5 and 50, and the means are mostly between 0 and 40. So far I've tried the naive solution of drawing a Gaussian around each $x_i$ and mixing all these Gaussians. The results don't look at all like the correct distribution of $\mu_i$'s. For example, imagine all $X_i$'s are standard normal RVs. Then my naive method would guess that the $\mu_i$'s lie in a pretty wide Gaussian around 0. However, given the samples $x_i$ and the standard deviations $\sigma_i$, a clever algorithm could clearly see that the best guess is that all $\mu_i$'s are equal to zero. Therefore, it should be possible to do much better than my naive algorithm, possibly with some clever use of Fourier Decomposition. **Answer to warm-up question:** Question 1 doesn't have a good solution; to get a good solution, one would need to assume a posterior distribution. As for question 2: when all standard deviations are the same, then to get an estimate for the distribution $\mathcal{M}$ we simply need to de-convolve the distribution of $x_i$'s with a Gaussian. I don't see how to generalize this to the case of differing $\sigma_i$'s, though. **Motivation:** I am a poker player playing a very swingy poker game (Pot Limit Omaha). We want to find out if the rake is too high, by figuring out the "true winrates" of the player pool. We have as data the winrates of all players in the player pool over a whole year (these are $x_i$), and the standard deviations of their winnings (this is $\sigma_i$) and we want to estimate the distribution of their "true winrates" (the $\mu_i$'s) in order to figure out if enough players are winning. This translates to the above problem. **Ongoing Research** I just found a paper of Bovy et al which seems to address a generalization of my question and suggests an algorithm. It seems pretty relevant. I'll read it and report any findings here.


title: Using dummies in random effects model

text: When using dummy variables in a random effects model, how do I specify the model? Assuming three dummies (e.g. for three industries):   * Do I treat it as three regressions, where I present a table with three columns for each of these dummies? (where all other variables are the same below in each column)   * Or do I have a table with only one column but the three dummy variables included below in the column?  EDIT: Particularly, what is the difference between the two ways of presenting the tables? Does it make sense to present first the table inlcuding all variables in one column, and present in a second step subsample analyses considering each dummy-category? Thanks a lot for any prompt answer bob


title: Estimate multinomial probit model with mlogit (R package)

text: From the document and help, probit model is supported by mlogit. But when I tried it with these R scripts, the estimation takes much longer time to run (than the logit verion) and the result is quite a bit different as well (argument probit=FALSE). Does the probit behave correctly? If so, how should I interpret the er.gc, er.gr, etc coefficients?               > require(mlogit)     > data(Heating)     > H <- mlogit.data(Heating, shape="wide", choice="depvar", varying=c(3:12))     > m1.probit = mlogit(depvar~ic+oc, H, probit=TRUE)     > summary(m1.probit)          Call:     mlogit(formula = depvar ~ ic + oc, data = H, probit = TRUE)          Frequencies of alternatives:           ec       er       gc       gr       hp      0.071111 0.093333 0.636667 0.143333 0.055556           bfgs method     37 iterations, 0h:4m:54s      g'(-H)^-1g = 0.011      last step couldn't find higher value           Coefficients :                       Estimate  Std. Error t-value Pr(>|t|)       er:(intercept)  2.5611e-01  3.6641e-01  0.6990  0.48457       gc:(intercept) -2.6944e-02  3.3211e-01 -0.0811  0.93534       gr:(intercept) -1.8439e+01  3.2798e+01 -0.5622  0.57398       hp:(intercept) -6.4231e-01  7.4214e-01 -0.8655  0.38677       ic             -1.1447e-03  5.3175e-04 -2.1528  0.03133 *     oc             -3.3779e-03  1.4011e-03 -2.4109  0.01591 *     er.gc           4.4987e-01  2.6880e-01  1.6736  0.09421 .     er.gr           5.8580e+00  1.1236e+01  0.5214  0.60212       er.hp           1.2613e+00  5.0231e-01  2.5109  0.01204 *     gc.gc           7.1013e-01  3.5489e-01  2.0010  0.04540 *     gc.gr          -8.4606e+00  1.6848e+01 -0.5022  0.61555       gc.hp           6.7245e-01  6.1475e-01  1.0939  0.27401       gr.gr           1.4085e+01  2.6034e+01  0.5410  0.58849       gr.hp           4.9476e-01  4.4568e-01  1.1101  0.26694       hp.hp           2.2620e-01  2.9062e-01  0.7783  0.43637       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Log-Likelihood: -1000.1     McFadden R^2:  0.021626      Likelihood ratio test : chisq = 44.213 (p.value = 1.4049e-05)          > m1.logit = mlogit(depvar~ic+oc, H, probit=FALSE)     > summary(m1.logit)          Call:     mlogit(formula = depvar ~ ic + oc, data = H, probit = FALSE,          method = "nr", print.level = 0)          Frequencies of alternatives:           ec       er       gc       gr       hp      0.071111 0.093333 0.636667 0.143333 0.055556           nr method     6 iterations, 0h:0m:0s      g'(-H)^-1g = 9.58E-06      successive function values within tolerance limits           Coefficients :                       Estimate  Std. Error t-value  Pr(>|t|)         er:(intercept)  0.19459102  0.20424212  0.9527 0.3407184         gc:(intercept)  0.05213336  0.46598878  0.1119 0.9109210         gr:(intercept) -1.35058266  0.50715442 -2.6631 0.0077434 **      hp:(intercept) -1.65884594  0.44841936 -3.6993 0.0002162 ***     ic             -0.00153315  0.00062086 -2.4694 0.0135333 *       oc             -0.00699637  0.00155408 -4.5019 6.734e-06 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Log-Likelihood: -1008.2     McFadden R^2:  0.013691      Likelihood ratio test : chisq = 27.99 (p.value = 8.3572e-07)


title: Difference of random variables subject to non-negative values

text: This is probably close to trivial for stats experts. Yet, although having searched the web for quite a while now, I have not been able to come up with a satisfactory answer. Given two (approximately) Gaussian random variables, $u$ and $v$, what is the PDF and the expectation value of $u – v$ given the prior knowledge that $u – v \ge 0$? I remember vaguely having seen something along this lines based on Bayesian inference. However, now I seem unable to locate this resource.


title: The bootstrapped confidence intervals for ORs are strangely big or small in my binary logistic regression. Are they valid in the first place?

text: This question is an extension to my previous question, although it is not repeated. My bootstrapped ORs are OK. They are equal to non-bootstrapped ORs. However, the confidence intervals of the bootstrapped and non-bootstrapped regressions differ vastly. While the non-bootstrapped regression reports for example OR's 95% CI ~= 1.2 to 2.3, the bootstrapped one reports OR's 95% CI ~= 1.1 to 25000. It is **strange** and **non-practical**. My questions are:   1. Is this bootstrapped OR CI **correct**?   2. Can I keep my bootstrapped model in my study for its advantages, but report the non-bootstrapped OR CI?   3. Can I report both regressions? Many thanks.


title: Multidimensional scaling pseudo-code

text: I am planning to write a program that performs MDS. Any pointers to where I can access the pseudo-code for MDS? Thanks!


title: How to explain how I divided a bimodal distribution based on kernel density estimation

text: I have a dataset of bimodal population. It contains a smaller peak, which is considered to be "bad", and a bigger peak. I try to separate the bad part of data from the rest of data. What I did was: first I did a kernel density estimation, then found the local maximum of this small peak, and the local minimum of the pit between two peaks, then I took the midpoint (arithmetic mean of x-coordinates) of them, and define it as a cutoff. Everything below this cutoff is considered "bad". The reason I took midpoint instead of pit is because I tried to be more conservative. Now I would like to ask: Is what I did reasonable? If yes, how can I explain my action in a statisticians-favored manner? If not, how can I change? (Any other methods are welcome, especially those implemented in R.) Thank you! Here's the figure. ![enter image description here](http://i.stack.imgur.com/C0KZv.png)


title: Number of features vs. number of observations

text: Are there any papers/books/ideas about the relationship between the number of features and the number of observations one needs to have to train a "robust" classifier? For example, assume I have 1000 features and 10 observations from two classes as a training set, and 10 other observations as a testing set. I train some classifier X and it gives me 90% sensitivity and 90% specificity on the testing set. Let's say I am happy with this accuracy and based on that I can say it is a good classifier. On the other hand, I've approximated a function of 1000 variables using 10 points only, which may seem to be not very... robust?


title: Data exploration tools for Excel

text: I have a sample of records, with string and numerical columns. My sample is currently hosted in an excel spreadsheet. I need a tool that wil produce discripte statistics for each colomn, such as max and min values, number of unique values, max string lenght. Are there any such tools available. Any tips? The ultimate goal is to derive a data dictionary.


title: Averaged estimators in stochastic versions of EM

text: Recently I've been working EM algorithms for MAP estimation in a problem where the expectation is intractable, but the maximization is easy. Further, draws from the distribution in the E-step are easily available through MCMC, so I've been experimenting with stochastic versions of EM. Let $X$ be the observed data, $Z$ be the missing data, and $\Theta$ be the parameters to be estimated. I'll use $\Theta_t$ for the current estimate of $\Theta$. Specifically I've looked at stochastic approximation EM and Monte Carlo EM, which estimate the E step like so: 1) MCEM: The $Q$ function is approximated with an average of $m_t$ MC draws $\hat Q(\Theta; \Theta_t) = \frac{1}{m_t}\sum_{i=1}^{m_t}\log(p(X,Z_i|\Theta)$ where $Z_i$ is drawn from $p(Z|X,\Theta_t)$ 2) SAEM: $\hat Q(\Theta; \Theta_t) = \gamma_t \hat Q(\Theta; \Theta_{t_1}) + (1-\gamma_t)\frac{1}{m_t}\sum_{i=1}^{m_t}\log(p(X,Z_i|\Theta)$. In the context of exponential families, this amounts to using a weighted average of the "new" and "old" sufficient statistics at each step, compared to MCEM which takes a completely new average (ie $\gamma_t=0$) Both have the same M step, with $\Theta_{t+1}$ chosen to maximize the approximate $Q$ function. My question is when does it make sense to average over the sequence of $\Theta_t$'s? Presumably this would further reduce the Monte Carlo error, and my understanding is that the convergence is (theoretically) unaffected since the $\Theta_t$'s are converging to a stationary point with appropriate conditions on $m_t$ and $\gamma_t$. What I'm looking for is more practical advice, since averaging too early will bias the (finite sample) estimates and averaging too late is wasteful. Also the related question of assessing convergence; if we consider convergence criteria like $|\Theta_t - \Theta_{t-1}| \leq \epsilon$ should the difference be taken instead between the two _averaged_ estimates? (Not the best criterion, I know, but super easy :) )


title: Review of case control matching algorithms?

text: Can you suggest a good review of case control matching algorithms? Algorithms that can be used to set up the matched _pairs of one case and one control_ , or _matched blocks of a case and multiple controls_.(Preferably a paper, book chapter or website discussing recent developments.)


title: How to compute the standard error of measurement (SEM) from a reliability estimate?

text: SPSS returns lower and upper bounds for Reliability. While calculating the Standard Error of Measurement, should we use the Lower and Upper bounds or continue using the Reliability estimate. I am using the formula : $$\text{SEM}\% =\left(\text{SD}\times\sqrt{1-R_1} \times 1/\text{mean}\right) × 100$$ where SD is the standard deviation, $R_1$ is the intraclass correlation for a single measure (one-way ICC).


I can see that there are a lot of formal differences between Kullback–Leibler vs Kolmogorov-Smirnov distance measures. However, both are used to measure the distance between distributions.   * Is there a typical situation where one should be used instead of the other?    * What is the rationale to do so?


title: Fitting an ARIMAX model with regularization or penalization (e.g. with the lasso, elastic net, or ridge regression)

text: I use the auto.arima() function in the forecast package to fit ARMAX models with a variety of covariates. However, I often have a large number of variables to select from and usually end up with a final model that works with a subset of them. I don't like ad-hoc techniques for variable selection because I am human and subject to bias, but cross-validating time series is hard, so I haven't found a good way to automatically try different subsets of my available variables, and am stuck tuning my models using my own best judgement. When I fit glm models, I can use the elastic net or the lasso for regularization and variable selection, via the glmnet package. Is there a existing toolkit in R for using the elastic net on ARMAX models, or am I going to have to roll my own? Is this even a good idea? edit: Would it make sense to manually calculate the AR and MA terms (say up to AR5 and MA5) and the use glmnet to fit the model? edit 2: It seems that the FitAR package gets me part, but not all, of the way there.


title: Spearman's or Pearson's correlation with Likert scales where linearity and homoscedasticity may be violated

text: I am wanting to run correlations on a number of measurements where Likert scales were used. Looking at the scatterplots it appears the assumptions of linearity and homoscedasticity may have been violated.   * Given that there appears to be some debate around ordinal level rating approximating interval level scaling should I play it safe and use Spearman's Rho rather than Pearson's r?    * Is there a reference that I can cite if I go with Spearman's Rho?


title: Proof of normal equation in regression using tensor notation

text: I'm struggling with a proof of the normal equation, so I posted a question which hopefully will get resolved soon. However, I mentioned there that I'm uncomfortable with the proofs dealing with matrix calculus, particularly when it comes to derivatives of product of vectors or product of matrices and vectors. Furthermore, there are a couple of "layout conventions" which are not even used consistently. Because of this, I thought I could get rid of all that insanity by using tensor notation instead. In relation to the normal equation question, I would be very interested in looking at a proof of this equation (here you can see the context): $$\beta = (X^{T}X)^{-1}X^{T}{\bf{y}}$$ using tensor notation. And as a more general request, I'd like to get an advice on the issue of dealing with matrix calculus and derivatives. Particularly, I want to hear about your approach when you have to calculate something of this kind given the inconsistencies in notations and non-intuitive nature of some identities. Do you consult Wikipedia or some handbook every time you want to calculate a derivative of some product of vectors or matrices? Do you memorize some recurrent identities? Do you know an alternative procedure?


title: Advice on explaining heterogeneity / heteroscedasticty

text: I am looking for any help, advice or tips in how to explain heterogeneity / heteroscedasticity to biologists in my department. In particular I want to explain why its important to look for it and deal with it if it exists, I was looking for opinions on the following questions.   1. Does heterogeneity influence the reliability of random effect estimates? I am pretty sure it does, but I couldn't find a paper.    2. How serious a problem is heterogeneity? I have found conflicting views on this, while some say that model standard errors etc. will be unreliable, I have also read that it is only a problem if the heterogeneity is severe. How severe is severe?    3. Advice on modelling heterogeneity. Currently, I focus largely on the nlme package in R and the use of variance covariates, this is pretty simple and most people here use R so providing scripts is useful. I am also using the MCMCglmm package as well, but other suggestions are welcome, particularly for non-normal data.    4. Any other suggestions are welcome.


title: Renewal Process Hypothesis Test

text: I have $n$ realisations $s_1,\, \dots , s_n$ of random variables $S_1,\, \dots, S_n$ which are assumed to be i.i.d. with unknown distribution. These measure the time between events. I want to calculate the probability that the data can be modelled by a renewal process: $(X_t)_{t\geq0}$ where $X_t$ is the total number of jumps by time $t$. In particular, I am using the elementary renewal equation $\lim_{t\rightarrow \infty} \frac{E[X_t]}{t} = \frac{1}{E[S_i]}$ and I want the probability that the observed values of the LHS and RHS agree for a given $(n, t)$. Do I want to some how construct a confidence interval for the statistic $\frac{E[X_t]}{t}$ and test whether $\frac{1}{E[S_i]}$ lies in this interval? I guess what I'm asking, is how can I calculate the error in the strong law of large numbers for a given $n$, I think? Apologies for the confusing description, my statistics isn't the best!


title: How to find parameters for ridge and lasso regularization when cost minimization does not converge?

text: In the Stanford ML course, we were taught to find good values for the lambda parameters of ridge/lasso by iterating for various lambda values on several cross-validation sets and picking the values which correspond to the hypothesis with the minimum CV error. The problem is, I am playing with a big data set (which might not even be well suited for logistic regression (??): Internet Ads set) and I can't use the method described above because, during optimization, the cost stays stable (changes only around the 8th decimal place between iterations) and doesn't seem to converge. I need good values for the regularization in order to converge, but I can't converge without the aforementioned good values. Any suggestions? Should I move on to using SVM, or is this data set solvable with logistic regression? NB: I am doing this for learning purposes, so I'm much more interested in explanations why my approach is bad than I am in black-box libraries which will give me a solution. **EDIT: Some relevant code snippets (occasionally pseudocode-ish, for clarity). The usual notations apply.** The function used to compute the cost:               def computeCost(theta, X, y):         global iter         iter += 1         if iter > 10:             raise TooManyIterationsException(iter) # Because the cost doesn't converge, I force interruption in order to jump to another combination of lambda values         m = y.size         h = sigmoid(X.dot(theta.T))         J = y.T.dot(log(h)) + (1.0 - y.T).dot(numpy.log(1.0 - h))         J_reg2 = theta[1:]**2         J_reg1 = theta[1:]         cost = (-1.0 / m) * (J.sum()) + LAMBDA2 * J_reg2.sum() + LAMBDA1 * J_reg1.sum()         print "Cost: ", cost         return cost      Invoking scipy.optimize.fmin_bfgs:               initial_thetas = numpy.zeros((len(train_X[0]), 1))     myargs = (train_X, train_y)          for LAMBDA1 in [0.01, 0.02, 0.04, ..., 10]:         for LAMBDA2 in my_range[0.01, 0.02, 0.04, ..., 10]:             try:                 iter = 0                 theta = scipy.optimize.fmin_bfgs(computeCost, x0=initial_thetas, args=myargs)                         except TooManyIterationsException as e:                 print '\n'      A typical output looks like this: ![enter image description here](http://i.stack.imgur.com/b9gdK.jpg) **EDITED AGAIN: Evolution of thetas!** ![enter image description here](http://i.stack.imgur.com/xzDCM.jpg)


title: Hypothesis testing in multiple regression

text: I have the following variables:   1. weight (ratio scale)   2. result obtained from psychological test (eg. iq) (ratio scale)   3. style1 (ratio scale)   4. style2 (ratio scale)   5. style3 (ratio scale) I want to test this hypothesis: People with high result on style 1 and high result on IQ scale are having bigger weight than people with high result on style 2 and low result on IQ scale. Should I make a categorical variable 'style' and recode values of style 1, style 2, and style 3, and then proceed regression analysis, or there is a way to test hipothesis without this conversion. I bet it should be multivariable regression analysis, but I don't know how to divide file in groups, when I have three variables.


title: Forecasting Extreme values

text: I am attempting to forecast "peak" values. I have a weekly data set ~10 years. I have weather metrics as explanatory variables. I am trying to produce a model that can accurately predict the peaks of the weekly data using the weather metrics. The peaks are not entirely accurate, should I be using another approach? I am currently using MLR. Pred vs actual shows good explanatory ability until you reach more extreme observations. Thanks


title: Assessing the consistency of an estimator of the average path length of a large network

text: I have found the following in a network-science article and would like clarification of what the authors are claiming: > We observe that the root mean square of the difference between the empirical > probability density distributions of sample networks with $h$ and $(h + 1)$ > seeds steadily decreases fitting to the $y = 1/x$ line, confirming that our > estimation would converge to the true average path length. The authors are trying to measure the average path length of a large network by sampling vertices at random (which they refer to as _seeds_ ; presumably the sampling is uniform) and then performing a breadth-first search from the selected vertex to determine the length of the path to every other node in the network. (Implicitly, the graph should be connected.) Can anybody guess what it is that the authors do here regarding their consistency claim? Is there a name for this method? And could you please explain it better so I can understand how to do the same thing in a similar situation? The excerpt above is from page 839, just above the section heading for Section 5.2, of the following article. > Y.-Y. Ahn, S. Han, H. Kwak, S. Moon and H. Jeong (2007), Analysis of > topological characteristics of huge online social networking services, > _Proc. WWW 2007 (Track: Semantic Web)_ , pp. 835-844.


title: Alternative argument in Fisher's exact test in R

text: Can anyone please tell me what the `alternative` argument is in `fisher.test()` in R? What do greater and less imply?               fisher.test(x, y = NULL, workspace = 200000, hybrid = FALSE,                 control = list(), or = 1, alternative = "two.sided",                 conf.int = TRUE, conf.level = 0.95,                 simulate.p.value = FALSE, B = 2000)      I tried the help section and all it says is that alternative hypothesis must be one of `"two.sided", "greater", or "less"`. When should I use each one? I'm looking at over-representation in my set against a control set.


title: How to customize axis labels in a boxplot?

text: How can I change the labels of the vertical y axis in a boxplot, e.g. from numbers to text? For example, I would like to replace {-2, -1, 0, 1, 2} with {0hour, 1hours, 2hours, ...}.