title: How to correctly measure effect on heavy-tailed distribution

text: I ran an experiment on my website where I randomly assigned users to either the treatment or control group, and have two questions about how to correctly compute significance of the results. Some context: The outcome variable is the total number of actions they took on my website. I am computing the mean number of actions taken by users in my control and treatment groups, respectively: $A_C = \frac{\sum_{i \in C} a_i }{|C|}$, $A_T = \frac{\sum_{i \in T} a_i}{|T|}$ (where $T$ and $C$ are my treatment and control groups and $a_i$ is the number of actions performed by user $i$). Is $A_T-A_C$ a reasonable outcome to measure? My concern is that since the $a_i$'s are distributed according to a heavy-tailed distribution, the mean is heavily skewed by a few power users. Usually I'd take the median, but the tails are so heavy that the median is very small. Another thing I have measured is the KS statistic on the distributions of number of actions for users in the control and treatment groups.   1. Is the difference in means or KS-statistic a more reasonable way to determine if the treatment had an effect?   2. I am calculating bootstrapped p values by computing the same measures (difference in means and KS statistic) for many random splits of the user population, then counting how many random splits give rise to "more significant" values than my experimental control/treatment split gave. Is this correct?


title: What GEE-exchangeable method can do that robust variance can't?

text: I asked a related question before here on the difference between GEE method with exchangeable varcov structure v. Robust standard errors known as Huber White method in group randomized trials. As Macro pointed out Freedman in his 2006 paper The American Statistician called On The So-Called “Huber Sandwich Estimator” and “Robust Standard Errors explains that with a little modification of the huber white method, we can get valid inference by both solving the issue with Heteroscedasticity and Correlation among clusters. I have a study with 34 schools randomly assigned to treament/control. The ICC in the study is 0.04 (although is small but is statistically significant, more info here) with the design effect of 17.7. The ICC indicates some degree of similarities among students within schools (i.e. clusters). Schools have different number of students varying from 150 to 800. My question is if you wanted to analyze the data, which method you would choose to get valid inference? Huber-White or GEE with exchangeable varcov matrix? I'm in favor of GEE because of two reasons: 1) we have unbalanced clusters 2) Our estimates with GEE are more efficient. If you choose GEE, could you explain why and what GEE can do that HuberWhite cannot do and the other way around? I appreciate your help.


title: How to average variables having given standard deviations

text: I'm a physicist who has some issues in data treatment, so sorry if I get the maths wrong. I have some diffractograms, which are displayed as counts vs. angle. The points have standard deviations estimated from a poisson statistic on the bins (an external program does that). Getting the average of different spectra is straightforward by averaging corresponding counts. But how do I get the correct standard deviation for the averaged counts? When I apply this method Is it possible to find the combined standard deviation? the errors actually get larger compared to single data points. This appears counterintuitive, because I would expect the errors to get smaller when I take into account more datapoints? Where is my fault?


title: How to handle unseen features in a Naive Bayes classifier?

text: I am writing a naive bayes classifier for a text classification problem. I have a bunch of words and an associated label:               [short,snippet,text], label1     [slightly,different,snippet,text], label2     ...      I am able to train the naive bayes fine. However, when I am classifying unseen data, sometimes there are unseen features (words). In that case, what happens to the naive bayes formula to determine the probability of a class $C$ given features $F_1,F_2,...$? $$P(C|F_1,F_2,...) = \frac{P(F_1,F_2,...|C)P(C)}{P(F_1,F_2,...)} = \frac{P(C)\prod_{i}P(F_i|C)}{P(F_1,F_2,...)}$$ Say feature $F_k$ never occured in the training data, then isn't $P(F_k|C)=\frac{0}{0}$? How is this typically handled in classification problems? One option is to simply ignore unseen features. However, I would not like to do that, since I am trying to calculate the actual probability score associated with classes. Probabilities should take a hit when there are unseen features, but I am not sure how to do that mathematically. Any insights, links to reseach articles, etc would really help! Thanks in advance.


title: How can I calculate covariance of quasi-Poisson-distributed variables in R?

text: I have two time series of count data which I believe are correlated. Is there an equivalent to `cov()` that assumes an overdispersed Poisson rather than a normal distribution? If there is a significance test for this (noting that both series increase over time but that doesn't prove correlation), so much the better. I wondered if `vcov()` might help but it focuses on parameters of a solved model. More broadly, I'm looking for a measure of the similarity between two variables which is appropriate for quasi-Poisson data.


title: How to determine if two sequences are significantly different?

text: I have two sequences of equal length, like `[1, 1, 2, 2, 2, 3, 3, 3, 4]` and `[1, 2, 2, 3, 3, 3, 3, 4, 4]`. I want to determine the odds that they were generated by the same random variable, or different ones. What is a good statistical test to do this? The data are not pairwise; nor are they time series. Each element in the sequence is independent of the others. (I'm sorry if I'm omitting necessary details; I'm new to stats.)


title: Test if course marks are better this semester compared to last semester

text: I am teaching a programming course at the university where I am working. I would like to do a statistical test to see if students marks (exam grades) are significantly better compared to the previous semester (e.g., summer semester vs. winter semester). There are three exams each semester, and the grades are scored from 1 (best mark) to 5 (fail). There are about 24 students that take the exam each time. What is the right test to use, in order to assess that the marks from the summer semester are significantly better than those from winter semester, considering 3 exams per semester? (I would like to run the test in R, so if someone happens to know how it's called in R that would be perfect. But I think I can find that out myself.) Background: I would like to apply for a teaching award. Although marks are probably not a good test of teaching performance if the course design improved, I would like to provide the statistical results as a indication of the merit of my nomination along with the student feedback from the instructor evaluation questionnaire. EDIT:   The students are different each semester (except possibly for those who fail multiple times). The instructor is different (I am the second instructor). The test is the same. Yes, the difference in the capabilities of the students is a potential confounding factor. I would like to run the test anyway and perhaps I can state that its the difference in the set of students in the classes is not a major factor affecting the difference in the distribution of sources and average scores.


title: PROC NLMIXED and PROC LIFEREG not arriving at the same answer for Log-normal survival function

text: Related question I have a project where, despite being able to implement some parametric models in LIFEREG, it is somewhat more convenient to do it in NLMIXED. Verifying that this technique works, I tried implementing a pretty bog standard Weibull AFT as follows:               data work.data;         call streaminit(123);         do i = 1 to 1000;             x = rand('BERN', 0.10); *10% exposed;             t1 = rand('Weibull', 1, 20*(exp(-0.500*x)));             t2 = t1;             weight = 1;             output;         end;     run;          PROC NLMIXED data=work.data fd;         parms alpha=1 f0=-1 f1=0;         bounds alpha>0;         *Rate of x;         lam=exp(-(f0*alpha+f1*alpha*x));                  *Density of x;         ff1=alpha*lam*t1**(alpha-1)*exp(-lam*t1**alpha);              *log Likelihood;         logl=log(ff1);              *Weighted log Likelihood;         wlogl=logl*weight;         model t1~general(wlogl);         ods exclude iterhistory parameters;     run; quit; run;          proc lifereg data=work.data;             model (t1, t2) = x / D=WEIBULL;          weight weight; run;      This works swimmingly. Both LIFEREG and NLMIXED produce the same estimates, and all is right in the world. Testing it for a Log-Normal AFT model however using the following modifications doesn't work so well:               PROC NLMIXED data=work.data fd;         parms sigma=1 g0=-1 g1=0;         bounds sigma>0;              *Rate of X;         mu=exp(g0+g1*x);              *Density of X;         fg1 = exp(-0.5*((log(t1)-mu)/sigma)**2)/((t1*(2*CONSTANT('PI'))**0.5)*sigma);              *log Likelihood;         logl=log(fg1);              *Weighted log Likelihood;         wlogl=logl*weight;         estimate "RT NDD" exp(g1);         model x~general(wlogl);         ods exclude iterhistory parameters;     run;          proc lifereg data=work.data;         model (t1, t2) = x / D=Lnormal;          weight weight;run;      Here, while both models estimate the scale as the same number (1.2676), the estimates for the intercept (g0) and x (g1) are wildly off. In NLMIXED they are 0.8747 and -0.1419 respectively, while in LIFEREG they are 2.3981 and -0.3172. Any idea what's going on?


title: How can you relate the coefficients of a multivariate regression of ${\bf Y} \sim {\bf X}$ to the coefficients of ${\bf X} \sim {\bf Y}$?

text: Is anyone aware of a orthogonal multiple regression library that is implemented in say R, Scipy, Matlab, Octave, etc.? (Or even fortran/C...) If I'm not mistaken, it would not be difficult to write but just wanted to check. Second question: my actual application is for a multivariate regression case where I have two matrices, $\mathbf{Y}$ and $\mathbf{X}$. If I am able to find a total least squares estimator for $\boldsymbol\beta$ and $\boldsymbol\beta'$, would the inverse of its transposed elements equal one another ($\hat{\beta}_{ji}^{-1} = \hat{\beta}'_{ij}$) if the two coefficient matrices are defined by $\mathbf{Y} = \hat{\boldsymbol{\beta}} \mathbf{X}$ and $\mathbf{X} = \hat{\boldsymbol{\beta}}' \mathbf{Y}$? It seems that this equality ($\hat{\beta}_{ji}^{-1} = \hat{\beta}'_{ij}$) does not necessarily hold as it does for the classic case, where the scalar values of $\beta^{-1}=\beta'$ when $\mathbf{Y}$ and $\mathbf{X}$ are simply vectors rather than matrices, but also wanted to check if this was a known fact. Thanks.


title: How to model a linear regression based on time?

text: I have some training set data variables $x_1$, $x_2$, $x_3$, $x_4$, and $x_5$ and a response variable $y$. But these are time series data. So the for the same set of values of $x_1$, $x_2$, $x_3$, $x_4$, and $x_5$ the response variable $y$ may be different in different observations. Can anyone please suggest me how to model my linear regression in this case? I am learning newly statistics. Please explain in the simplest way possible. I am uploading a picture explaining the same question for more clarity ![enter image description here](http://i.stack.imgur.com/DnKXs.png)


title: Applying log-transformation when comparing two populations

text: I have a survey data containing two samples drawing from two populations (or the same population but at different time series). Based on these two samples, I would like to estimate the total and its confidence intervals for each population. Since my sample distributions are highly skewed and both sample sizes are small. The recommended method is to first apply log-transform to the original data and compute CI in the log domain. However, since log-transform is nonlinear, i noticed that in the original domain sum(sample_1)>sum(sample_2), while in the log domain sum(log(sample_1)) is smaller than sum(log(sample_2)). Can someone help me with this situation? Shall I still use the data in log domain. If the purpose of the study is to compare the total, or even further to use the data in prediction, using the original and log domain would give different results? Thanks!


title: Drawing conclusion from fixed significance level or p-value in a two-sample test

text: This was an example done in class, However I was sick An experiment was performed to determine whether the average nicotine content of brand A cigarette exceeds that of brand B cigarette by 0.20 milligram. If 50 cigarettes of brand A had a sample mean of 2.61 milligrams whereas 40 brand B cigarettes had an average nicotine of 2.38 milligrams. The population standard deviations of the nicotine content for the two brands of cigarettes are known to be 0.12 and 0.14 for brand A and B, respectively. (a) Based on a significance level of 5%, what can you conclude about the difference between the two brands of cigarettes? (b) Base on a p−value, what can you conclude about the difference between the two brands of cigarettes? My Attempt: (a) $H_{0} :\mu_{A}-\mu_{B} =0.2$ $H_{1} :\mu_{A}-\mu_{B} \ne 0.2$ Significance Level : $\alpha = 0.05$ Rejection Region : $|z| >1.96$ Test Statistic : $ z = \frac{2.61-2.38 -0.2}{\sqrt{\frac{0.12^2}{50}+\frac{0.14^2}{40}}} =1.08$ Conclusion : Since $ 1.08 <1.96 $ I fail to reject $H_{0}$ at 5% I really need Help with B


title: Techniques for analyzing ratios

text: I am looking for advice and comments that deal with the analysis of ratios and rates. In the field in which I work analysis of ratios in particular is widespread but I have read a few papers that suggest this can be problematic, I am thinking of: > Kronmal, Richard A. 1993. Spurious correlation and the fallacy of the ratio > standard revisited. _Journal of the Royal Statistical Society Series A_ > 156(3): 379-392 and related papers. From what I have read so far it seems ratios can generate spurious correlations, force regression lines through the origin (which is not always appropriate) and modelling them may violate the principle of marginality if not done correctly (Using Ratios in Regression, by Richard Goldstein). However, there must be occasions when the use of ratios is justified and I wanted some opinions from statisticians on this topic.


title: Automatically choose a bandwidth in R

text: I have a function that I have added noise to. I want to test a few smoothers, but for each smoother I want to minimise MSE. I know the underlying function, so how can I make R automatically choose the bandwidth that minimises the MSE?


title: How to derive mean of chi square variable function using the MGF?

text: I am working through past examination questions from the Royal Statistical Society and came across this one from 2009 in Module 5 (Question 2(i) and Solution): > The random variable $X$ has a $\chi^{2}_{k}$ distribution ($k=1,2,3, ...$) > which has the moment generating function (mgf) $m(t) = (1 − 2t)^{−k/2}$ for > $t < \frac{1}{2}$. > > Using the mgf, find the mean and variance of $X$. I know the mean will be $k$ and variance $2k$ but I can't derive it. Here is my best attempt which is incorrect: $$ E[X] = d/dx[M_X(0)] = -2 \times \frac{-k}{2}(1-2X)^{-k/2 -1} = k^{\frac{-k}{2} -1} $$ I'm pretty sure I should be doing an expansion at $0$ but I can't see how to do it.


title: Cross-validation of panel data

text: When I perform k-fold cross validation, the observations are selected randomly into each fold. This works fine for cross-sectional data. However, when I have panel data (and fit with hierarchical model), does it make sense if the training set has observations that happen before those in the test set?


title: Minimum sample size for kernel density estimation

text: I was wondering if there is any recommendation for a "decent" minimum sample size when estimating the density/distribution/quantiles of a random variable using KDE. Is it always the case that "the more data the better"?


title: Time Series on Oil Filter Pressure

text: I am not really strong with time series but I have a project I am working on.. I have a problem where I am trying to model a time series of the difference in pressure before and after oil has passed through an oil filter. I am trying to see if this oil filter can go longer without being changed. I also have data on the temperature of this oil. The temperature of the oil is a strong indicator of Difference in Pressure, 'PD'. I was thinking that I would like to remove the effect of the Temperature by modeling the residuals from a model with Temperature as a predictor and PD as the response. Would this work? Or is there a better approach?


title: Cohort analysis for media site

text: I'm trying to find a reasonable way to measure **retention** of a site (news site). I'm trying to do cohort analysis: I've grouped all visits by year/week, and measured the percent of users that kept coming back in the following weeks. **Clarification** : these users are first-time visitors, who have been tracked along the weeks following their first time visit. I've got the following data:               yyyy/ww 1       2       3       4       5       6       7       8       9       10      11      12      13     201348  27.33   7.51    5.71    6.61    3       5.41    4.8     4.5     4.2     5.41    4.2     3.9     2.1     201349  7.88    5.23    3.03    2.53    4.02    3.16    3.13    2.83    2.63    2.16    2.06    1.25         201350  7.91    3.79    3.01    4.59    3.57    3.48    3.08    2.93    2.45    2.25    1.4          201351  6.64    4.1     5.69    4.41    4.13    3.68    3.36    2.8     2.56    1.56                 201352  6.09    6.99    5.05    4.79    3.99    3.79    3.24    2.94    1.84                     201401  11.66   6.69    6.28    5.45    4.67    3.74    3.4     1.99                         201402  8.45    6.14    5.11    4.61    3.49    3.1     1.77                             201403  9.45    6.27    5.38    4.12    3.57    2.08                                 201404  9.78    6.78    4.9     4.18    2.29                                     201405  8.96    5.37    4.45    2.47                                         201406  8.88    5.85    3.07                                             201407  8.46    3.76                                                 201408  6.29                                                           The row labels are a week in time, the column labels are the number of weeks that passed since the first visit. The cell values are the **percent** of inital users that came back. The two funky values at`week+1` account for christmas and new year. My question is: how should I interpret these values? I don't have a very deep statisical knowledge beyond the basics; I tried to calculate the slope of each row, and I got this (oldest to newest):               -4.956     -1.042     -0.788     -0.471     -0.64     -1.522     -1.145     -1.391     -1.758     -2.039     -2.905     -4.7      My rationale is that a good performing site would have an `increasing`slope value over time. In this case it seems the opposite. **Clarification** : What I mean with increasing slope is that the slope value would be less farther away from zero for the next first-time visitors; over time, the percentage of people who kept coming back would be higher (but still < 100% of course) than the previous cohorts. Something like this, where you can see the fraction of retained users is bigger than the previous month: ![enter image description here](http://i.stack.imgur.com/s33qF.jpg) Any suggestions and/or critiques are welcome. Thanks!


title: Truncated SVD: how do I go from [Uk, Sk, Vk'] to low-dimension matrix?

text: I have a large word-frequency matrix (~6m unique words X ~4k documents) and I'm trying to use truncated singular value decomposition (SVD) to project it onto a matrix with fewer dimensions. I know how to get to $[U_{k}, S_{k}, V_{k}']$, but I don't know what to do then. Each tutorial I find gives me a different answer. Here, for instance, I'm told to retain only the first $k$ rows of $S_{k}$ and $V_{k}'$ and then multiply $S_{k}V_{k}'$. In the scikit-learn implementation, however, they do it by multiplying $U_{k}S_{k}'$ (formula and code). (The book chapter they mention as source doesn't say anything about that.) In this previous answer they suggest $X'U_{k}$ instead (though they are talking about PCA, so maybe that doesn't apply here). Other answers haven't helped (here, for instance, they say we need to apply truncated SVD to the reduced matrix - which doesn't sound right, we use truncated SVD precisely to get to the reduced matrix, no?). So, bottom line: how do I go from $[U_{k}, S_{k}, V_{k}']$ to a reduced data matrix? Which of the above formulas is the correct one (and why)? Just to be clear: I'm not looking to reduce rank only (that I can do by multiplying $U_{k}S_{k}V_{k}'$ but it gives me the exact same dimensions I had before, which doesn't help me); I want to reduce the actual dimensions of the data matrix.


title: Need reasons or references on small p-values with large data sets

text: What are reasons or references for a statement such as "When you have a lot of data the statistical problem you run into is that even a tiny difference will be statistically significant." I see this quite often and I need references as to why this is true. Any help is appreciated.


title: Standard error of the mean in a repeated measures study

text: I am having difficulty figuring out how to calculate the SEM and report it in a table. I have 5 different time periods, and for those time periods 6 different measures of digestibility for each animal:               horse   period   adfdig'y   ndfdig'y     1         1         45.1      70.2                 2         1         48.2      69.8              3         1         55.6      59.0           4         1         49.1      57.0     5         1         42.0      66.3               6         1         46.0      60.2     1         2         56.3      70.6     2         2           etc, etc,       So the issue I'm having is building an ANOVA table like this with the SEM:               item, unit            treatment (time period)                   p1   p2  p3   p4    p5   SEM     ADF, %                           NDF, %                                           How do I get one SEM for ADF digestibility over the time periods?


title: The most optimal way to solve this set of non-linear equations in high dimensions

text: So I have a series of non-linear equations which I wish to solve as fast as possible, to illustrate for the case of $n = 4$, I have the following equations: \begin{gather*} c_0\exp(\lambda_0-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = 1 \\\ c_1\exp(\lambda_0+\lambda_1-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_1 \\\ c_2\exp(\lambda_0+\lambda_2-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_2 \\\ c_3\exp(\lambda_0+\lambda_3-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_3 \\\ c_4\exp(\lambda_0+\lambda_4-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_4 \end{gather*} where \begin{gather*} f = \int_{d_4}^{\infty} \int_{d_3}^{\infty} \int_{d_2}^{\infty} \int_{d_1}^{\infty} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\\ c_0 = \int_{-\infty}^{d_4} \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\\ c_1 = \int_{-\infty}^{d_4} \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{d_1}^{\infty} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\\ c_2 = \int_{-\infty}^{d_4} \int_{-\infty}^{d_3}\int_{d_2}^{\infty} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\\ c_3 = \int_{-\infty}^{d_4} \int_{d_3}^{\infty}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\\ c_4 = \int_{d_4}^{\infty} \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\\ \end{gather*} The definitions are as follows:   * $P_1, P_2, P_3, P_4$ are known values.   * $d_1, d_2, d_3, d_4$ are known values.   * $q(x_1, x_2, x_3, x_4)$ is a multivariate normal distribution with mean vector $0$ and a known variance-covariance matrix.   * Thus $f, c_0, c_1, c_2, c_3, c_4$ are all simply real numbers. * * * For the general case I have: \begin{gather*} c_0\exp(\lambda_0-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = 1 \\\ c_1\exp(\lambda_0+\lambda_1-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_1 \\\ c_2\exp(\lambda_0+\lambda_2-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_2 \\\ c_3\exp(\lambda_0+\lambda_3-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_3 \\\ \vdots \\\ c_n\exp(\lambda_0+\lambda_n-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_n \end{gather*} where \begin{gather*} f = \int_{d_n}^{\infty} \cdots \int_{d_3}^{\infty} \int_{d_2}^{\infty} \int_{d_1}^{\infty} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\\ c_0 = \int_{-\infty}^{d_n} \cdots \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\\ c_1 = \int_{-\infty}^{d_n} \cdots \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{d_1}^{\infty} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\\ c_2 = \int_{-\infty}^{d_n} \cdots \int_{-\infty}^{d_3}\int_{d_2}^{\infty} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\\ c_3 = \int_{-\infty}^{d_n} \cdots \int_{d_3}^{\infty}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\\ \vdots \\\ c_n = \int_{d_n}^{\infty} \cdots \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\\ \end{gather*} The definitions are as follows:   * $P_1, P_2, P_3, \cdots, P_n$ are known values.   * $d_1, d_2, d_3, \cdots, d_n$ are known values.   * $q(x_1, x_2, x_3, \cdots, x_n)$ is a multivariate normal distribution with mean vector $0$ and a known variance-covariance matrix.   * Thus $f, c_0, c_1, c_2, c_3, \cdots, c_n$ are all simply real numbers. * * * At most $n$ will be 35, so there will be 36 equations to solve simultaneously. What is the best method to solve this as quickly as possible using whatever programming language possible? Also I will need to change parameter inputs for $P_i$, $d_i$ and the variance covariance matrix and solve for a new set of $\lambda$'s, I will need to do this at least 2000 times, which means effectively I will need to solve a set of 35 equations for 2000 times.


title: importance of predictor variables in multiple linear regression

text: I am running multiple linear regression with R.               mod=lm(varP ~ var1 +var2+var3+var4)      The table is:               all:     lm(formula = varP ~ var1 + var2 + var3 + var4)          Residuals:         Min      1Q  Median      3Q     Max          -4.9262 -0.6985  0.0472  0.7319  4.3305           Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept)  0.700823   0.084737   8.271 1.45e-15 ***     var1      1.080172   0.175348   6.160 1.59e-09 ***     var2     -0.057803   0.007777  -7.432 5.25e-13 ***     var3     -9.924772   4.268235  -2.325   0.0205 *       var4     -0.015104   0.001290 -11.710  < 2e-16 ***     ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1          Residual standard error: 1.139 on 460 degrees of freedom     Multiple R-squared:  0.657, Adjusted R-squared:  0.654      F-statistic: 220.3 on 4 and 460 DF,  p-value: < 2.2e-16      it means that my model explains 65.4% of the variance. But now, I would like to determine the importance of each predictor. I was using:               lm.sumSquares(mod)       Is dR-sqr relevant to interpret this importance ?                             SS       dR-sqr pEta-sqr  df        F p-value     (Intercept)   88.73054 0.0510   0.1294   1  68.4015  0.0000     var4         177.88026 0.1022   0.2296   1 137.1262  0.0000     var2          71.65234 0.0412   0.1072   1  55.2361  0.0000     var1          49.22579 0.0283   0.0762   1  37.9477  0.0000     var3           7.01377 0.0040   0.0116   1   5.4069  0.0205          Error (SSE)  596.71237     NA       NA 460       NA      NA         Total (SST) 1739.76088     NA       NA  NA       NA      NA


title: How to calculate the standard deviation from a normalised value?

text: I have 2 data sets (I have more but for simplicity lets say I have 2).               Data 1     0.15239652     0.145848036     0.175981261          Data 2     0.417902092     0.342696648     0.354871141      I average these values               Data 1 : 0.158075272     Data 2 : 0.371823294      Then normalise               Data 1: 1     Data 2: 2.352191387      Now how do I calculate standard error? I'm using Excel for these calculations, by the way.


title: Trend: Linear or exponential

text: I am trying to forecast data using ARIMA. I have read articles which say that I should first take the log of the series, if my series follows an exponential trend. My question is: Is there a mathematical approach to find what kind of trend (linear or exponential) my data is following? Im working with MATLAB and R. So any builtin functions would be really helpful.


title: Incorporating random effects in the logistic regression formula in R

text: I'm trying to find the best model based on AIC using the stepwise (`direction = both`) model selection in R using the stepAIC in MASS package. This is the script i used:               stepAIC (glmer(decision ~ as.factor(Age) + as.factor(Educ) + as.factor(Child), family=binomial, data=RShifting), direction="both")      however I got this error result:               Error in lmerFactorList(formula, fr, 0L, 0L) :        No random effects terms specified in formula      I tried to add `(1|town)` to the formula since town is the random effect (where the respondents are nested) and ran this script):               stepAIC (glmer(decision ~ as.factor(Age) + as.factor(Educ) + as.factor(Child) + (1|town), family=binomial, data=RShifting), direction="both")      The result is this:               Error in x$terms : $ operator not defined for this S4 class      I hope you could help me figure out how to solve this problem. Thanks a lot.


title: can I compute the product of these two multivariate normal densities in a simple way

text: I have two multivariate normal densities and I parameterise them by the _precision_ matrix (i.e. the inverse covariance matrix). One of them has a sparse precision matrix and the other one has a simple precision matrix form with non-zero elements only along the diagonal matrix. So, we have as follows: $$ A = N(\mu_1, \Sigma_1) \\\ B = N(\mu_2, \Sigma_2) $$ the $\Sigma$'s are the precision matrices. Now I can compute a new MVN density as the product of these two with the new precision $\Sigma$ as: $$ \Sigma = \Sigma_1 + \Sigma_2 $$ So, the precisions can just be added. However, to compute the mean of this new Gaussian, I need to invert $\Sigma$, which is problematic because it is quite sparse. My question is that considering one of the matrices is a diagonal matrix, is there some simplification/trick I can use to compute the mean of the new Gaussian density. Also, to compute the normalisation of this new density, is it simply: $$ Z = \frac{det(\Sigma)^{\frac{1}{2}}}{2\pi^{\frac{n}{2}}} $$ where $n$ is the number of components in the MVN.


title: Symbolic computation in R?

text: I was wondering if it is possible to do symbolic computation in R? For example, I was hoping to get the inverse of a symbolic covariance matrix of 3D Gaussian distribution. Also can I do symbolic integration and differentiation in R? Thanks and regards!


title: Random Forest, is it a boosting algorithm?

text: Short definition of boosting: > Can a set of weak learners create a single strong learner? A weak learner is > defined to be a classifier which is only slightly correlated with the true > classification (it can label examples better than random guessing). Short definition of Random Forest: > Random Forests grows many classification trees. To classify a new object > from an input vector, put the input vector down each of the trees in the > forest. Each tree gives a classification, and we say the tree "votes" for > that class. The forest chooses the classification having the most votes > (over all the trees in the forest). Another short definition of Random Forest: > A random forest is a meta estimator that fits a number of decision tree > classifiers on various sub-samples of the dataset and use averaging to > improve the predictive accuracy and control over-fitting. As I understand Random Forest is an boosting algorithm which uses trees as its weak classifiers. I know that it also uses other techniques and improves upon them. Somebody corrected me that Random Forest is not a boosting algorithm? Can someone elaborate upon this, why Random Forest is not a boosting algorithm?


title: Is Hoeffding's bound tight in any way?

text: The inequality: $$\Pr(\overline X - \mathrm{E}[\overline X] \geq t) \leq \exp \left( - \frac{2n^2t^2}{\sum_{i=1}^n (b_i - a_i)^2} \right)$$ Is this bound (or any other form of hoeffding) tight in any sense? e.g. does there exist a distribution for which the bound is no more than a constant multiple of the true probability for every n?


title: LARS - LASSO with weights

text: I am interested in solving the following problem $$ \min_{\boldsymbol{\beta}} \left( \mathbf{y}-\mathbf{X}\boldsymbol{\beta} \right)^T W \left( \mathbf{y}-\mathbf{X}\boldsymbol{\beta} \right) + \lambda \left|\boldsymbol{\beta}\right|_1 $$ $\mathbf{y}$ is the vector of observations for each dataset. $\mathbf{X}$ is the matrix of predictors. $\boldsymbol{\beta}$ is the set of regression coefficients. $\mathbf{W}$ is a diagonal matrix filled with positive real numbers. using the LARS-Lasso approach. Is there any existing package that does this, meaning that it accepts weights for each as an input? With all the normalizing of the datasets (centering and scaling) that need to be done, I am hesitant to pre-process my data and multiply both the observation and the predictors with $\sqrt{W}$ and feed it in to the algorithm.


title: Using a regression model to make prediction: When to stop?

text: I computed a simple linear regression model from my experiment measures in order to make predictions. I have read that you should not calculate predictions for points that depart too far from the available data. However, I could not find any guidance to help me know how far I can extrapolate. For example, if I calculate the reading speed for a disk size of 50GB, I guess the result will be close to the reality. What about a disk size of 100GB, 500GB? How do I know if my predictions are close to the reality? The details of my experiment are: > I am measuring the reading speed of a software by using different disk size. > So far I have measured it with 5GB to 30GB by increasing the disk size of > 5GB between experiments (6 measures in total). My results are linear and the standard errors are small, in my opinion.


Can someone help me in understanding Sequential signed rank test? how do we relate sequential test with signed ranking test and how do we decide H0 and H1 in SSRT?


title: How to check linear independence

text: How can I check the linear independence of my variables? I have this system $Ax=b$ where A is a $N\times 4$ matrix. I want to check the linear independence between the 4 variables in matrix $A$.


title: Resources about functional distributions

text: Is there a article/textbook that treats probability distributions on functions just like basic textbooks cover the classical distributions for scalar variables? Suppose the random function $f$ has a gaussian process distribution with some mean function $m$ and covariance function $w$. I know how to generate random (functional) deviates from that distribution. What I'd like to know, for example, is given another function $g$ (that can be evaluated at an arbitrary number of points over the same range than where $f$ is defined), what is the probability that $g$ or a more extreme function was drawn from $f$'s distribution. There must be a way to define a cumulative distribution function and a probability density function for functions!


title: Joint model with interaction terms vs. separate regressions for a group comparison

text: After gathering valuable feedback from previous questions and discussions, I have came up with the following question: Suppose that the goal is to detect effect differences across two groups, male vs. female for example. There are two ways to do it:   1. running two separate regressions for the two groups, and employ Wald test to reject (or not) the null hypothesis $H_0$: $b_1-b_2=0$, where $b_1$ is the coefficient of one IV in male regression, and $b_2$ is the coefficient of the same IV in female regression.   2. pool the two groups together, and run a joint model by including a gender dummy and an interaction term (IV*genderdummy). Then, the detection of the group effect will be based on the sign of interaction and the t-test for significance.  What if Ho is rejected in case (1), i.e. group difference is significant, but the coefficient of interaction term in case (2) is statistically insignificant, i.e. group difference is insignifant. Or vice versa, Ho is not rejected in case (1), and interaction term is significant in case (2). I have ended up with this outcome several times, and I was wondering what outcome would be more reliable, and what is the reason behind this contradiction. Many thanks!


title: Does this relationship have a name?

text: Suppose we have two finite discrete probability distributions, say $A$ and $B$. For some event $e$ in both $A$ and $B$, $f_A(e) = \frac{P_A(e)}{P_A(e) + P_B(e)}$ and $f_B(e) = \frac{P_B(e)}{P_A(e) + P_B(e)}$. Does this have a name?


title: First order conditions for maximum likelihood estimator

text: Given a distribution $f(x_i,y_i;\alpha,\beta)$ and joint probability function $F(x_1,...,x_N,y_1,...,y_N;\alpha, \beta)$ The first order derivatives are $\frac{\partial F}{\partial \alpha}$ and $\frac{\partial F}{\partial \beta}$ Is the first order condition > $\frac{\partial F}{\partial \alpha}=0$ and $\frac{\partial F}{\partial > \beta}=0$ OR > $\frac{\partial F}{\partial \alpha}=0$ or $\frac{\partial F}{\partial > \beta}=0$ Do they need to satisfy at the same time or just one of them?


title: What can I do beyond Pearson correlation?

text: While checking to see whether two variables were correlated, I observed that applying Pearson correlation yielded numbers as low as 0.1 indicating no correlation. Is there anything I can do to strengthen this claim? The dataset (subset because of the posting limits) I am looking at is this:               6162.178176 0.049820046     4675.14432  0.145022261     5969.056896 0.47210138     5357.506176 0.052263122     33.796224   16.45154204     6162.178176 0.064262991     6725.448576 0.419005508     3247.656192 0.867394771     5357.506176 0.052263122     3612.97728  0.091337414     6162.178176 0.053065652     867.436416  0.129116092     556.833024  1.01107509     1517.611392 168.1484478     1517.611392 35.11570899     4675.14432  0.053902079     4182.685056 0.070289777     2808.30528  0.071929502     5969.056896 0.47193385     3247.656192 0.896646636     4387.071744 0.056985619     6273.222912 0.046547047     4387.071744 0.034875199     7946.940672 0.074997414     6612.794496 0.062001019     4675.14432  0.083205988     6612.794496 0.03326884     3321.686016 0.090917684     6194.365056 0.041166447     556.833024  1.738402642     2929.00608  0.126664128     5969.056896 0.233872792     3247.656192 0.419995196     3247.656192 0.429232627     807.890688  0.341630377     2336.767488 0.204556081     2969.23968  0.125284598     2929.00608  0.157732687     1517.611392 1666.591338     6725.448576 0.219167537     3247.656192 0.419687282     5357.506176 0.039383996     3719.193984 5.892405746     3070.628352 5.482916872     3070.628352 5.482916872     4675.14432  3.599033281     6612.794496 0.03931772     2336.767488 0.089439793     2929.00608  0.133833795     2969.23968  0.107434911     6725.448576 0.202365684     3247.656192 0.425229741     5357.506176 0.05076989     4182.685056 0.06837713     2484.827136 6.860034548     4675.14432  3.588766218     3070.628352 5.47966021     3247.656192 0.910502783     5357.506176 0.527484226     6009.290496 0.484250179     3247.656192 0.867394771     3247.656192 0.449247061     2969.23968  0.155932175     3564.69696  0.115577847     3247.656192 0.447707489     6009.290496 0.227647507     5969.056896 0.228009219     5357.506176 0.05076989     4182.685056 5.23945736     8445.837312 2.601991867     3247.656192 6.727313086     6273.222912 2.677571041     4675.14432  3.588766218     4675.14432  3.599033281     7946.940672 2.112259383     6194.365056 0.243447066     4182.685056 0.062400108     2969.23968  0.047150118     5357.506176 0.044796962     5969.056896 0.228176749     4387.071744 0.057669447     6823.61856  2.473321135     4675.14432  3.588766218     6162.178176 0.240499375     4274.417664 0.046789999     3627.461376 0.526263357     4675.14432  0.084703268     7946.940672 0.049578827     2929.00608  0.117104571     8445.837312 0.055885519     5969.056896 0.242584386     6471.172224 0.215725985     3247.656192 0.429848456     3247.656192 0.419071453     807.890688  0.313161179     4182.685056 5.231089529     7946.940672 2.753260771     4387.071744 4.981454892     6405.18912  2.622405004     4675.14432  0.097750993     3070.628352 0.092814879     3070.628352 6.443957956     6162.178176 0.046899001     6405.18912  0.196715503     5357.506176 0.104152936     867.436416  4.666624464     4675.14432  0.191651838     589.019904  1.169739758     7946.940672 10.02209571     4675.14432  0.088339519     2803.477248 0.078830674     3310.420608 0.033228406     6009.290496 0.468940552     2518.62336  1.156187164     5357.506176 0.052263122     4387.071744 0.059948871     4182.685056 0.058813895     6823.61856  0.148454957     589.019904  1.745272092     1517.611392 1429.682204     3564.69696  0.095940834     3564.69696  0.115577847     3564.69696  0.109406214     2518.62336  0.506625969     8078.90688  0.180965076     6009.290496 0.232306959     5357.506176 0.254409413     8078.90688  0.174281004     5357.506176 0.046850156     8445.837312 2.606491125     7946.940672 2.744955688     4274.417664 0.055680099     6273.222912 2.686976094     3070.628352 5.46663356     4675.14432  0.030159497     4182.685056 0.375356973     6162.178176 0.254617759     2484.827136 6.755399503     6405.18912  2.600079356     4675.14432  0.021817508     6162.178176 0.239038852     4274.417664 0.090772599     4274.417664 0.02362895     3564.69696  0.051056233     6162.178176 0.051442849     4182.685056 5.23945736     3923.580672 0.060404008     7946.940672 2.109742691     4816.766592 3.487194092     4816.766592 3.493214728     4182.685056 0.339733922     6162.178176 0.25494232     6612.794496 0.063362017     6965.240832 0.275367363     3627.461376 0.537566027     2336.767488 0.06975448     2177.442432 0.129050484     3070.628352 0.085650222     4675.14432  3.588766218     4182.685056 0.344276459     4182.685056 0.040643749     6725.448576 0.164896064     4675.14432  3.590477395     6194.365056 0.258299275     4274.417664 0.067845499     3041.66016  0.635836977     3070.628352 0.085650222     6373.00224  0.032010031     6162.178176 0.044140236     3070.628352 0.074903236     3627.461376 1.892783765     556.833024  6.99850733     589.019904  0.499134236     6162.178176 0.044951638     8078.90688  0.356112534     3247.656192 0.867702686     5357.506176 0.543536471     3247.656192 0.891104178     3247.656192 0.076978592     6725.448576 0.065125767     3564.69696  0.135775917     7946.940672 0.027683609     6725.448576 0.173073957     6725.448576 0.089064691     6725.448576 0.086239601     2484.827136 7.955885427     6823.61856  2.905496523     3070.628352 6.477501579     6162.178176 0.047710402     7946.940672 0.078646617     3321.686016 0.204414263     7946.940672 0.078646617     3299.1552   0.29431777     7946.940672 0.063420632     3041.66016  0.626960245     6162.178176 0.051442849     3070.628352 0.094117544     3564.69696  0.037029796     4885.968384 4.066747559     6273.222912 3.154518862     6162.178176 0.024179762     3321.686016 0.191770082     3321.686016 0.298643519     3627.461376 0.526263357     3627.461376 0.525712007     7946.940672 0.009185925     3070.628352 0.085650222     6373.00224  0.03436371     3564.69696  0.34028138     3564.69696  0.166914609     6725.448576 0.090848959     2808.30528  0.084392534     2336.767488 0.103989807     4675.14432  4.245858233     3564.69696  0.345611426     3564.69696  0.221337188     2808.30528  0.060178643     556.833024  1.828196167     2336.767488 0.087728026     2336.767488 0.093291267     4675.14432  4.250136175     4675.14432  4.239227421     585.801216  1.7736392     3321.686016 0.190565874     3321.686016 0.213746873     3070.628352 0.028007297     3564.69696  0.183465806     589.019904  0.329360687     6273.222912 3.165996216     4675.14432  4.221901753     6373.00224  0.036089741     4387.071744 0.223839512     3321.686016 0.165879616     7946.940672 0.124827911     3321.686016 0.174008018     7946.940672 0.134013836     2336.767488 0.076601545     556.833024  0.441784142     4816.766592 4.125173936     6194.365056 0.037453395     2808.30528  0.35216969     3299.1552   0.167012452     3299.1552   0.279465483     3299.1552   0.168527992     3321.686016 0.29412774     3321.686016 0.141795461     6373.00224  0.045347544     3627.461376 0.321712592     4387.071744 5.00151383     4387.071744 5.005844737     6405.18912  3.097332433     3299.1552   0.193079731     3321.686016 0.286601441     3321.686016 0.077370347     3299.1552   0.180652308     7946.940672 0.126589595     3299.1552   0.081232917     3041.66016  0.395178928     4182.685056 5.254041293     7946.940672 2.73979647     4274.417664 0.062464649     2929.00608  4.474213997     2808.30528  4.58034249     2808.30528  4.553992079     4675.14432  0.041282148     7946.940672 0.078646617     3299.1552   0.215206608     7946.940672 0.078646617     3321.686016 0.329049764     3321.686016 0.11981867     1517.611392 0.125855671     6965.240832 3.154951929     3247.656192 6.768265697     2929.00608  4.343794329     6373.00224  0.044719896     4885.968384 4.056514173     4675.14432  4.239227421     3070.628352 0.090860882     3321.686016 0.191770082     7946.940672 0.12507958     3321.686016 0.06532827     3321.686016 0.167384875     3321.686016 0.324232933     7946.940672 0.053731369     7946.940672 2.753260771     3719.193984 5.894556749     4182.685056 0.059292057     4675.14432  3.590477395     3299.1552   0.07789873     3321.686016 0.185447991     3321.686016 0.324232933     3321.686016 0.089713476     6965.240832 0.169986944     3627.461376 0.337150385     4387.071744 5.00151383     2803.477248 4.763370207     8445.837312 1.506304174     2848.53888  4.464745098     8445.837312 0.046768602     69.201792   7.962221556     4182.685056 0.087981762     5969.056896 0.170211144     1517.611392 0.801259141     2336.767488 0.124103062     8078.90688  1.652946395     5969.056896 2.148245564     8078.90688  1.650470812     5969.056896 2.150758524     3310.420608 0.057394519     69.201792   3.655974689     7946.940672 0.025418587     5357.506176 0.04237046     5357.506176 0.04330373     5357.506176 0.050209928     1733.263488 0.218662657     8078.90688  0.298679021     2177.442432 1.065010935     1757.403648 0.145669437     3070.628352 5.394335654     2929.00608  4.484114966     8445.837312 1.517197115     6471.172224 1.975221731     2929.00608  0.121201524     3310.420608 0.05648829     75.639168   3.344827907     1733.263488 0.133274601     1757.403648 0.169568329     2177.442432 0.118028379     3321.686016 0.164073304     3299.1552   0.296136417     3299.1552   0.14276382     75.639168   2.578029415     3041.66016  0.384000821     2969.23968  4.311878252     2929.00608  4.371107348     2848.53888  4.524775874     2929.00608  0.089791551     5357.506176 0.026878177     8445.837312 0.018470638     2177.442432 0.130887502     6823.61856  0.048507987     6273.222912 0.052923354     2808.30528  4.589600743     2969.23968  0.088574864     8445.837312 0.042624548     5357.506176 0.04237046     7946.940672 0.036366196     6273.222912 0.394534043     1733.263488 0.218662657     6273.222912 0.384650766     5357.506176 0.456929011     1733.263488 0.147698259     1733.263488 0.147698259     6194.365056 0.050691233     4834.469376 0.045920241     2969.23968  4.283251395     2808.30528  0.114303812     72.42048    3.452062179     75.639168   7.284585679     5357.506176 0.052449776     1733.263488 0.188661448     8445.837312 0.020483463     1757.403648 0.106406972     5357.506176 0.478954184     6273.222912 0.012114985     5357.506176 0.041997152     4816.766592 0.052524862     3070.628352 0.08499889     4885.968384 0.059558306     6373.00224  0.036089741     3321.686016 0.086702957     3627.461376 0.321988266     2803.477248 4.761230008     2808.30528  4.686100223     6471.172224 0.068611989     7946.940672 0.036114527     5357.506176 0.018665401     1733.263488 0.107888963     6162.178176 0.050469167     1757.403648 0.194036242     1733.263488 0.171353059     2177.442432 0.084043554     2177.442432 0.109302545     5357.506176 0.474847796     2177.442432 0.152013204     8078.90688  0.019061985     5357.506176 0.056929472     2177.442432 0.116191361     6405.18912  0.051832974     2484.827136 0.117110762     6194.365056 0.041327884     6162.178176 0.032618336     7946.940672 0.059645594     6823.61856  0.066973263     6273.222912 0.074602801     6162.178176 0.047548122     6373.00224  0.047387399     2969.23968  4.33545331     2969.23968  4.33545331     2808.30528  4.638028527     2969.23968  0.093626662     2177.442432 0.104710001     2177.442432 0.136857809     2177.442432 0.124917195     1485.424512 3.495297107     6273.222912 0.405533174     1757.403648 0.195743306     8078.90688  0.011882796     5357.506176 0.04237046     8078.90688  0.303011291     5357.506176 0.028744717     5357.506176 0.036397532     4816.766592 0.052940078     3070.628352 0.065458915     6273.222912 0.073168132     4816.766592 0.097160614     4675.14432  0.08598665     6162.178176 0.05209197     3070.628352 0.073600571     3070.628352 0.142316148     3070.628352 0.142316148     6273.222912 0.052923354     4675.14432  0.073152822     3070.628352 0.085650222     3070.628352 0.034846288     6162.178176 0.02872361     1517.611392 0.3380312     6965.240832 0.981157747     6162.178176 0.02921045     6405.18912  0.064166724     4675.14432  0.047912959     3070.628352 0.021819638     6162.178176 0.031157814     3070.628352 0.088255552     1485.424512 0.34535582     3041.66016  2.255676058     6273.222912 0.055155062     4885.968384 0.066721676     3070.628352 0.139059486     3070.628352 0.111377855     4675.14432  0.078928045     3070.628352 0.071320907     3070.628352 0.112029188     7946.940672 0.013464301     6273.222912 0.040330147     6273.222912 1.089392183     3041.66016  2.268826771     2336.767488 0.034663269     3070.628352 0.103236199     6405.18912  0.069162673     4675.14432  0.08598665     3070.628352 0.103236199     2484.827136 0.148501276     4675.14432  0.073366719     4816.766592 0.064150918     4675.14432  0.054115976     2484.827136 0.059159045     3627.461376 1.902432386     3627.461376 1.864389252     5357.506176 0.049649966     6612.794496 0.054591141     7946.940672 0.045426286     7946.940672 0.036617865     4816.766592 0.06165962     6162.178176 0.054850735     4675.14432  0.08598665     3070.628352 0.087278553     4885.968384 0.033360838     6273.222912 0.069661162     3627.461376 1.892783765     6405.18912  1.066947419     6405.18912  1.064605568     1517.611392 0.03031079     5357.506176 0.425571138     2177.442432 1.042048215     6273.222912 0.023273523     2177.442432 0.130887502     3070.628352 0.08499889     6273.222912 0.056111508     4675.14432  0.072083336     4675.14432  0.067805393     3070.628352 0.073600571     6373.00224  0.057586674     4675.14432  0.08598665     4816.766592 0.068925906     7946.940672 0.039889564     4675.14432  0.061816274     2484.827136 0.043061346     1485.424512 0.34535582     6273.222912 1.089073367     6273.222912 1.094015006     2177.442432 0.14007259     3070.628352 0.108121193     4816.766592 0.068925906     6162.178176 0.048197243     6373.00224  0.038600332     2484.827136 0.150111046     6162.178176 0.059556863     6405.18912  0.073065758     4675.14432  0.073152822     4675.14432  0.067805393     4675.14432  0.073366719     4675.14432  0.022886994     4675.14432  0.054115976     6823.61856  1.001228298     867.436416  47.76834271     5357.506176 0.041623844     5357.506176 0.052823084     1517.611392 1.775157998     75.639168   2.181409505     2929.00608  38.44000215     7946.940672 0.036995369     7946.940672 0.632696305     4675.14432  0.21838898     2177.442432 0.129968993     3310.420608 0.086695932     1517.611392 0.187136181     69.201792   1.835212591     5969.056896 18.86244376     7946.940672 0.035233684     6823.61856  0.150067005     6612.794496 0.04612271     5357.506176 0.041623844     1517.611392 0.198996925     8078.90688  2.260701883     7946.940672 0.007424241     6373.00224  0.048328871     6405.18912  0.149722355     6405.18912  0.159870377     4675.14432  0.212185963     7946.940672 0.103813535     8078.90688  2.260330546     5969.056896 18.85071661     2843.710848 39.53179701     2969.23968  59.65533911     8445.837312 20.85690187     8445.837312 20.86826841     6009.290496 29.34722495     2969.23968  59.37681663     33.796224   294.3819996     3247.656192 0.086523937     1757.403648 2.972566949     3070.628352 0.028007297     3020.738688 0.097327187     114.263424  944.37044     6009.290496 3.027645279     2969.23968  6.110992023     6009.290496 18.71252523     2929.00608  38.43897791     6009.290496 29.34273191     2929.00608  60.36484567     2969.23968  59.39635025     6823.61856  0.271996446     6405.18912  0.294448761     6405.18912  0.290233429     3247.656192 0.034794323     4885.968384 0.044617563     3070.628352 0.0511296     7946.940672 0.734118982     4675.14432  0.12192137     7946.940672 0.130993805     6273.222912 0.148567971     2969.23968  0.067020524     2177.442432 5.664902924     6009.290496 0.034280253     5969.056896 0.02831268     33.796224   173.2146171     6965.240832 1.690537382     6823.61856  0.151825603     4387.071744 0.254383804     1517.611392 48.99805075     1757.403648 0.16103301     1517.611392 1.614379025     5969.056896 3.059779848     3020.738688 0.097327187     867.436416  192.18815     2929.00608  6.235562338     8078.90688  13.93641512     3020.738688 1.112641757     1733.263488 0.397516018     2929.00608  6.2345381     2929.00608  6.218833114     2969.23968  6.110992023     867.436416  345.2783333     5969.056896 3.059779848     5969.056896 18.86244376     3020.738688 0.091037335     1733.263488 83.43047725     8445.837312 2.15668374     6009.290496 3.027645279     6009.290496 3.015663831     2929.00608  6.194934222     114.263424  1242.602357     5969.056896 3.048052702     2969.23968  6.092131976     2808.30528  0.10184078     2808.30528  0.089733834     2969.23968  37.87905731     2929.00608  38.42327292     5969.056896 29.67822272     6009.290496 29.32742228     8078.90688  21.81557513     3321.686016 0.55875239     3299.1552   0.571661497     4387.071744 0.418502388     3321.686016 0.034018869     3923.580672 0.051228716     2969.23968  6.091458403     2848.53888  0.100402351     6009.290496 18.72450667     2843.710848 62.00876581     2929.00608  60.23203612     3070.628352 0.066110247     4387.071744 0.423061237     6823.61856  0.016560128     4387.071744 0.054934137     33.796224   12.33865653     3070.628352 0.083370558     6273.222912 0.163233479     5969.056896 3.048052702     6471.172224 2.797020288     2808.30528  0.1039773     2969.23968  37.87905731     8445.837312 13.32514419     2843.710848 39.53144536     2929.00608  60.16272933     5969.056896 29.61255071     6009.290496 29.34456241     5969.056896 0.054112401     33.796224   237.3342063     2177.442432 2.527736173     6273.222912 0.24819778     3247.656192 39.86351767     6373.00224  0.028714881     7946.940672 0.04416794     6162.178176 0.049982326     4675.14432  0.195715883     2929.00608  0.067941136     3247.656192 0.040028868     5969.056896 1.514979696     33.796224   140.1635875     2177.442432 4.328013389     3247.656192 1.878585552     6471.172224 0.775284574     7946.940672 0.11740367     4387.071744 0.649636059     3627.461376 0.777403178     4387.071744 0.642797785     3020.738688 0.087064797     114.263424  1656.94317     4387.071744 0.423061237     6965.240832 0.270773121     3321.686016 0.567783948     2864.63232  0.039446598     3299.1552   0.034251193     4885.968384 0.041138211     5969.056896 21.63038521     33.796224   352.7908917     3627.461376 3.23035831     4834.469376 0.048816112     6405.18912  0.181415408     7946.940672 0.136782196     585.801216  4.832697377     585.801216  4.883909288     3627.461376 0.763619433     4182.685056 0.676837955     585.801216  4.817333804     6162.178176 0.042192873     3247.656192 3.312542758     3041.66016  3.842309589     8445.837312 0.138766585     6965.240832 0.166828402     3247.656192 0.343632433     4182.685056 0.658189647     3070.628352 0.256625     3627.461376 0.777403178     3247.656192 0.856309854     6162.178176 0.02888589     3070.628352 6.975119599     3070.628352 0.172603109     6373.00224  0.082378757     3070.628352 0.16120479     2843.710848 6.36105461     3310.420608 0.086393856     2969.23968  37.85952369     8078.90688  13.93604378     5969.056896 29.6026664     2864.63232  0.658374196     3247.656192 0.57610778     6965.240832 0.016223416     3923.580672 0.058110185     33.796224   3865.550187     6273.222912 0.163074071     6162.178176 0.049170925     6009.290496 1.860119761     2177.442432 0.1621168     5969.056896 1.21861797     3247.656192 0.073283619     3627.461376 3.23035831     4387.071744 0.236148406     3247.656192 0.351638207     4675.14432  0.176465141     2843.710848 6.361406263     2808.30528  0.1039773     6009.290496 18.72450667     8445.837312 13.32680181     6009.290496 18.72217695     2969.23968  59.6455723     2929.00608  60.34811645     2843.710848 62.00841416     6471.172224 27.25703998     5969.056896 29.49527925     3247.656192 2.339841273     114.263424  57.53372138     1757.403648 0.166154202     4675.14432  0.117643427     8445.837312 2.158341361     6009.290496 3.025315553     2848.53888  0.088466407     5969.056896 18.85071661     6471.172224 17.37351381     2929.00608  38.39937403     2969.23968  37.86019726     8078.90688  21.80851972     33.796224   6.272890131     6823.61856  0.228178053     3020.738688 0.001655224     114.263424  6417.215364     3070.628352 0.097048541     3070.628352 0.060573921     3247.656192 0.580726496     6373.00224  0.042209306     4387.071744 0.429899512     3321.686016 0.034018869     6373.00224  0.031539295     114.263424  91.70038524     114.263424  1235.758522     3070.628352 0.066110247     3321.686016 0.567783948     3321.686016 0.567783948     4387.071744 0.0257575     6405.18912  0.017641946     6373.00224  0.044092249     1733.263488 75.69131924     114.263424  2.608008666     3041.66016  3.841323286     3627.461376 3.241385305     7946.940672 0.122311219     6273.222912 0.185231741     3070.628352 6.947763635     4675.14432  0.112723793     3627.461376 3.124499154     4387.071744 0.25210438     6273.222912 0.174551425     7946.940672 0.131245475     4182.685056 0.676837955     3070.628352 0.896559168     3070.628352 0.921961135     4182.685056 0.665840235     3070.628352 0.881252854     4834.469376 0.056676334     4675.14432  0.031015085     2484.827136 0.226172675     3020.738688 1.086158168     1733.263488 0.171930005     3070.628352 0.053083598     2864.63232  0.669195829     3299.1552   0.565599339     4387.071744 0.0257575     3321.686016 0.019869428     114.263424  1126.89604     6725.448576 1.729252684     6273.222912 0.165146371     8445.837312 0.135688145     6405.18912  0.15955813     4387.071744 0.649636059     4182.685056 0.66727472     585.801216  4.754172446     3070.628352 0.881252854     3070.628352 0.080765228     6162.178176 0.091201517     6373.00224  0.075945368     4675.14432  0.105237393     2484.827136 0.183916214     6162.178176 0.082925223     3070.628352 6.476850247     2484.827136 7.998141888     1517.611392 1.908920831     3020.738688 0.461145483     5357.506176 0.117218717     5357.506176 0.766961318     4675.14432  0.13689417     556.833024  1.111643838     6612.794496 0.042947048     4675.14432  0.088339519     4834.469376 0.037439476     4387.071744 0.058353274     1517.611392 1311.730401     4182.685056 0.085830034     589.019904  1.769040389     556.833024  0.375336934     3564.69696  0.115577847     5357.506176 0.054129662     729.032832  0.73659234     5357.506176 0.27288816     3247.656192 0.43015637     3247.656192 0.419379368     3247.656192 0.428924713     5357.506176 0.046850156     3627.461376 0.337150385     3070.628352 5.482916872     6273.222912 2.690323656     2336.767488 0.088155968     8445.837312 0.039901313     6162.178176 0.051442849     3070.628352 0.108121193     6612.794496 0.045517822     6162.178176 0.071403323     3070.628352 0.921961135     585.801216  4.883909288     585.801216  4.883909288     3247.656192 0.866778943     7946.940672 0.066692331     6273.222912 0.083689039     6273.222912 0.083051409     4675.14432  0.105237393     7946.940672 0.056751399     6162.178176 0.035539381     589.019904  6.453092628     3070.628352 0.773131662     4387.071744 0.634819798     4834.469376 0.053366767     6273.222912 0.083689039     3070.628352 0.147852474     4675.14432  0.097750993     3070.628352 0.095745876     6405.18912  3.100767148     7946.940672 2.48988898     3070.628352 6.452750945     556.833024  6.763248295     589.019904  0.483854617     3020.738688 0.420426965     6162.178176 0.022881519     5357.506176 0.09034054     4834.469376 0.108595165     4816.766592 0.105049724     3070.628352 0.160227792     4675.14432  0.093259153     3070.628352 0.090860882     3070.628352 6.477501579


title: Homogeneity of variance in regression

text: I have carried out simple linear regression and I am now checking the models meets the assumption of homogeneity of variance: • am I correct in concluding that the Levenes tests which gave a p<0.05 indicates a violation of homogeneity of variance? • some models contain 2 or 3 explanatory variables, but in some models, only 1 explanatory variable gave p<0.05 in the Levenes test. Therefore, should only those explanatory variables which gave p<0.05 be corrected? • is transformation the best place to start when correcting for violation of homogeneity?


title: GDA and LDA terminology

text: Can the terms LDA (Linear Discriminant Analysis) and GDA (Gaussian Discriminant Analysis) be used interchangeably? Do they often refer to the same thing?


title: Literature on the algorithm for optimal splitting in the growing of classification trees

text: In ESL, Section 9.7, there is a paragraph stating that the computation time of a split in the growing of a classification (or regression) tree _typically_ scales like $p N \log N$ where $p$ is the number of predictors and $N$ is the number of samples. A naive approach results in a $pN^2$ scaling, and I haven't been able to find any references to the literature that explains the details for the splitting part of the algorithm and how one achieves a _typical_ $p N \log N$ scaling. In the naive approach the optimal split for a given variable, after an initial ordering of the observed values, is sought among the $N-1$ midpoints between the observed values, and computing the loss for each split can be done in time that scales like $N$. I could (and will probably) study the source code for some of the implementations I know, but a literature reference would be nice $-$ in particular regarding the time complexity.


title: Does my multilevel modelling create false replication and spurious precision at the highest level?

text: I am fitting a 3-level model in Bugs. The data are in long format. The data consists of 885 Treatments in 1:3 COURSEs of treatment in 108 Patients (MRNs).               model{          for(DMRN in 1:NMRN) {            RMRN[DMRN] ~ dnorm(meanMRN[DMRN], TMRN)       meanMRN[DMRN] <- intercept            for(DCOURSE in SMRN[DMRN]:(SMRN[DMRN+1]-1)){              RCOURSE[DCOURSE] ~ dnorm(meanCOURSE[DCOURSE], TCOURSE)         meanCOURSE[DCOURSE] <- RMRN[DMRN]              for(Dobservations in SCOURSE[DCOURSE]:(SCOURSE[DCOURSE+1]-1)){                totalEEG[Dobservations] ~ dpois(meanobservations[Dobservations])           log(meanobservations[Dobservations]) <- RCOURSE[DCOURSE] +inprod(betaobservations[] , Xobservations[Dobservations,])              }#observations       }#COURSE     }#MRN               # priors          intercept ~ dflat()     betaobservations[1] ~ dflat()     betaobservations[2] ~ dflat()     betaobservations[3] ~ dflat()          TMRN <- pow(SDMRN, -2)     SDMRN ~ dunif(0, 100)     TCOURSE <- pow(SDCOURSE, -2)     SDCOURSE ~ dunif(0, 100)          } # model      I am worried that my covariates in the model - one of which (Age) varies only at the top level - are replicated in long format and therefore "counted more" in BUGS and giving me spuriously precise posterior distributions. However, obviously it is not contributing any variation at the lowest level (within patient).               bgs.toteeg<-glmmBUGS(formula = totalDuration ~ AnaDose+ AgeYrs + log.workDose, effects=c("MRN", "COURSE"), family="poisson", data=book,  modelFile="model.bug")     source("getInits.R")     startingValues=bgs.toteeg$startingValues     myResult = bugs(bgs.toteeg$ragged, getInits ,model.file = "model.bug", n.chain = 3, n.iter = 10000,  n.burnin = 1000, parameters.to.save = names(getInits()),n.thin = 10, debug=T)      AgeYrs is included as a covariate here but obviously doesn't vary across patient-treatment, but only between patients. All of the other covariates are different values for each treatment. The results back in R are:               > mysummary$scalars[, c("mean", "2.5%","97.5%")]                     mean        2.5%    97.5%     intercept 3.36308519 3.220175000 3.473550     SDCOURSE  0.04677978 0.001696725 0.094267     SDMRN     0.31222407 0.265697500 0.357700     > signif(mysummary$betas[, c("mean","2.5%", "97.5%")], 3)                    mean     2.5%   97.5%     MHXKG        0.0605 -0.01360 0.14200     AgeYrs       0.0027  0.00162 0.00407     log.workDose 0.0457  0.03170 0.06000      How does the model prevent spurious replication of repeated measures which in this case, are exactly the same across patient treatments, but only differ between patients, therefore at the COURSE" level there will be a low or no variability, at the Patient level there will be, variability according to : ![Histogram of Age \(with replication\)](http://i.stack.imgur.com/4qMdT.png) Also here is the same model fit with ML in lmer:                   Linear mixed model fit by REML          Formula: totalEEG ~ AgeYrs + MHXKG + log.workDose + (1 | MRN/COURSE)             Data: book            AIC  BIC logLik deviance REMLdev          7015 7048  -3500     7006    7001         Random effects:          Groups     Name        Variance   Std.Dev.            COURSE:MRN (Intercept) 2.7923e-11 5.2842e-06          MRN        (Intercept) 1.8464e+02 1.3588e+01          Residual               2.4688e+02 1.5713e+01         Number of obs: 817, groups: COURSE:MRN, 114; MRN, 103              Fixed effects:                      Estimate Std. Error t value         (Intercept)  31.12003    7.35545   4.231         AgeYrs        0.10522    0.05946   1.769         MHXKG         4.57064    3.71167   1.231         log.workDose  1.14654    0.82651   1.387              Correlation of Fixed Effects:                     (Intr) AgeYrs MHXKG          AgeYrs      -0.504                       MHXKG       -0.562  0.002                log.workDos -0.671 -0.006  0.093      And here is a model that has -as far as I can see Age varying only at the Patient level. Fit with Stata:               . xtmixed totalEEG MHXKG workDose, || MRN: AgeYrs, covariance(exchangeable)          Performing EM optimization:           Computing standard errors:          Mixed-effects REML regression                   Number of obs      =       817     Group variable: MRN                             Number of groups   =       103                                                          Obs per group: min =         1                                                                    avg =       7.9                                                                    max =        31                                                               Wald chi2(2)       =     10.60     Log restricted-likelihood = -3504.9961          Prob > chi2        =    0.0050          ------------------------------------------------------------------------------         totalEEG |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------            MHXKG |   2.434717   3.560351     0.68   0.494    -4.543442    9.412877         workDose |  -.0102969   .0033165    -3.10   0.002    -.0167971   -.0037966            _cons |   48.81418   4.001097    12.20   0.000     40.97217    56.65618     ------------------------------------------------------------------------------          ------------------------------------------------------------------------------       Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]     -----------------------------+------------------------------------------------     MRN: Exchangeable            |                 sd(AgeYrs _cons) |   .1967401   .0184842      .1636517    .2365186               corr(AgeYrs,_cons) |   .9999981   .0007988            -1           1     -----------------------------+------------------------------------------------                     sd(Residual) |   15.86585   .4234722       15.0572    16.71794     ------------------------------------------------------------------------------     LR test vs. linear regression:       chi2(2) =   160.00   Prob > chi2 = 0.0000          Note: LR test is conservative and provided only for reference.      I'd be delighted if somebody could clarify how this works .. or whether I'm getting it wrong. edit( +2hrs) Having looked back over PinHeiro & Bates I wonder whether:               check.fit<-lmer(totalEEG ~ MHXKG + workDose +(AgeYrs|MRN), data=book, family=gaussian(link = "identity"))      is more correct? It produces the same parameter values as above STATA output. However, how would one specify this model in BUGS - or in glmmBUGS to get into BUGS??


title: Zero inflated models - "true zero" vs. "excess zero"

text: I am trying to decide if zero inflated poisson is appropriate for my data vs. a Poisson hurdle model. In background reading between the two I've run across a statement saying that a zero inflated model attempts to distinguish between true zeros and excess zeros. I'm having a problem understanding what is the different between those two zeros. Can anyone explain what those two types of zeros mean in the context of zero inflated modeling?


title: At what value of mean and variance should I throw data away?

text: I have some score values that are output from a program. There are about 10 such values. The data set is a measure of the "quality" of a speech waveform received over a mobile phone and landline channel. The waveform is passed through an algorithm and a score is received of its quality relative to a "golden" waveform (that gets a 100). My task is to make adjustments to the algorithm that bring the scores from the mobile and landline channels closer. I am afraid that is the most detail I can give out about the task. Please find below some of these scores:               Mobile: 52 66 69 54 88       Landline: 60 57 72 49 75        When I compute the mean and variance of this small data-set, I get very high variance (expected from a small data-set). My question is:   1. Are data-sets with very high variances usually rejected?    2. If so, is there a relation between the number of elements in the data-set, its mean and variance such that I can take a look and say "Ahhh... the variance is too high in that one (according to some relation that I do not know of), I must reject(?) this data". P.S: Please let me know in case my question does not make sense. I will try to elaborate.


title: Testing "closeness of data"

text: I have data from two different sets (each with fields named _Mobile_ and _Landline_ ). There are about 5 elements from each field in each set. They go something like so: Set A:               Mobile: 45 61 78 65 71       Landline: 56 71 75 98 12      Set B:               Mobile: 67 71 75 64 23       Landline: 9 34 56 21 22      I need to prove/disprove that the figures in the _Mobile_ field are closer to the figures in the _Landline_ field in _Set A_ as against _Set B_. How do I go about doing that? Would a distance measure between the two suffice? If so, what distance measure? P.S: There may not be enough data in the two sets, but that is all I have at the moment. _**Edit 1_** : Here is what I have tried in the interim. I made a single vector out of the `Mobile` and `Landline` scores for `Set A` and did the same for `Set B`. I then compute mean and variance for the two vectors. In case the variance/standard deviation for the vector from `Set A` is lesser than that of `Set B`, I can say that the `Mobile` scores in `Set A` are "closer" to `Landline` scores in `Set A` than `Set B`. Is this approach correct? Why or why not?


title: Plotting n number of lists along the same domain in R

text: I generate 3 or 4 lists (really tables? but read.csv gives me lists) all with a similar Time column on the far left. I have to format the time column since they all format it their own way (some including dates etc.), now I'm wondering if I should format it in # seconds instead of HH:MM:SS.mm. I thought I would mention the formatting as it could affect the question. Now how do I get the 3 lists to plot the rest of their columns along the time column, the 3 lists don't necessarily have the same times within the rows. They certainly don't have the same number of columns either. Some columns are numeric, others are not. So it's sort of an infographic where it's not (I was thinking hover over show text info or just print it), and a xy plot (for now) if the column is numeric. Does this make any sense? Here are list examples. list 1 Time Event 12:00:01.00 Some event List 2 Time Info More Columns 12:01:11.45 Some info More info at this time List 3 Time Page Faults/sec 12:05:37.88 120 I'm new to stats so I just want a decent visual of these lists.


title: Averaging Variables to create general equation

text: My data best fit is the 3 parameter inverse gamma distribution(thereby giving a, b, y variables with each event having a specific value for each variable) but I am not sure how to create a general equation for these equation...do I average each variable? or is there some other method? I've seen one researcher use MLE but he doesn't fully explain what he did to create the generalized equation. Is it possible to use MLE to find the mean of an equation with multiple variables and multiple values for each variable?


title: Difference between SUR and SEM

text: SUR means Seemingly unrelated regression, and SEM means simultaneous equation model. What is the difference between them?


title: What is the Dirichlet equivalent of a Beta (1,1) distribution?

text: If parameter p ~ Beta(1,1), this would reflect we know nothing about parameter $p$. Generalizing to the multivariate case, how would the same be said about a vector $P$ of probability parameters $p_i$?


title: How to adjust for a mid-study change in diagnostic protocol?

text: The setup for this question is as follows. Lets say we have count and denominator data for six different study sites, by month for 12 months. They're all using a uniform diagnostic standard at the time, so theoretically their counts are all roughly the same. The goal of this is to get an estimate of the typical monthly case count per person-time each study site gets, though we allow it to vary by month in the event that it changes dramatically. So essentially, a pretty standard Poisson regression estimation of an incidence-density. _Except_ two of the sites change their diagnostic procedure to a more accurate one somewhere in the middle of the study period. Lets say Site 1 changes 2 months in, and Site 2 changes 8 months in, for the sake of the example. The new procedure is believed to be of improved accuracy, but there's no direct data from the sites to produce our own accuracy measures, we'd have to go off published sensitivity/specificity figures from the two tests in the literature. I'm trying to figure out a way to get around just out and out discarding the data from after the protocol change. Any ideas? An update: One of the sites in question is about 3% of the total study person- time but about 10% of the cases, both because of the better test, and by the looks of it, an outbreak at the site at the time. So while "toss it" is still an option, it's not my favorite option.


title: Logistic regression: if only some classes of a categorical variable appear significant

text: I am performing (rare events) logistic regression analyses in R and want to test several categorical variables consisting of more than two classes. I understood that I can do this by using factor(). However, I am not sure what to do in case some, but not all categories are significant. Can I (1) just leave out the ones that are not significant, assuming their coefficients equal 0, or (2) do I have to include them all in the model, or (3) do I somehow have to recalibrate the logistic regression model with only the significant dummies? Thanks in advance!


title: R: power calculation for beta coefficient

text: I am interested in calculating the power for the beta effect size, for example               Call:     lm(formula = log1p(y) ~ x)          Residuals:         Min      1Q  Median      3Q     Max      -0.5684 -0.1881 -0.0413  0.1494  1.2312           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  0.59725    0.02460  24.279   <2e-16 ***     x           -0.06087    0.05514  -1.104     0.27         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.2551 on 1667 degrees of freedom     Multiple R-squared:  0.0007306, Adjusted R-squared:  0.0001312      F-statistic: 1.219 on 1 and 1667 DF,  p-value: 0.2697      So, the beta (effect size) for x here is -0.06087. Can I calculate the power of detecting current effect size with two-tailed significant level of 0.05 in the following equation               Z(power)=abs(beta)/se-1.96=0.06087/0.05514-1.96=-0.85     Power=pnorm(Z(power))=0.196      So can I say with current sample size, I have 20% power to detect the significant associations between x and y with current effect size.


title: computing a p-value for bootstrapped-mean difference

text: I have read a number of times about computing p-values in bootstrapping. Unfortunately, I could not use none of them. Maybe you can help me, please? I obtained 10000 bootstrap resamples of the difference between the means of two groups of data. I also got the confidence interval for the bootstrap-mean difference. Now I would like to compute a p-value too. How can I do it?


title: Covariance with a common term

text: So I've been digging into the math of correlations and covariances and one of my articles mentioned finding the covariance by using the common component, specifically with regards to, say, $cov(A + B, C + B)$. I know that $$Cov(A + B, C + D) = cov(A, C) + cov(A, D) + cov(B, C) + cov(B, D)$$ but it seems like there would be a simplification in this case, but I don't understand how just using the "common component" would come into play. Any pointers?


title: Estimation accuracy of precision matrix

text: I have a couple of questions related to estimation of high-dimensional precision matrix (inverse of the covariance matrix) in the case where p is close to 100 and n < p. As a measure of estimation accuracy I am calculating a 'relative squarred error' defined in the following manner: $$ RSE(\Sigma^{-1},\hat{\Sigma^{-1}}) = \frac{||{\Sigma^{-1}-\hat{\Sigma^{-1}}}||}{||\Sigma^{-1}||} $$ where, $||A|| = \sum_{i,j}a_{ij}^2 $ is the Frobenius norm of $A$. Now, my questions are, a) Is there any known result for lower bound of Frobenius norm of a precision matrix? or for that matter, any general matrix? b) What are the common ways of showing the estimation accuracy of a method. I think we can of course calculate various norms such as the spectral norm etc, but, is there a graphical way of doing this? I appreciate any kind of help and apologies in advance if this is not the right place or way to ask such questions.


title: Estimating finite sample bias for Instrumental Variables

text: Are there ways to **_estimate_** the finite sample bias with instrumental variables? I guess this would be conditional on assuming some structure to the problem and also would involve simulation, but, at least in applied econometrics papers, I have never seen people do it. _PS: If possible, the answer could be illustrated with R code._


title: Calculating sampling bias when estimating population proportions

text: Assume a true population of size N consisting of red and blue balls, where p represents the frequency of red in the population and q represents the frequency of blue in the population. Let's say that N=1000, p=0.7 and q=0.3. If I take a sample from the population and observe the frequencies in the sample to infer the true population proportions, the smaller the sample the more variance there will be. Is there a way to calculate the extent to which different sample sizes will over or underestimate the true proportions?


title: Random effect on the intercept in models with categorical predictor

text: I'm doing a mixed model analysis of a model with a categorical predictor (using glmer, binomial model). To my surprise, the variance I obtain for the intercept does not change when I change the reference level of this categorical predictor. This goes against my understanding that as the intercept is the log-odd for the reference level of the categorical predictor, then the variance for the random intercept should change when changing the reference level. What am I missing here ? Many thanks for your assistance. Update: @charles Many thanks for your reply. I still don't get it. The following code, in which I simulate data similar to those that I'm analyzing, might help pursuing the discussion.               require(lme4)     set.seed(123)          # Create 6 groups of 100 observations each     group <- gl(6,100)     # Create a categorical predictor     category <- as.factor(round(runif(600,1,5))) # unbalanced design per group     # Create a binary response variable     used <- runif(600)      used[1:100] <- used[1:100]+runif(100,0.1,0.3) # to increase the number of 1s in the first category, just to create a bit of variability      used[101:200] <- used[101:200]-runif(100,0.1,0.3) # to decrease the number of 1s in the second category, just to create a bit of variability      used <- round(used)     used[used>1] <- 1; used[used<0] <- 0          # Fit two models with different reference categories     mod.ref1 <- glmer(used~category+(1|group),family="binomial")     category <- relevel(category,"3") # change reference category     mod.ref3 <- glmer(used~category+(1|group),family="binomial")          # Compare coefficients and random effects for both models     coef(mod.ref1); ranef(mod.ref1)     coef(mod.ref3); ranef(mod.ref3)      And the results are:               coef(mod.ref1)     $group       (Intercept)   category2  category3  category4   category5     1  0.72458362 -0.02654181 -0.2808615 0.02515119 -0.09999195     2 -0.57971784 -0.02654181 -0.2808615 0.02515119 -0.09999195     3  0.14790658 -0.02654181 -0.2808615 0.02515119 -0.09999195     4  0.14356195 -0.02654181 -0.2808615 0.02515119 -0.09999195     5 -0.02681879 -0.02654181 -0.2808615 0.02515119 -0.09999195     6  0.11810340 -0.02654181 -0.2808615 0.02515119 -0.09999195          ranef(mod.ref1)     $group       (Intercept)     1  0.63660977     2 -0.66769169     3  0.05993273     4  0.05558810     5 -0.11479264     6  0.03012955          coef(mod.ref3)     $group       (Intercept) category1 category2 category4 category5     1   0.4436543 0.2808285 0.2543196 0.3060796 0.1809765     2  -0.8605623 0.2808285 0.2543196 0.3060796 0.1809765     3  -0.1329837 0.2808285 0.2543196 0.3060796 0.1809765     4  -0.1373277 0.2808285 0.2543196 0.3060796 0.1809765     5  -0.3076979 0.2808285 0.2543196 0.3060796 0.1809765     6  -0.1627864 0.2808285 0.2543196 0.3060796 0.1809765          ranef(mod.ref3)     $group        (Intercept)     1  0.63656425     2 -0.66765237     3  0.05992621     4  0.05558221     5 -0.11478796     6  0.03012352      and the variance of the random intercepts is               Random effects:      Groups Name        Variance Std.Dev.      group  (Intercept) 0.18119  0.42567      Number of obs: 600, groups: group, 6      for both models. So, (1) variance for the random intercepts is the same in both models, (2) random intercepts appear to be the same in both models, (3) intercepts vary between models as they refer to the reference category. I don't understand how (2) and (3) can be true simultaneously, and I'm wondering if I'm misunderstanding the meaning of random intercepts when categorical predictors come into play. Could anyone put down the equation underlying this model ?


title: Asymptotic distribution for t-test with time series

text: I once read the following statement: Given the following equation $ Y_t= c + a Y_{t-1} + \theta_t$   1. If $c\neq 0$, then the t-ratio for testing $H_0 : a=1$ is asymptotically normal.    2. If $c=0$, then the t-ratio for testing will converge to another nonstandard asymptotic distribution. But I do not know how to get the above two conclusions.


title: Using Mutual Information for Binary Logistic Regression Variable Selection

text: In addition to `proc varclus`, `randomForest`, and assessing multicollinearity among potential predictor variables, I am seeking other methods of variable selection in lieu of using stepwise methods for building more parsimonious binary logistic regression models from a wide array of potential predictor variables. I have done some research into other methods such as Mutual Information (MI), and I have two questions in regards to its use: 1) Has anyone used MI for binary logistic regression variable selection? If so, what are your thoughts on its application? 2) Does anyone know how to calculate MI using either Base SAS or R for potential predictor variables in reference to the outcome of interest? Any help or references in this area would be greatly appreciated! Thanks!


title: How to deal with high correlation among predictors in multiple regression?

text: I found a reference in an article that goes like: > According to Tabachnick & Fidell (1996) the independent variables with a > bivariate correlation more than .70 should not be included in multiple > regression analysis. **Problem:** I used in a multiple regression design 3 variables correlated >.80, VIF's at about .2 - .3, Tolerance ~ 4- 5. I cannot exclude any of them (important predictors and outcome). When I regressed the outcome on the 2 predictors which correlated at .80, they remained both significant, each predicted important variances, and these same two variables have the largest part and semipartial correlation coefficients amongst all 10 variables included (5 controls). **Question:** Is my model valid despite high correlations? Any references greatly welcomed! * * * Thank you for the answers! I did not use Tabachnick and Fidell as a guideline, I found this reference in an article dealing with high collinearity amongst predictors. So, basically, I have too few cases for the number of predictors in the model (many categorical, dummy coded control variables- age, tenure, gender, etc) - 13 variables for 72 cases. The Condition Index is ~ 29 with all the controls in and ~23 without them (5 variables). I cannot drop any variable or use factorial analysis to combine them because theoretically they have sense on their own. It is too late to get more data. Since I am conducting the analysis in SPSS perhaps it would be best to find a syntax for ridge regression (although I haven't done this before and interpreting the results would be new to me). If it matters, when I conducted stepwise regression, the same 2 highly correlated variables remained the single significant predictors of the outcome. And I still do not understand if the partial correlations which are high for each of these variables matter as an explanation for why I have kept them in the model (in case ridge regression can't be performed). Would you say the "Regression diagnostic: identifying influential data and sources of collinearity / David A. Belsley, Edwin Kuh and Roy E. Welsch, 1980" would be helpful in understanding multicollinearity? Or might other references be useful?


title: Incorporating systematic uncertainty in interval estimation

text: I am trying to fit some binomially distributed data to a model. There is a reasonable systematic error associated with the measurements. This uncertainty comes from imperfect theoretical knowledge of one of the parameters assumed by the model. I can calculate this number by 3 or 4 different methods and model the uncertainty. The model is quite complex and itself is a function of 10 or so paramters of the experiment as well as the calculation method. What I don't know is how to incorporate this into the binomial maximum likelihod estimation expression. So my model without systematics is $$ L(\theta | x, k) = \Pi_i^N p(\theta, x_i)^{k_i}(1-p(\theta, x_i))^{(n_i- k_i)}$$ I assume that it is some kind of muliplicative factor so I end up with $$ L(\theta | x, k) = \Pi_i^N p(\theta, x_i)^{k_i}(1-p(\theta, x_i))^{(n_i- k_i)}\times L(syst)$$ Is there a general way of building a systematic uncertainty model into an MLE minimisation/maximisation? More specifically is there a way of weighting different contributions to the systematic error model? Whatever the answer to my above question I suspect I will end up with a binomial maximum likelihood expression multipled by a bayesian systematic error model. Would such a hybrid approach be invalid in anyway?


title: Modelling aggregated binary outcomes (proportions) with individual-level covariates

text: Suppose the set-up is that of typical logistic regression with individual- level covariates and binary outcomes, but instead of observing those outcomes, one only observes the proportion of 'successes' among known groups of the individuals. So the outcome is an aggregation of individual-level outcomes, but the individual-level covariates are all observed. How does one model in this framework?


title: Is there a minimum data set for making classification decisions with Naïve Bayes?

text: I have objects from two class -- 11 objects in each class. Where each object has 3 features. with such small dataset, does it make any sense to train & use Naïve Bayes to classify new objects (test cases) ?


title: Expectation of Truncated & Random Variable

text: I have what appears to be a relatively simple question, but am struggling to understand how to go about answering it. The general question is as follows: What is the expected value of $S_{I}$, where: $S_{I} = S$ if $S <3000 $ $S_{I} = 3000$ if $S >3000$ where $S$ is a compound distribution (details not necessary for my problem here) My initial attempt is as follows: $E[S_{I}] = E[E[S_{I}|S]] = E[S * P(S\leq3000) + 3000*P(S>3000)] = E[S]*P(S\leq 3000)+3000*P(S>3000).$ Now, from the definition of $S_{I}$, it is clear that we should have $E[S_{I}]<3000$. For this particular problem $E[S] = 4000$, and hence my proposed method for solving is wrong. So I thought that my $E[S]$ in the line above should actually be $E[S|S\leq 3000]$ Is this assumption correct? If so, is it true that $E[S|S\leq 3000] = E[S]*P(S\leq 3000)$ ?? Thanks for any help.


title: GLM after model selection or regularization

text: I would like to pose this question in two parts. Both deal with a generalized linear model, but the first deals with model selection and the other deals with regularization. **Background:** I utilize GLMs (linear, logistic, gamma regression) models for both prediction and for description. When I refer to the " _normal things one does with a regression_ " I largely mean description with (i) confidence intervals around coefficients, (ii) confidence intervals around predictions and (iii) hypothesis tests concerning linear combinations of the coefficients such as "is there a difference between treatment A and treatment B?". Do you legitimately lose the ability to do these things using the normal theory under each of the following? And if so, are these things really only good for models used for pure prediction? **I.** When a GLM has been fit via some model selection process (for concreteness say its a stepwise procedure based on AIC). **II.** When a GLM has been fit via a regularization method (say using glmnet in R). My sense is that for I. the answer is technically that you should use a bootstrap for the " _normal things one does with a regression_ " but no one really abides by that. Thank you. **Add:**   After getting a few responses and reading elsewhere, here is my take on this (for anyone else benefit as well as to receive correction).      **I.**   **A) RE: Error Generalize.** In order to generalize error rates on new data, when there is no hold out set, cross validation can work but you need to repeat the process completely for each fold - using nested loops - thus any feature selection, parameter tuning, etc. must be done independently each time. This idea should hold for any modeling effort (including penalized methods).      **B) RE: Hypothesis testing and confidence intervals of GLM.** When using model selection (feature selection, parameter tuning, variable selection) for a generalized linear model and a hold out set exists, it is permissible to train the model on a partition and then fit the model on the remaining data or the full data set and use that model/data to perform hypothesis tests etc. If a hold out set does not exist, a bootstrap can be used, as long as the full process is repeated for each bootstrap sample. This limits the hypothesis tests that can be done though since perhaps a variable will not always be selected for example.      **C) RE: Not carrying about prediction on future data sets** , then fit a purposeful model guided by theory and a few hypothesis tests and even consider leaving all variables in the model (significant or not) (along the lines of Hosmer and Lemeshow). This is small variable set classical type of regression modeling and then allows the use of CI's and hypothesis test.      **D) RE: Penalized regression.** No advice, perhaps consider this suitable for prediction only (or as a type of feature selection to then apply to another data set as in B above) as the bias introduced makes CI's and hypothesis tests unwise - even with the bootstrap.


title: Estimating R-squared and statistical significance from penalized regression model

text: I am using the R package _penalized_ to obtain shrunken estimates of coefficients for a dataset where I have lots of predictors and little knowledge of which ones are important. After I've picked tuning parameters L1 and L2 and I'm satisfied with my coefficients, is there a statistically sound way to summarize the model fit with something like R-squared? Furthermore, I'm interested in testing the overall significance of the model (i.e. does R²=0, or do all the =0). I've read through the answers on a similar question asked here, but it didn't quite answer my question. There's an excellent tutorial on the R package that I'm using here, and the author Jelle Goeman had the following note at the end of the tutorial regarding confidence intervals from penalized regression models: > It is a very natural question to ask for standard errors of regression > coefficients or other estimated quantities. In principle such standard > errors can easily be calculated, e.g. using the bootstrap. > > Still, this package deliberately does not provide them. The reason for this > is that standard errors are not very meaningful for strongly biased > estimates such as arise from penalized estimation methods. Penalized > estimation is a procedure that reduces the variance of estimators by > introducing substantial bias. The bias of each estimator is therefore a > major component of its mean squared error, whereas its variance may > contribute only a small part. > > Unfortunately, in most applications of penalized regression it is impossible > to obtain a sufficiently precise estimate of the bias. Any bootstrap-based > cal- culations can only give an assessment of the variance of the estimates. > Reliable estimates of the bias are only available if reliable unbiased > estimates are available, which is typically not the case in situations in > which penalized estimates are used. > > Reporting a standard error of a penalized estimate therefore tells only part > of the story. It can give a mistaken impression of great precision, > completely ignoring the inaccuracy caused by the bias. It is certainly a > mistake to make confidence statements that are only based on an assessment > of the variance of the estimates, such as bootstrap-based confidence > intervals do.


title: Detecting a given face in a database of facial images

text: I'm working on a little project involving the faces of twitter users via their profile pictures. A problem I've encountered is that after I filter out all but the images that are clear portrait photos, a small but significant percentage of twitter users use a picture of Justin Bieber as their profile picture. In order to filter them out, how can I tell programmatically whether a picture is that of Justin Bieber?