title: Mean of a function of $K$ Gaussian random variables

text: Consider $K$ Gaussian random variables $Z_1, Z_2,\ldots,Z_K$ that are mutually independent with $Z_1 \sim \mathcal{N}(\mu,\sigma^2)$ and the remaining random variables $\sim \mathcal{N}(0,\sigma^2)$. Let $m_1({\bf Z})$ and $m_2({\bf Z})$ be functions of ${\bf Z} = (Z_1,\ldots,Z_K)$ such that $m_1({\bf Z})$ is the random variable with the largest magnitude and $m_2({\bf Z})$ is the random variable with the second largest magnitude among the $K$ random variables. I am trying to find the mean of the ratio $R = \frac{|m_2({\bf Z})|}{|m_1({\bf Z})|}$. **Q1:** The event $m_1({\bf Z})=0$ is equivalent to the vector ${\bf Z}= {\bf 0}$ and hence it occurs with zero probability. So, I have assumed that $R$ is well-defined. The range of $R$ is $(0,1]$. So, the mean of $R$ must exist. Am I correct in assuming that the mean exists? **Q2:** If the mean exists, the expression for the mean is given by $$ E[R]= \int\limits_{-\infty}^{\infty}\ldots\int\limits_{-\infty}^{\infty} \frac{|m_2({\bf z})|}{|m_1({\bf z})|} \frac{1}{(2\pi \sigma^2)^{K/2}} \exp \\{-\frac{(z_{1}-\mu)^2+\sum_{i=2}^{K}z_{i}^2}{2\sigma^2}\\} dz_{1}\ldots dz_{K}. $$ In general, for $ K>2$ this integral cannot be evaluated in closed form and numerical integration is required. Is numerical integration even feasible for values of $K >2 $? **Q3:** In my attempt to answer Q2, I expressed the $K$-dimensional integral in terms of repeated one-dimensional integrals. My approach was to consider the following three cases.   1. $Z_1$ has the largest magnitude i.e., $m_1({\bf Z}) = Z_1$. Let $P_1$ be the conditional probability of $m_1({\bf Z}) = Z_1$ and $E_1$ be the conditional mean of symbol ratio when $m_1({\bf Z}) = Z_1$.   2. $Z_1$ has the second largest magnitude i.e., $m_2({\bf Z}) = Z_1$. Let $P_2$ be the conditional probability of $m_2({\bf Z}) = Z_1$ and $E_1$ be the conditional mean of symbol ratio when $m_2({\bf Z}) = Z_1$.   3. $Z_1$ has neither the largest nor the second largest magnitude. Let $P_3$ be the conditional probability that this is the case and $E_3$ be the conditional mean when this occurs. Then, $E[R] = E_1P_1+E_2P_2+E_3P_3$. I was able to express each of these terms as repeated one-dimensional integrals with appropriate limits. **Is this the right approach?** The expressions for each of the terms are given below. $E_1 = \int\limits_{-\infty}^{\infty} \frac{1}{|z_1|} \frac{1}{\sqrt{2\pi \sigma^2}} e^{-\frac{(z_1-\mu)^2}{2 \sigma^2}} (K-1) \int\limits_{-|z_1|}^{|z_1|} \frac{|z_2|}{\sqrt{2\pi \sigma}} e^{-\frac{z_2^2}{2 \sigma^2}} \left[ \int\limits_{-|z_2|}^{|z_2|} \frac{1}{\sqrt{2\pi \sigma^2}} e^{-\frac{z^2}{2\sigma^2}} dz\right]^{K-2} dz_{2} dz_{1}$, where $\Phi(.)$ is the standard Gaussian distribution function. Using change of variables $x=\frac{z_1-\mu}{\sigma}$ and $y=\frac{z_2}{\sigma}$, we have $E_1 =\int\limits_{-\infty}^{\infty} \frac{1}{|x+\frac{\mu}{\sigma}|} \phi(x) (K-1) \int\limits_{-|x+\frac{\mu}{\sigma}|}^{|x+\frac{\mu}{\sigma}|} |y| \phi(y) \left[ 2\Phi( |y|) -1 \right]^{K-2} dy dx$ Similarly, $P_1 = (K-1) \int\limits_{-\infty}^{\infty} \phi(x) \int\limits_{-|x+\frac{\mu}{\sigma}|}^{|x+\frac{\mu}{\sigma}|} \phi(y) \left[ 2\Phi( |y|) -1 \right]^{K-2} dy dx $, $E_2=\int\limits_{-\infty}^{\infty} \frac{1}{|x|} (K-1) \phi(x) \int\limits_{-|x|-\frac{\mu}{\sigma}}^{|x|-\frac{\mu}{\sigma}} |y+\frac{\mu}{\sigma}| \phi(y) \left[ 2\Phi\left( \Bigg|y+\frac{\mu}{\sigma} \Bigg|\right) -1\right]^{K-2} dy dx$, $P_2 = (K-1)\int\limits_{-\infty}^{\infty} \phi(x) \int\limits_{-|x|-\frac{\mu}{\sigma}}^{|x|-\frac{\mu}{\sigma}} \phi(y) \left[ 2\Phi\left( \Bigg|y+\frac{\mu}{\sigma} \Bigg|\right) -1 \right]^{K-2} dy dx$, $E_3=(K-1)(K-2)\int\limits_{-\infty}^{\infty} \frac{1}{|x|} \phi(x) \int\limits_{-|x|}^{|x|} |y| \phi(y) \left[ \Phi\left( |y|-\frac{\mu}{\sigma} \right) -\Phi \left( -|y|-\frac{\mu}{\sigma} \right)\right] \left[ 2\Phi\left( |y|\right) -1\right]^{K-3}dydx$, and $P_3=(K-1)(K-2)\int\limits_{-\infty}^{\infty} \phi(x) \int\limits_{-|x|-\frac{\mu}{\sigma}}^{|x|-\frac{\mu}{\sigma}} \phi(y) \left[ \Phi\left( |y|-\frac{\mu}{\sigma} \right) -\Phi \left( -|y|-\frac{\mu}{\sigma} \right)\right] \left[ 2\Phi\left( |y|\right) -1\right]^{K-3} dy dx$. I was able to find the sample average by generating large number of random samples of $\bf Z$ according to the distribution. I am trying to figure out if numerical integration is feasible at all. On using matlab to evaluate the expressions for $E_1$ through $E_3$ and $P_1$ through $P_3$, the answers do not match with the sample average. This is the last puzzle in my research and your inputs will be highly appreciated.


title: Cumulative Hazard Function where "status" is dependent on "time"

text: This question is in relation to a dataset I've previously discussed here. I'm trying to determine whether a treatment affects not only the number of visits a patient has to the doctor, but also how long it takes for the injury to heal. The two are inherently interlinked. Now the structure of the cohort changes over time because of censoring; when a wound is healed, the patient doesn't have to see the doctor and those left in the cohort are those with longer healing times and thus more doctor visits. Because of this "intuitive" interpretation of the data, I thought a cox proportional hazards model (problem below) and/or "inverse" Kaplan-Meier Curves would be good to show how the initial treatment affects outcome. First I looked at the median and mean number of visits all patients had, which was around three. Then I stratified the whole cohort into those patients requiring $>3$ visits and $\leq 3$ visits. I then used the following function in R               library(survival)     km <- km(Surv(Time, Visits>3)~Treatment, data=mydata)     plot(km, fun="event")      This produced the following plot ![KM Plot](http://i.stack.imgur.com/Ou7Wt.png) I then wanted to basically do the same for a coxph hazards model, but have realized that its interpretation is a bit wonky as Treatment 2 was resulting in low HR, which getting your head around is a bit of task and I honestly don't think is right, because I'm trying to look at cumulative hazards. The code used in R was:               cox <- coxph(Surv(Time, Visits>3)~Treatment, data=mydata)     summary(cox)     cox(formula = Surv(Time, Visits>3) ~ Treatment, data=mydata)            n= 4302, number of events= 1514         (41 observations deleted due to missingness)                                             coef exp(coef) se(coef)      z Pr(>|z|)         Treatment 2                    -0.36705   0.69278  0.05318 -6.902 5.12e-12 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1                                         exp(coef) exp(-coef) lower .95 upper .95     Treatment 2                       0.6928      1.443    0.6242    0.7689          Concordance= 0.541  (se = 0.008 )     Rsquare= 0.011   (max possible= 0.99 )     Likelihood ratio test= 48.43  on 1 df,   p=3.419e-12     Wald test            = 47.64  on 1 df,   p=5.119e-12     Score (logrank) test = 48.13  on 1 df,   p=3.986e-12      So I'm wondering   * Is there a cumulative hazards function in R for this?    * Is it even ok to stratify patients into two groups based on a time-dependent variable?    * How would you try to actually interpret (in "lay-terms") the hazard ratio resulting from this model?    * Have I totally got this wrong?  Cheers for your opinions and help.


title: Repeated measures paired test

text: I have the following data:   2 conditions - experimental & control.   20 participants, each having 3 results for each condition. Like in the following table:               | Condition    | PARTICIPANT #1 | PARTICIPANT #2 | PARTICIPANT #3 |     |--------------|----------------|----------------|----------------|     | EXPERIMENTAL | Result #1      | Result #1      | Result #1      |     |              | Result #2      | Result #2      | Result #2      |     |              | Result #3      | Result #3      | Result #3      |     |--------------|----------------|----------------|----------------|     | CONTROL      | Result #1      | Result #1      | Result #1      |     |              | Result #2      | Result #2      | Result #2      |     |              | Result #3      | Result #3      | Result #3      |     |--------------|----------------|----------------|----------------|        I need to compare between the 2 conditions in a way that takes into account that the data is paired regarding participants (take into account that the same participant took part in both conditions). It seems to me like the best approach is a 2 factor repeated measures anova, but i would like to be sure and to know if there are other ways to analyze this type of data. Also, i would like to know which test should i choose if the order of the results (1-2-3) matters. I'll appreciate any input! Thank you!


title: Integration of Weka and/or RapidMiner into Informatica PowerCenter/Developer

text: Are there any existing integration packages (or similar) of Weka and/or RapidMiner into Informatica PowerCenter/Developer? I know there exists java- transformations in PowerCenter 9.x, so I could use the Weka API directly as code, but I just wondering if someone knows of existing projects like this so that I don't reinvent the wheel, or could stand on someone else's shoulders.


title: Calculation of relative information gain for multiple splits in decision trees

text: As it well known, we can calculate Relative Information gain (RIG) as follows: $RIG = \frac{H(x) - H(x|a)}{H(x)}$. In binary decision trees we calculate $H(x|a)$ for univariate split for variable $x_i$ as follows: $-(p(x_i>a)H(x|x_i>a) + p(x_i>a)H(x|x_i>a))$. So that means that for each point we look whether this point is from the left or from the right of current decision boundary. Now, my question is: is it possible to calculate $RIG$ for complex condition in binary tree nodes? Like $(x_i<3)\;or\;(x_i > 7)$? Lets say, that all points for which $x_i < 3$ will be in zone $I$, points that have $x_i$ will be in zone $II$, and every point with $x_i>7$ will be in zone $III$. Using this notation, is it mathematically correctly to calculate relative entropy in this case as follows: $-(p_{I}H_{I}+p_{II}H_{II}+p_{III}H_{III})$? Thank you.


title: Do significant control variables need to be eliminated from statistical analysis in general and SmartPLS specificly?

text: I have been told that we need to eliminate control factors that shows a significant impact on the dependent variable from the model. As I know that control variables are variables that have testified to have an impact on DV; therefore, we need to include them in the study. Thus, I have a confusion why we need to exclude them from the model in case they shows a significant impact on the DV; whereas they are included since they have high potential to affect the DV?


title: Estimating distribution of traffic across all Web pages

text: The web is estimated to have approx 5 billion web pages ( http://www.worldwidewebsize.com/ ). How would you go about estimating a CDF to know how many pages are required that make up, say, 90% of traffic? I'm trying to apply Zipf's law but I'm not sure how to do it. Of course, Zipf's law requires a parameter. I'd like to estimate that. As a starting point: What are typical values of that param in other domains? How consistent are they? Then: Given that certain data is available (eg logs of requests of small sets of users) - how can we use those to estimate the right Zipf param? Can we use PageRank data (which has lots published and can be calculated with enough computing power)? Can we use PageRank of a subset of the web? How would you, as a statistician, go about trying to get some answers to this? Put the question concretely: What % of web pages get 90% of the traffic?


title: spatial-location/time-series prediction models

text: How efficient it is to build predictive model. However, every crime is dependent on three factors: Time, Spatial Location and people behavior. Statistically, we can't measure people behavior (we can't use these private data). Hard to find good and accurate results using spatial-location/time- series predictive systems. Statistically, Is it valid to use only one dimension (spatial-location or time-series) to analyze and predict events (like crimes) ?


title: Logistic regression vs MLE of conditional probability

text: I am trying to build a classification model from a large dataset (50M instances) with a categorical independent variable (`pred`) and a binary dependent variable (`class`). The dependent variable is very skewed - only 2K positive instances. What confuses me is the very low predicted probability of logistic regression model. `(glm(class~pred,family=binomal)),` compared to standard conditional probability estimate               (prop.table(table(class,pred),2)[2,]):                          1            2            3            4            5            6      ...     7.898106e-03 5.133151e-03 3.332895e-03 2.162637e-03 1.402706e-03 9.095645e-04 (glm)...     1.029593e-01 1.572959e-02 4.906542e-03 1.695033e-03 7.622454e-04 5.922985e-04 (mle)...      Is this expected differences? Should I use some method other than logistic regression (results are similar with non-parametric locfit)? Is it OK to report conditional probabilities only?


title: How can I compute covariances between parameters in WinBUGS?

text: I am wondering if there is an easy way to compute the covariances between parameters in WinBUGS/OpenBUGS. It is easy to obtain the variances, but for subsequent analysis, I need the covariances. Consider the Salm example:               model     {         for( i in 1 : doses ) {             for( j in 1 : plates ) {                 y[i , j] ~ dpois(mu[i , j])                 log(mu[i , j]) <- alpha + beta * log(x[i] + 10) +                      gamma * x[i] + lambda[i , j]                 lambda[i , j] ~ dnorm(0.0, tau)                  cumulative.y[i , j] <- cumulative(y[i , j], y[i , j])             }         }         alpha ~ dnorm(0.0,1.0E-6)         beta ~ dnorm(0.0,1.0E-6)         gamma ~ dnorm(0.0,1.0E-6)         tau ~ dgamma(0.001, 0.001)         sigma <- 1 / sqrt(tau)     }      How would you obtain the covariance between alpha and beta? Alternately, is there anyway to save the last 1000 iterations and then manually compute the variances outside?


title: larger sample size or more accurate measurement?

text: This question is completely theoretical in the sense that I don't have any real data. Suppose I am designing a study, in which I will go measure N pairs of (X,Y), then the goal is to build a predictive model (say linear regression) between X and Y. i.e. I want good predictive accuracy, but it will also be good to have a prediction interval. Given the limited resources, I have two choices:   1. Increase N, but X will be measured less accurately; or    2. Decrease N but measure X more accurately. The question is, how can I find the balance between the 2 cases? I understand I haven't provided real data, which is why I am not asking for specific answer. I just want to know what's the line of thought you would follow? What's a good way to approach this question? etc.... **Clarification:** Don't know why I put it the other way around. X in the dot points above should be Y. I am predicting Y using X. My choice is between increasing N or measuring Y more precisely. **Added:** It will be good if your response considers both small sample size and adequate sample size situations.


title: How to find the level curves of a multivariate normal?

text: I'm trying to get some intuition for multivariate analysis. Thought it would be good to visualise, say, the .95 probability contour for a multivariate normal. I'm looking for suggestions on an efficient way to do this in R. My attempt was to generate a bunch of vectors $\boldsymbol v=(cos(x), sin(x))'$ for evenly spaced x in $[0,2\pi]$ and then try to manipulate the cdf to find what constant $k$ made $F(k*\boldsymbol v)=.95$ (when $\boldsymbol \mu = \boldsymbol0$ and $\boldsymbol \Sigma$ arbitrary). Couldn't get this working and wondered if there was a better approach. I didn't have much luck with the R package I tried - mvtnorm \- not sure what they mean by 'equicoordinate' quantile?


What is the sample standard error formula? I know only $s$ but I guess this is not it. I am confused about its formula. Please help me. Thank you.


title: Properties of bivariate standard normal and implied conditional probability in the Roy model

text: Sorry for the long title, but my problem is quite specific and hard to explain in one title. I am currently learning about the Roy Model (treatment effect analysis). There is one derivation step at my slides, which I do not understand. We calculate the expected outcome with treatment in the tretment group (dummy D is treatment or not treatment). This is written as \begin{align} E[Y_1|D=1] \end{align} since $Y_1=\mu_1 + U_1$ this can be rewritten as \begin{align} E[Y_1|D=1] &= E[\mu_1+U_1|D=1]\\\ &=\mu_1+ E[U_1|D=1] \end{align} before we also said, that $D=1$ if $Y_1>Y_0$ so it follows: $Y_1-Y_0>0$ $\mu_1+U_1-(\mu_0-U_0)>0$ $(\mu_1+U_1)/\sigma-(\mu_0-U_0)/ \sigma >0$ $Z-\epsilon>0$ so $D=1$ if $\epsilon<Z$ Therefore it holds, that \begin{align} E[Y_1|D=1] &=\mu_1 + E[U_1|\epsilon<Z] \end{align} It is further known, that \begin{align} \begin{bmatrix} U_1 \\\ U_0 \\\ \epsilon \end{bmatrix}=N\left( \begin{bmatrix} 0 \\\ 0 \\\ 0 \end{bmatrix}, \begin{bmatrix} \sigma_1^2 & \sigma_{10} & \sigma_{1\epsilon} \\\ \sigma_{10} & \sigma_{0}^2 & \sigma_{0\epsilon} \\\ \sigma_{1\epsilon} & \sigma_{0\epsilon} & \sigma_{\epsilon}^2 \end{bmatrix}\right) \end{align} therefore it follows: $P(D=1)=P(\epsilon<Z)=\Phi(Z)$ So now comes my question, the slides say, that \begin{align} \mu_1 - E[U_1|\epsilon<Z] =\mu_1 - \sigma_{1\epsilon} \frac{\phi(Z)}{\Phi(Z)} \end{align} And I do not understand why? I know, that if two random variables follow a standard bivariate normal distribution: $E[u_1|u_2)=\rho u_2$ so $E[u_1|u_2>c)=E[\rho u_2|u_2>c]=\rho E[u_2|u_2>c)=\rho\frac{\phi(c)}{1-\Phi(c)}$ Therefore I would have expected a "plus" and not a minus sign? Also why do we use the covariance $\sigma_{1\epsilon}$ and not the correlation $\rho$? So I would have expected something like \begin{align} \mu_1 - E[U_1|\epsilon<Z] =\mu_1 + \rho \frac{\phi(Z)}{\Phi(Z)} \end{align} I am aware of the fact, that if I do the truncation from above the $1-\Phi(c)$ becomes a $\Phi(c)$.


title: General theoretical properties of empirical Bayes estimates

text: I was wondering if someone could provide reference (if such exists) for the theoretical properties of empirical Bayes(EB) point estimates, in the sense of what can we say about their risk under certain regularity conditions, for general case.   For example, assume we have a prior which is structurally correct, but has an unknown hyperparameter. So using EB approach, we estimate the hyperparameter and plug it in, finally getting, let's say, the EB "MMSE" for the needed parameter. Can we say it'll always be better than just using the maximum likelihood (in the mse sense) ?


title: Estimating parameters and latent variables in a split Poisson distribution using R and JAGS

text: I am trying to estimate parameters and latent variables in a split Poisson model that describes observable and unobservable counts in time assuming the split probability is $\pi$. An observable event in time is denoted as $O_t$ while an unobservable event is denoted $U_t$. The distributions of $O_{t+1}$ and $U_{t+1}$ are assumed to follow the conditional distribution: $O_{t+1}|O_{t},U_{t}$ ~ Poisson$[\pi(\mu_{t}O_{t}+\mu U_{t})]$ $U_{t+1}|O_{t},U_{t}$ ~ Poisson$[(1-\pi)(\mu_{t}O_{t}+\mu U_{t})]$ $\mu$ and $\mu_r$ are type specific rates where $\mu_{t} = \mu$ if $t \le C$ and $\mu_t = \mu_r$ if $t \gt C$; $t=1,2,\ldots,C-1 , C, C+1, \ldots,N$. Data is available for all $O_t$. Each latent variable $U_t$ is to be updated using the relation above and assuming that $U_{t} = 0$. It is desired to estimate parameters $\pi$, $\mu$, $\mu_r$ and the latent variables $U_{t}$ for all $t$. Bayesian method is chosen to carry out these estimations and my code is given below:                   library(R2jags)         # sample observable data for 30 time points         o <- c(1, 1, 2, 11, 15, 30, 60, 46, 60, 54, 51, 44, 27, 20, 16, 10, 7, 7, 4, 4, 6, 10, 8, 10, 2, 4, 4, 1, 0, 0)         # setting the input for jags model function         oo <- data.frame(o)         datain <- as.list(oo)         datain$n <- nrow(oo)         datain$C <- C         model <- function() {             u[1] <- 0    #bringing in the data through loop and minding the rate conditions             for (i in 2:C) {                 o[i]~dpois(pi*(mu*o[i-1]+mu*u[i-1]))             }             for (i in (C+1):n){                 o[i]~dpois(pi*(mur*o[i-1]+mu*u[i-1]))             }             pi~dunif(0,1)     # setting non-informative priors             mu ~ dunif(0,10)             mur ~ dunif(0,1)                  for(j in 2:C){    # the priors for the latent variable are updated via the conditional relation                 u[j]~dpois((1-pi)*(mu*o[j-1]+mu*u[j-1]))             }             for(j in (C+1):n){                 u[j]~dpois((1-pi)*(mur*o[j-1]+mu*u[j-1]))             }              }              parameters <- c('pi','mu','mur','u[10]','u[15]','u[20]') # usually we obtain distribution for all latent variables but only decide to choose these three.         jagsfit <- jags(datain, model=model,parameters=parameters,                 n.chains=2,n.iter=5000000,n.burnin=100000,n.thin=5000)         jagsfit                       mu.vect     sd.vect         2.5%          25%          50%                  75%         u[10]    1.402007e+04    2037.291 1.103313e+04 1.213150e+04 1.397350e+04 1.561675e+04         u[15]    1.251555e+06  130601.609 1.055833e+06 1.132489e+06 1.251471e+06 1.352963e+06         u[20]    1.136881e+08 7039334.938 1.028741e+08 1.072936e+08 1.138424e+08 1.191725e+08         mu       2.466000e+00       0.021 2.431000e+00 2.449000e+00 2.465000e+00 2.485000e+00         mur      1.860000e-01       0.157 8.000000e-03 6.500000e-02 1.470000e-01 2.570000e-01         pi       0.000000e+00       0.000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00         deviance 1.512097e+04     143.825 1.487507e+04 1.499773e+04 1.510491e+04 1.524422e+04                         97.5%  Rhat n.eff         u[10]    1.771087e+04 0.998   200         u[15]    1.488133e+06 0.997   200         u[20]    1.262587e+08 0.998   200         mu       2.499000e+00 0.997   200         mur      6.220000e-01 1.004   200         pi       0.000000e+00 1.001   200         deviance 1.534904e+04 0.997   200      I simulated the observed data above by fixing parameters and writing function below:                   set.seed(123)         C =7;mu=2;mur =0.4;pi=0.7         o <- numeric(); u <- numeric()         o[1] <- 1; u[1] <- 0         f <- function(){             for(i in 2:30){                 if(i <= C){             o[i] <<- rpois(1,pi*(mu*o[i-1]+mu*u[i-1]))             u[i] <<- rpois(1,(1-pi)*(mu*o[i-1]+mu*u[i-1]))         }else{             o[i] <<- rpois(1,pi*(mur*o[i-1]+mu*u[i-1]))             u[i] <<- rpois(1,(1-pi)*(mur*o[i-1]+mu*u[i-1]))         }     }     return(list(o,u))     }          f()      In this case, the unobservable is:               u <- c(0, 1, 3, 0, 7, 12, 31, 32, 26, 20, 22, 15, 13, 8, 4, 2, 4, 4, 3, 4, 4, 4, 5, 5, 4, 2, 1, 0, 0, 0)         I want to recover the parameters and the latent variables back but the results I got were not OK and was thinking I have mis-specified some of the distributions. Can anyone please help look this problem? Thanks!


title: Violation of assumption of two sample Hotelling's test

text: I would like to compare the nutrient intake of men and women. But the assunption of same variance covariance matrix is violated and therefore two sample hotelling's test cannot be applied. How elsed can I compare the intakes?


title: Outlier detection using k-means in a binary classification problem

text: I'm using `k-means` in every class of a binary classification problem and remove samples that have high distance from center of my features (21 features so 21 dimensions problem) before inserting data set to a neural network. After designing neural network model, now i want use this model for a new data set (out sample). As you know we must use outlier detection parameters in per-process stage for out sample data (like normalization `x-min(x)/max(x)-min(x)` that will use `max(x)` and `min(x)` for normalization of out sample). what parameters of `k-means algorithm` should i use for out sample and how can I do that in `MATLAB` ? Thanks.


title: Using SVM coupled with RBF kernel for EEG data classification

text: I want to classify the significant features obtained from nonlinear methods using SVM-RBF kernel function. I have used 10-fold cross validation that divides 180 samples as training set and remaining 20 as testing set. Can you please guide me through the steps on how to implement SVM with RBF kernel to classify the data into two groups?


title: Epanechnikov kernel and GCV

text: Does anyone know how to select bandwidth with GCV for Epanechnikov kernel in R or Matlab? Thanks a lot!


title: How can 8 dimensions be reduced to 3?

text: I can't seem to understand how PCA works. My (lack of) math knowledge don't help me either. I have read that the new set of variables has to be a linear combination of the old set. What does that mean exactly? That there should be a way to multiply numbers a, b, c, ... with dimensions/variables x, y, z... and get the old bigger(!) set? Can you answer my example: if I have 8 variables/dimension can they be reduced to 3? Or a vector of 3 components (sorry for my lack of proper terminolgy, English is not my native language) is not a linear combination of 8 and thus no? Thanks.


title: How to test for temperature differences between years and sites?

text: I have a dataset with monthly temperatures for 4 different years and I would like to test if they are significantly different between the years. As I collected the data for different sites I wondered if there were some possibility to perform the calculation at once or if I have to repeat it for each site. I am using the statistical software R. Alternatively, can you indicate me some way to plot this? My data look like this (I have some missing values): Code:               > head(MeanTemp)       Site Year    Month     Temp     1    1 2008  January       NA     2    1 2008 February       NA     3    1 2008   March        NA     4    1 2008    April       NA     5    1 2008      May 16.22594     6    1 2008     June 18.53764      Thank you very much


title: Relatively normalizing values for collaborative filtering

text: I am trying to derive a formula for my collaborative algorithm problem to calculate popularity rating of an item. I am considering three factors to calculate rating for an item based on three different ratings $a$, $b$ and $c$, each with weight $w_1$, $w_2$ and $w_3$. The cumulative rating is: $$w_1 a + w_2 b + w_3 c$$ In this equation, I want to give $b$ less weight than $a$, i.e. $w_1 > w_2 > w_3$. But the problem is that how should I decide/calculate the weights such that they maintain this relative ordering when normalized (I am not looking for linear normalization) i.e. $$w_1 a > w_2 b > w_3 c$$ Is there a way to calculate $w_1$, $w_2$, $w_3$ dynamically and what normalization process I should use or is recommended? It would be really helpful to get comments on the approach as well.


title: Multinomial chi square with small expected values

text: I'm studying extinction in Austronesian languages, and am trying to find out if a subset of 384 languages is randomly selected with respect to extinction risk from a population of 1249 languages. "Extinction risk" can take 10 different values. But the first two expected values are very small (both 0.40) and make up most of the chi-square:                  Observed | 10   | 6    |  6   | 22    | 113   |  64    | 70    | 20    | 5     | 19        Expected | 0.40 | 0.40 | 5.63 | 10.86 | 59.54 | 129.54 | 90.11 | 25.74 | 12.07 | 13.68     Chi square  | 229.0| 77.9 | 0.0  | 11.4  | 48.0  |  33.1  |  4.5  |  1.3  |  4.1  |  2.1     contribution      ----------- | --------------------------------------------------------------------------      Chi square:  411.44          Critical value (alpha = 0.01): 21.67 with 9 degrees of freedom      I've heard that it is best if the expected values stay above 5 (and that one might consider Yates's correction if they aren't). What correction would you suggest?


title: Logistic regression for self-regulation

text: My dataset contains samples of the following variables:     * $X_0$: the state of the system at time 0 (a continuous scalar)     * $X_1$: the state of the system at time 1 (a continuous scalar)     * $Y$: some binary variable describing the system at time 1      * $\boldsymbol{C}$: a bunch of covariates I would like to investigate whether the system "auto-regulates" by investigating whether larger $X_0$ tends to lead to larger $Y$ which, in turn, leads to smaller $X_1$. If this is the case, then, through the variable $Y$, the state of the system tends to be regulated. The way I have been looking at it for now is through a logistic regression: $$ \text{logit } Y = \beta_0 + \beta_{x,0}X_0+\beta_{x,1}X_1+\boldsymbol{\beta_c}\boldsymbol{C} $$ I realized that whenever $\beta_{x,0}$ is positive and significantly different from zero, $\beta_{x,1}$ is negative and significantly different from zero (Wald test). Regardless of the details of the estimation (accounting for multi- collinearity, etc.), I am wondering whether the logistic regression is the right way to approach such a problem and, if not, what other types of model could I investigate.


title: Extracting fixed effects in System GMM (xtabond2), Stata?

text: In fixed effect, you can use the option "predict u" to see the fixed effects. But it is not possible in System GMM in xtabond2... Somebody can help me?


title: The main effect will be non-significant if the interaction is significant?

text: I am using linear mixed models to identify important factors, and it turns out that:   * `A`: significant    * `B`: not significant   * `A`×`B`: significant Does it mean that because `A`×`B` shows that the effect of `A` depends on the effect of `B`, only the effect of `A` is not actually significant? I have read many sources, and they seem to suggest that if the effect of `A`×`B` is significant, then we cannot interpret that the effect of `A` is significant on our dependent variable. Am I understanding right?


title: Is it possible that the interaction effect is significant while the individual effect(s) is not?

text: Consider a simple model: Y = A + B + A*B Is it possible that coefficient for A*B is significant while that for A or B or both of them is insignificant? Either yes or no, can you tell what's the reason behind it?


title: Calculate similarity of waiting times of users

text: Let's say I have waiting times(seconds) of users in web pages.               userid page waitingtime     1      p1   10     1      p2   20     1      p3   5      2      p1   2      2      p3   5     3      p1   10      3      p5   2     3      p6   5      I need to compare users based on the waiting time in each page. What measure can I use to calculate the similarity between users based on waiting time?


title: What is the type 1 error of rejecting samples with low p-values?

text: Consider an amateur statistician who will discard any sample which does not lead him to reject his null hypothesis at an alpha=.05 significance level with a two-sided test. What is this statisticians Type 1 Error? (He continues drawing independent samples until he is able to reject the null hypothesis.)


title: Median + MAD for skewed data

text: I am trying to figure out what happens if you apply Hampel's outlier detection technique based on the median and the MAD to data that is skewed. Apparently, the advantage of Hampel's method over z-scores is that it is much less influenced by outliers itself. However, several papers and websites say that this method should not be applied when the data distribution is skewed, so when the data is not normally distributed. However, I did not find any literature about what happens if you apply this method to skewed data. Does it not detected any outliers at all? Or does it detected false positives? I found several questions in this forum about whether using z-scores or Hampel's approach and even when data is skewed but no one gave an answer of what the outcome of Hampel's method is when it is applied to skewed data. The closest comment I found in this forum is the following: "Using the MAD amounts to assuming that the underlying distribution is symmetric (deviations above the median and below the median are considered equally). If your data is skewed this is clearly wrong: it will lead you to overestimating the true variability of your data." Mean$\pm$SD or Median$\pm$MAD to summarise a highly skewed variable? It says "it will lead you to overestimating the true variability of your data" but what does that actually mean? Does it lead to the identification of too many or too less outliers? In addition, can anyone see a problem of applying this technique to studies of small sample sizes compared to z-scores?? Can anyone help to shed light on that?


title: What is the mathematical model formula corresponding to this gam model fit in R?

text: If I have a generalised additive model how would I write this in the correct manner in a scientific paper? For example:               set.seed(10)     RandData <- rnorm(8760*2)     America <- rep(c('NewYork','Miami'),each=8760)          Date = seq(from=as.POSIXct("1991-01-01 00:00"),                 to=as.POSIXct("1991-12-31 23:00"), length=8760)          DatNew <- data.frame(Loc = America,                          Doy = as.numeric(format(Date,format = "%j")),                          Tod = as.numeric(format(Date,format = "%H")),                          Temp = RandData,                          DecTime = rep(seq(1, length(RandData)/2) / (length(RandData)/2),                                        2))     require(mgcv)     mod1 <- gam(Temp ~ Loc + s(Doy) + s(Doy,by = Loc) +       s(Tod) + s(Tod,by = Loc),data = DatNew, method = "ML")      So, if I wanted to write the Temperature of both locations can be described by the following model... How would I write this?


title: Inter-rater reliability: multiple raters and multiple replications

text: This is really a two part question. I have been asked to assess inter-rater agreement on a study in which 4 raters were asked to rate microarrays from patient samples. Ratings were based on a 5 point ordinal scale (0-4). I have experience with analyzing rater agreement for 2 raters but not for multiple raters so would appreciate insight as to the best way to go about that. I have done some research and have found two SAS macros that can be used for multi-rater agreement (MKAPPA and MAGREE). Would using one of these methods be the best approach? To make this problem a bit more complex, each sample was rated 3 times by each of the 4 raters (studies assessing cell lines often involve doing things in triplicates). How do you account for that? Do you just view each sample as independent? I’m not sure if this is the best approach, though, since they do want to compare all 12 ratings for one sample. Do you instead consider there to actually be 12 raters for each sample since each of the 4 raters rated each sample 3 times? I would appreciate any help with tackling this problem. If something here isn’t clear, please let me know. Thanks!


title: What model to use for a 'spraying experiment'?

text: I separate my plants into four groups of 16, and I spray them to keep the bugs away:               group 1: I spray daily     group 2: I spray every 2 days     group 3: I spray every 3 days     group 4: I spray every 4 days      at the end of the 12 days, I check number of plants that's bug free               group 1 is: 12     group 2 is: 10     group 3 is: 7     group 4 is: 6      How should I fit a model to this?


title: Help using maxFratio() in R (Hartley's test)

text: I found documentation at: http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/SuppDists/html/maxFratio.html I am not sure how really to use this though and am coming up short finding help online. If anyone has a link or tutorial on this in R please let me know.


title: How to sample a statistic?

text: Disclaimer: _I am a software developer and I like stats but I'm not a professional statistician, I already experienced that my wording is not always the correct jargon. Please keep in mind._ I would like your opinion on how to build a concrete statistical model: **The Goal** Why am I doing this: I need to assess the reliability of a food database. Is the food present that most of the people eat? With the help of _CrossValidated_ I made an arbitrary benchmark that represents what I want to know:   1. The top 30% food items    2. by amount a normal person eats in one meal   3. sold in Austria annually **The Target Group** The goal is to understand what are the most eaten foods in Austria (citizens: 8m), but more by People that are Internet savvy, therefore I tend to say urban, educated. There is a special segment _Mother with Children_ and _Athletes_ that I will address in a second study. **The Data** It is about branded food because the database is exhaustive with natural food. Each Brand has to be recorded separately. I would like to know what you think of my approach: I go to a supermarket and observe the cash desk for an amount of time and record every food that has been bought. This is one sample the population is the overall consumption of items per food. e.g. 23000 liters of coca-cola, 600000 apples   1. Is this a wise design?    2. What would be the appropriate sample size?   3. How would you do it? Take into account this is currently a hobby we are doing in 3 persons. I'm comfortable to spent a few hundred Euros on the inquiry but only if spent wisely.


title: How to show whether there's a significant higher rate?

text: Maybe someone can help me or at least give me some clues. I have 1057 patients all with different types of prostheses. One of them seems to have a higher revision rate than the others.               total numbers:     prosthesisA 662     prosthesisB 162     prosthesisC 151     other       82          revision rate:     prosthesisA 9   1,36%     prosthesisB 11  6,79%     prosthesisC 3   1,99%     other       4   4,88%      Could anyone tell me how I can show whether the revision rate of prothesisB (6,79%) is significantly higher or not compared to the other prostheses (incl. "other")? Thank you & kind regards Edit: Is it possible to compare only the prosthesisB with all others at once? I've created the following contingency table:                       prostB  other   total     Rev     11      16      27     noRev   162     868     1030             173     884     1057          Odds-Ratio = (11*868)/(16*162) = 3.68     p-Value (Fisher Test) = 0.001      Can I say now that the chances for a revision with prosthesisB are ~3.5 times higher compared to all others? Is this plausible since the p-Value of the Fisher's exact test is highly significant?


title: Kullback-Leibler divergence - interpretation

text: I have a question about the Kullback-Leibler divergence. Can someone explain why the "distance" between the blue density and the "red" density is smaller than the distance between the "green" curve and the "red" one? ![Graph of three pdfs](http://i.stack.imgur.com/Zrm76.jpg) Best regards, Marco


title: Optimizing the ending rule of gradient method

text: My intuition says that the third equation must be "the length of the gradient squared less than epsilon". $x_{k+1} = x_k - f(x_k)$   $x_{k+1} = x_k + 1$   $|f(x_k)|^2 < \epsilon$ However, I am not sure whether it is the standard form. **How would you write the standard form of the gradient method, and particularly its ending rule?**


title: How to calculate Sidak correction by hand?

text: I can calculate Sidak correction of the p value in Excel but is there a quick way to do it by hand? I am looking for something similar to calculating the margin of error by hand, i.e. dividing 1 with the square root of the sample size.


title: advantage of euclidean distance for classification

text: Has euclidean distance any advantage in compare to another distance based methods like Manhatan distance or Maximum difference metric?


title: Active learning using SVM Regression

text: I have trained an **SVM Regression model** using training data, $x_1,x_2,\dots,x_N$. I want to perform **active learning** to improve the model; i.e., I want to add more samples to the training data and relearn a better model, and to choose these new samples in such a way as to maximize the resulting model performance. For an SVM classifier, a useful heuristic for active learning is to choose samples that fall close to the decision boundary; i.e., for a particular sample, the 'confidence' ( _c_ ) can be computed using the SVM. Samples that have small |c| are more likely to improve the decision boundary in the retrained SVM. Any suggestions on how to do this for SVM regression? (I can generate samples at will, but it is costly to label them, so I want to know if I can use the 'already-trained' regression-SVM to help me decide which ones to label)


title: Plain english explanation of the Rayleigh distribution?

text: I need to understand the Rayleigh distribution for a homework assignment in computer networks. Unfortunately, I lack the background knowledge in the field of statistics and probability theory to understand the descriptions that are given elsewhere about the Rayleigh Distribution, as it leads me from one page to another and I get more confused. If you could use an example to illustrate the case it would be even more useful.


title: Suggestions for identifying key features

text: I have a large set of customer data. For these customers, I have devised a **_customer loyalty score_** which is a measure of the loyalty of the customer. I want to find the features that are strongly associated/correlated with this score. Features could be _number of purchases_ at various merchant types. One obvious answer would to be just to calculate the **_correlation_** for each feature with the customer loyalty score and see which have the highest correlations. Is this preferred way of doing this or are there better techniques?


title: Corrected multiple comparison: is the preceeding ANOVA necessary at all?

text: To compare multiple groups, one usually runs ANOVA. If ANOVA rejects the null hypothesis, one goes through the multiple comparisons. The multiple comparisons are necessary to establish which group is different. Multiple comparisons are performed applying Bonferroni correction (or more recent types of correction, such as Holm's or Hommel's). Such corrections guarantee the family wise error rate (FWER) to be at most alpha. FWER is the probability to return at least one wrong significance statement when comparing populations whose mean is actually the same. I wonder why is the preceeding ANOVA necessary at all. I think one should run directly the multiple comparison with adjusted p-values. In fact: a) if all means are equal: the multiple comparison will have the same FWER of ANOVA, because of the correction. b)Some means are different. Then multiple comparisons are anyway necessary. The preceeding ANOVA is often done, but I do not see a compelling reason for it. Thanks.


title: Multiplicative errors for linear model

text: I am trying to figure out the 'standard' way of handling multiplicative error in a linear model, i.e. my model reads: $$ Y_i = (ax_i + b)\varepsilon_i , \quad \varepsilon_i\sim\mathcal{N}(1, \sigma^2) $$ How do I fit this (in R)? My idea was to use a log transform on both sides: $$ Z_i = \log(Y_i) = \log(ax_i + b) + \log(\varepsilon_i) $$ The additive errors are now lognormally distributed, but how do I fit this? I would say this is a GLM with link function $g=\exp$ and lognormal errors, but R does not seem to know such a model. Am I missing something? This seems to be an extremely 'standard' model... Edit: I obviously mixed something up here ;) $\varepsilon_i$ should be lognormally distributed, then it is valid to take logs. But still, the link function would be the exponential - or not? So I tried               glm(z ~ x, data = data.frame(x, z = log(y)), family = gaussian(link = 'exp'))      but the link function does not exist...


title: Can someone shed light on linear vs. nonlinear mixed-effects?

text: I'm new here, so I hope this hasn't been covered already, but my first few searches didn't find anything. I am about to dive into learning R and my learning project will entail applying mixed- or random-effects regression to a dataset in order to develop a predictive equation. I share the concern of the writer in this post How to choose nlme or lme4 R library for mixed effects models? in wondering whether NLME or LME4 is the better package to familiarize myself with. A more basic (hopefully not dumb) question is: what's the difference between linear and nonlinear mixed-effects modeling? For background, I applied M-E modeling in my MS research (in MATLAB, not R), so I'm familiar with how fixed vs. random variables are treated. But I'm uncertain whether the work I did was considered linear or nonlinear M-E. Is it simply the functional form of the equation used, or something else?


title: PCA on original variables vs. PCA on residuals

text: I have 10 measures of snakes that are strongly correlated (body length, tail length and 8 measures of head size). My dataset consists of different snake sizes (but excluded non-adults) so smaller (younger) snakes have smaller heads and shorter tails (because they depend on body length). Usually there are two approaches how to implement PCA on measures in biology in this case. One of them is to regress all measures on body length and then preform PCA on residuals + body length. The second approach is to perform PCA on original variables. I don´t know why should I regress all measures on body length, because in PCA it is not problem, when my varables are correlated. Only sense I can see in this procedure is in reducing potential multicollinearity. If I always perform PCA using correlation matrix and oblique rotation in SPSS, does it make sense to do it on residuals? What is the difference between results of these two approaches? Aim of these PCAs is to reduce multidimensionality and then plot the graph of new components... Can one think, that performing PCA on residuals is something like "size- standardizing" and performing PCA on original variables is not "size- standardized PCA"? Does it make sense? Do I have to standardize all measures on the same size before performing PCA? I don´t think so, but I want to be sure. Thanks EDIT: some biological sources to this topic PCA on residuals: http://www.pnas.org/content/94/8/3828.full http://www.ispa.pt/ui/uie/pdf/OliveiraAlmada1995.pdf PCA on original variables: http://swfsc.noaa.gov/uploadedFiles/Divisions/PRD/Publications/Jefferson_VW2004(82).pdf http://www.zool.uzh.ch/static/research/oekologie/literatur/pdf05_01/2001Evolution.pdf These two approaches are most common in biology before PCA, few authors use for size-standardization ratio of dependent variable and body length. I don´t know if this question is unclear, or you don´t know answer so please, let me know how could I improve it.


title: How to compare 2 non-stationary time series to determine a correlation?

text: I have two data series that plot median age at death over time. Both series demonstrate an increased age at death over time, but one much lower than another. I want to determine if the increase in the age at death of the lower sample is significantly different than that of the upper sample. Here are the data, ordered by year (from 1972 through 2009 inclusive) as rounded to three decimal places:               Cohort A    70.257  70.424  70.650  70.938  71.207  71.263  71.467  71.763  71.982  72.270  72.617  72.798  72.964  73.397  73.518  73.606  73.905  74.343  74.330  74.565  74.558  74.813  74.773  75.178  75.406  75.708  75.900  76.152  76.312  76.558  76.796  77.057  77.125  77.328  77.431  77.656  77.884  77.983     Cohort B    5.139   8.261   6.094   12.353  11.974  11.364  12.639  11.667  14.286  12.794  12.250  14.079  17.917  16.250  17.321  18.182  17.500  20.000  18.824  21.522  21.500  21.167  21.818  22.895  23.214  24.167  26.250  24.375  27.143  24.500  23.676  25.179  24.861  26.875  27.143  27.045  28.500  29.318      Both series are non-stationary - how can I compare the two please? I'm using STATA. Any advice would be gratefully received. ![Data plots](http://i.stack.imgur.com/DMkM4.png)


title: Is centering a valid solution for multicollinearity?

text: Let's assume that $y = a + a_1x_1 + a_2x_2 + a_3x_3 + e$ where $x_1$ and $x_2$ both are indexes both range from $0-10$ where $0$ is the minimum and $10$ is the maximum. I found by applying VIF, CI and eigenvalues methods that $x_1$ and $x_2$ are collinear. Can these indexes be mean centered to solve the problem of multicollinearity?


title: one sample t-test: result is 'infinity' - what do I report?

text: I have been comparing the outcomes from a number of groups to an expected value. The outcomes are a simple scale with three possible options: -1, 0, +1 The expected value is the '0' and one of my groups has answered with all '1's with the result that my t value is reading 'infinity'. Do I report this as infinity? Or is there an alternative? I am also calculating effect sizes based upon the t-value which is obviously impossible in this case so what do I write instead? thanks


title: Optimal Number of Entries in Contest of Skill

text: **Objective:** I'm looking for the optimal number of unique entries in a contest of skill with monetary prizes. **Description:** The contests vary from as little as 20 entries up to 10,000+ entries. You can enter an unlimited amount of entries. Within each entry, there are 9 slots that can be filled with variables. Each slot has a unique set of variables that cannot be used in any other slot. The unique variables for each slot may number from about 2-30 depending on the contest. Once you enter your chosen variables in the slots, the variables score points based upon real-life events. Some slots traditionally score many more points than other slots. The entry with the highest amount of points wins. Prize money is split in the event of a tie score. **Assumptions:** I am an above average player. When facing similar competition, assume I can win 57% of the time (if playing in a fifty/fifty contest where half of the entrants win prize money). The contest has exactly 10,055 entries. The entry fee is $2. The prize distribution is detailed as follows (minus 10.5% entry fee going to contest organizers): **Prizes:** **1st:** $2,000; **2nd:** 1,200; **3rd:** 800; **4th:** 500; **5th:** 350; **6th:** 300; **7th:** 250; **8th:** 200; **9th:** 180; **10th** : 160; **11th:** 150; **12th:** 140; **13th:** 130; **14th:** 120; **15th:** 110; **16th:** 100; **17th - 18th:** 90; **19th - 20th:** 80; **21st - 22nd:** 70; **23rd - 25th:** 60; **26th - 28th:** 50; **29th - 35th:** 40; **36th - 50th:** 30; **51st - 60th:** 20; **61st - 75th:** 18; **76th - 100th:** 16; **101st - 200th:** 14; **201st - 400th:** 12; **401st - 600th:** 10; **601st - 800th:** 8; **801st - 1050th:** 6 **Discussion:** I have an above average understanding of probability (but that's not really saying much), so I've gotten as far as running a simulation (the paper and pencil kine) with 10 entries and a 10 dollar entry fee (1st-$50, 2nd-30, 3rd-20). Internet research suggests, and my paper simulations concur, that using 2 or more of the exact same entries will always result in a negative expected value (EV). I then tried to run the simulation with two different entries. There seem to be 44 unique possibilities. I added up the prize money; then, I multiplied by 1/44. I took this answer then subtracted the $20 in entry fees. The result was an expected value of 45 cents. Is this correct? Assuming my math is correct, I have no clue how to apply it to the contests with a much higher number of entries to find the optimal number of unique entries I should submit. I don't really know how to translate my assumed 57% success rate in one contest to another contest. I also don't know if it's best to have completely unique entries or just slightly unique entries (it seems that I should diversify my entries and make them as unique as possible while still maximizing points scored). Any help, guidance, or tips are very much appreciated. Have a great day. P.S. If you know the contests I'm talking about, please don't refer to them by name. P.S.S. I wasn't sure whether to post this question here or on the math.stackexchange


title: How is Support Vector Machine related to nonlinear programming?

text: I have read in few papers that that Support Vector Machine (SVM) is a NP hard problem, I cannot relate these two together. I was wondering if somebody can explain the relation between SVM and nonlinear programming.


title: Is Outlier detection in two separate databases is equal to one combined database?

text: Suppose that we have two databases : `Database_1` and `Database_2` . `Database_1` has 300 samples and `Database_2` has 700 samples. `Database_all` is combination of two databases.   * Is finding outliers using `abs(X-mean(X))>=1.9*std(X)` in `Database_all` is equal to finding outliers separately in `Database_1` and `Database_2`? In both cases we will remove same samples or not?   * Suppose that we will remove outliers with column 1 on above databases. In column `2` and `3` we have data that we want create a label for every sample like this `(column_2_value-column_3_value)/mean(column_2)` (`mean(column_2)` is average of all values in column `2`) We want calculate this label for all samples after removing outliers. In this case we have same label values if we remove outliers from `Database_all` compare to calculate after removing outliers separately in `Database_89` and `Database_90`? Please mathematically prove this problem if you have time. Thanks.


title: simply demostration on conditioned probability

text: I don't find the answer of this simply problem: how can I algebrically demonstrate that, with 3 variable A,B,C $ \sum\limits_C P(B|C)P(C|A)=P(B|A)$ under condition that $P(A,B,C)=P(A|C)P(B|C)P(C)$ it seems trivial, but I miss the algebric steps.... thanks


title: OLS standard error log log regression

text: I am estimating the following Power Law relationship: $$\ln(\text{Rank}) = \text{constant} + \alpha \ln(\text{Size})$$ where $\text{Rank}$ is $1,~2,~3,~...,~n$, and $\text{Size}$ is the raw value. What are the standard errors of the alpha coefficient in this instance? I read in an article that you cannot use the normal formula because Rank is autocorrelated (by construction) and taking logs alters some properties. If it helps, the formula I have seen for the standard error in this instance is $\text{SE}(\alpha) = \alpha (2/n)^{\frac{1}{2}}$.


title: Formula for calculating sample size for hypergeometric distribution

text: Suppose you have an urn with "N" marbles. All marbles are either black or white. You take a sample of size "n" without replacing them to the urn. With this one sample you would like to be able to make one of the following two statements:   1. Y% or greater of the marbles are white.   2. Less than Y% of the marbles are white. **Thoughts:**   * This approach uses the hypergeometric distribution.   * In practice, my "N" will be large (1*10E5 to 1*10E6). **Questions:**   1. What is the size of "n"?    2. What is the formula for calculating the required sample size "n"?    3. How to estimate the confidence interval for Y? Thanks in advance.


title: Blackwell-MacQueen Urn Scheme

text: I'm having some difficulty understanding the steps in the proof in the bottom half of the second page of this paper: http://www.cs.berkeley.edu/~jordan/courses/281B-spring04/readings/blackwell- macqueen.pdf. Specifically, let: \begin{equation} m_n(A)=\frac{\alpha}{\alpha+n}\mu(A)+\frac{1}{\alpha+n}\sum_{i=1}^n{\mathbb{I}(X_i\in A)}, \end{equation} where $A$ is any subset of the probability space of $\mu$, and $\mu$ is any probability "base" measure. We call a sequence of random variables $\\{X_n\\}_{n=1}^{\infty}$ a Polya sequence with parameters $\alpha>0$ and $\mu$ being any distribution, we have: \begin{equation} \mathbb{P}\left(X_1\in A\right)=\mu(A), \end{equation} \begin{equation} \mathbb{P}\left(X_{n+1}\in A|X_1,\dots,X_n\right)=m_n(A). \end{equation} Also, let $I_j$ be the indicator of the event that $X_j$ is different from all $X_i$ with $i<j$ and define: \begin{equation} \begin{split} &f_{nj}=I_jm_n(X_j),&\text{for }1\leq j\leq n,\\\ &f_{nj}=0,&\text{otherwise.} \end{split} \end{equation} We want to prove that: \begin{equation} \sum_{\substack{1\leq j\leq r\\\x_j\in A}}{f_{nj}}\leq m_r(A), \end{equation} and that: \begin{equation} m_r(A)\leq \sum_{\substack{1\leq j\leq r\\\x_j\in A}}{f_{nj}}+1-\sum_{1\leq j\leq r}{f_{nj}}. \end{equation} I would really appreciate your help! Thanks,   Anthony


title: Help understanding regression models with dlm in R

text: In the "Dynamic Linear Models with R" book, the Regression models section reads: "The static regression linear model corresponds to the case where $W_t = 0$ for any $t$, so that $\theta_t = \theta$ is constant over time." I am not understanding what this means, because when I fit a `dlmModReg` with `dW = 0` (default values) and then plot the predicted values of the state vectors, my regression coefficients vary over time, and eventually stabilize to a value similar to that of standard least-squares regression. Plot of regression slope coefficient with `dW = 0`: ![dW = 0](http://i.stack.imgur.com/I5MLB.jpg) However, if I use `dlmMLE` to find the MLE of `dW` before fitting the model, plotting the predicted values of the state vectors results in non-sensible values with multiple large discontinuities. Plot of regression slope coefficient with `dW` MLE and no intercept: ![dW != 0 no intercept](http://i.stack.imgur.com/VzdIB.jpg) Plot of regression slope with `dW` MLE with intercept: ![dw != 0 with intercept](http://i.stack.imgur.com/2Yu7P.jpg) I've yet to find any literature on how the intercept coefficient gets set, but I'm observing a near perfect linear relationship between the intercept and slope coefficient. The inclusion of an intercept term in the parameter MLE also causes my `dV` to drop from 16 to 0.003. Can anyone point me to any references which discuss how the intercept term is set on each update? Plot of intercept and slope with `dW` MLE: ![intercept and slope](http://i.stack.imgur.com/TjvTW.jpg) **My questions are:**   1. Why are the regression coefficients dynamic/time-varying even with `dW = 0` in the first example? If the regression coefficient changes at every time step regardless of whether `dW` is 0 or not due to the measurement update stage $\beta_t = \hat{\beta_t} + K_t(y_t - \alpha_t - \hat{\beta_t}x_t)$ ($K$ is the Kalman gain), what's the difference between simple linear regression and dynamic linear regression except for some random variation added in the time update equation $\hat{P_t} = P_{t-1} + dW$ ($P$ is estimate error covariance)?   2. Why does the plot of regression coefficients have so many discontinuities when my `dW != 0`?   3. What is the relationship between the intercept and the slope coefficients? Plotting them reveals a near-perfect linear relationship, but I can't find any literature explaining this. I haven't found the intercept term included in any formulation of the Kalman update equations. Would anyone be willing to take a look at my data/code?


title: Combining two covariance matrices

text: I'm calculating the covariance of a distribution in parallel and I need to combine the distributed results into on singular Gaussian. How do I combine the two? Linearly interpolating between the two almost works, if they are similarly distributed and sized. Wikipedia provides a forumla at the bottom for combination but it doesn't seem right; two identically distributed distributions should have the same covariance, but the formula at the bottom of the page doubles the covariance. Is there a way to combine two matrices?


title: Unpaired t-Test with equal sample sizes and unequal variances

text: I have two samples of $n=10$ with the following values Sample 1:   * Mean = $3$    * $s_{d} = 0.4$ Sample 2:   * Mean = $3.35$    * $s_{d} = 0.3$ but the observations of both samples are unknown and I want to know the common $s_{d}$ to calculate $t_{0}$ * * * Note: this is no homework is only an exercise from the tutorial, which has not been solved.


title: getting the variance of a set composed of a group of subsets

text: I have a data set $\mathcal{S}$, which is composed of several subset as $\mathcal{S}=\\{S_1,...S_i,...S_n\\}$ Each $S_i$ has its own mean $\mu_i$, variance $\sigma_i$, and size $\|S_i\|$. Is there a way to get the variance of the whole data set $\mathcal{S}$.


title: How do I find peaks in a dataset?

text: If I have a data set that produces a graph such as the following, how would I algorithmically determine the x-values of the peaks shown (in this case three of them): ![enter image description here](http://i.stack.imgur.com/tyLDG.png)


title: What to do with very low Durbin-Watson?

text: For 100 companies, I have collected (i) `tweets` and (ii) corporate website `pageviews` for `148` days. The tweetvolume and pageviews per day are two independent variables corpaired against the stock `trading volume` for each company, resulting in 100 x 148 = 14,800 observations. My data is structured like this:               company  date  tweetVol  pageviewVol  tradingVol     ------------------------------------------------     1        1     200        150          2423325     1        2     194        152          2455343     1        3     214        199          3100429     .        .      .          .              .     .        .      .          .              .     1       148    205        233          2563463     2        1     752        932          7434124     2        2     932       2423          7464354     2        3     600       1435          5324323     .        .      .          .              .     .        .      .          .              .     .        .      .          .              .     100      148     3         155           32324      Because there is much difference in company-size (some companies only receive 2 tweets per day, where others like Apple get over 10,000 per day), all variables are logged to smoothen distribution. (This is in line with previous research - this is for my thesis). I just performed a linear regression on this data, including both independend variables. R-Squared is .411 but Durbin-Watson only .141 (!) Without looking for the exact bounderies, I know this directly means my residuals are non- linear, eg. auto-correlated, right? My question is: how can I solve this? When I think about it, this data should not be autocorrelated, so I don't really understand. Is it due to this actually being a timeseries analysis? I wouldn't think that either, since for instance trading volume today is independent of yesterdays trading volume. Can somebody explain this to me? P.S. At my university, we use SPSS/PASW without additional modules, so I am unable to perform a timeseries analysis on this like you could in STATA or R.


title: Real World Challenge: Large difference between training and testing set accuracy

text: I have a classification dataset of ~100,000 rows and ~200 features. Within the dataset my predictor variable (Y) is an integer value between 0-55, therefore I am trying to predict 1 of 56 possible classes. I split my training/testing set into 80/20% and performed an extensive 10-fold cross validation exercise to tune the parameters and fit the final model. I end up a very high accuracy (and F1) score in the **training set (~90%)** but a very slow accuracy (and F1) score in the **testing set (~10%)**. **EDIT 1:** The training and testing set are split randomly. About 150 of the features are binary features (0 or 1) and the others are continuous values which I center and scale to be between 0 and 1. I have tried numerous learning algorithms (SVM, NN, logistic regression, PCA + SVM) and believe that the 10-fold CV should have eliminated overfitting as best as possible. However, nothing that I try seems to yield any meaningfully different results. Can anyone suggest new ways to increase accuracy of the testing set? **Caveats:** 1) This is real world data so getting more of it is very expensive and time consuming. 2) We need these 56 classes so cannot simply eliminate any. Thanks and any suggestions are appreciated.


title: How do you test for trend in a generalized linear model, particularly for odd numbers of levels of a variable?

text: I'm trying to test for a trend across levels of a categorical variable using a contrast statement in a model and can not figure out how to do it when there are odd numbers of levels. Looking at SAS's support page for contrast statements (ie http://support.sas.com/kb/22/912.html) it's easy to set up for even numbers of levels using coefficients balanced around 0, but when you have an odd number the value in the middle would be assigned the coefficient 0, meaning it would not be counted in the contrast at all. I've contacted SAS support directly and they suggested I ask on a statistical forum such as this.


title: LASSO in R for variable selection: how to choose the tuning parameter

text: I apologize in advance if this question is basic. I am trying to use LASSO for variable selection, with an implementation in R. I currently have 15 predictors, and looking to reduce the variable space and select the best predictors only, to be included in my final factor model. Some have advised me to use LASSO for this purpose. However, after reading some documentation about the subject, I am still unsure of how to choose the **tuning parameter** $\lambda$. As for the implementation in R, I attempted to use the `glmnet` package: **Get value of tuning parameter $\lambda$**               library(glmnet)     ans <- cv.glmnet(data, return[,1], standardize = TRUE)     par(mfrow=c(1, 2))     plot(ans$glmnet.fit, "norm", label=TRUE)     plot(ans$glmnet.fit, "lambda", label=TRUE)     dev.new()     plot(ans)      However, I am not sure how to interpret the results. I then run the following:               model <- cv.glmnet(data, return[,1], standardize = TRUE)     coef(model)      Could anyone please clarify the results of LASSO in R? Thanks,


title: What problems do non-normality in predictor variables cause for a multiple regression analysis?

text: I am talking about a situation in which I have several continuous predictor variables predicting a continuous outcome. One of the predictors has a very non-normal distribution and has some wild outliers. I intend the generalize the regression model to a wider population. Does the non-normality and/or the presence of outliers preclude my analysis, or cause difficulties in interpretation? What should I do?


title: can Dirichlet prior distribution be larger than 1?

text: This question is related to my quest of clustering the sequences using mixture Markov modeling. I have trouble understanding Dirichlet priors in the context of MAP-estimate (Mixture Markov Models). Namely, my priors end up being (much) larger than one.   I have non-informative priors defined as following: $$ p(\theta_n^{j}|a_n^{j})=\frac{\Gamma(\sum_{m=1}^{M}(a_{nm}^{j}+1))}{\prod_{m=1}^{M}\Gamma(a_{nm}^{j}+1)}*\prod_{m=1}^{M}(\theta_{nm}^{j})^{a_{nm}^{j}}, $$ where each $a_n^{j}$ is a _M_ -vector with components $a_{nm}^{j}>0$. _j_ denotes the jth component in the mixture and _n_ is the number of the row of the TPM ( transition probability matrix). Then I use the sum of log(Dirichlet Priors) across each row and each component. The most confusing aspects are: 1) In all literature the Dirichlet prior formula is given as: $$ p(\theta_n^{j}|a_n^{j})=\frac{\Gamma(\sum_{m=1}^{M}(a_{nm}^{j}))}{\prod_{m=1}^{M}\Gamma(a_{nm}^{j})}*\prod_{m=1}^{M}(\theta_{nm}^{j})^{a_{nm}^{j}-1}, $$ (notice the -1 in the exponent term). Is there possibly a typo in the article or can the -1 term be skipped? 2) The article sets $a_n^{j}$ equal to 10% of the corresponding relative frequencies of the TPM of the original counts across all sequences ( if our model has just 1 component). Then, let us consider such example: there are 3 possible states in the Markov Chain, and the _n_ -th row transition probabilities are 0.1, 0.8, 0.1. Let $a_n^j$ be equal to (0.01,0.03,0.06). Then, following the formula, my prior will be ( calculated in R) : according to the first version of formula: **1.961152** according to the second version of the formula: **245.144** This has been computed by defining a function in R:                   > dirichletPrior<-function(matrix_row,alpha_row){             dirichl<-((gamma(sum(alpha_row+1)))/(prod(gamma(alpha_row+1))))*prod(matrix_row^(alpha_row))             dirichl          }      The results seem a bit like a nonsense to me, and I do not understand where is my logic faulty. 3) Is "Dirichlet prior" a probability density function or a likelihood function? What would it mean if the result for each row is larger than 1? Are alpha multinomial parameters like I gave in the example, meaningful at all? What do they have in common with concentration parameter? The article I have been referring to is located under the following link: http://www.cs.uoi.gr/~kblekas/papers/C19.pdf.


title: Compare influence of same set of independet variables on two different dependet variables

text: I'm currently doing two multiple linear regressions. Each of them with the same set of predictors (measurements for real estate quality) $X_1,...,X_n$, but with different dependent variables (one of them the purchasing price, the other one the yearly rent) $Y_1$ and $Y_2$. $Y_1= a_1X_1+a_2X_2+...$   $Y_2= b_1X_1+b_2X_2+...$ What I am interested in is the influence of the independent variables on a third dependent variable $Y_3$ (a real estate investor's return assumption), which is approximately the quotient of the first two dependent variables $Y_1/Y_2$. So what I want to do is to find out which of the independent variables could possibly have an influence on the third dependent variable. As predictors for the regression on $Y_3$ I only want to use the independent variables out of the original set, of which I think that they have influence on $Y_3$. To do so I want to compare the influence of the independent variables on $Y_1$ and $Y_2$. If the influence points in a different direction (e.g. coefficient $a_1$ is negative and $b_1$ is positive) it is obvious that this independent variable will probably have influence on $Y_3$. But what if the direction is the same? It could happen, that the independent variable $X_n$ has strong influence on both $Y_1$ and $Y_2$ but with the same magnitude, so that $Y_3$ is not determined by this independent variable. _So my question is_ , is there a way to find out (e.g. by comparing standardized regression coefficients?), how large the influence of one predictor is relatively on $Y_1$ and $Y_2$? So I can say for example, $X_n$ determines both $Y_1$ and $Y_2$ in a positive way, but $Y_1$ is determined stronger, so that the quotient $Y_1/Y_2$ and therefore probably $Y_3$ is influenced by $X_n$ in a positive way, that's why I use it as a predictor in the regression on $Y_3$." I don't want to use the quotient of the predicted $Y_1$ and $Y_2$ as a estimator for $Y_3$, but do a fully new regression on $Y_3$.


title: Weights argument in glmer(), proportion data

text: What do the `weights` argument in `glmer` refer to? I used sample sizes as `weights` with `glm`, but here I am not sure. The variance of sample sizes is quite low, but including it or not in `glmer` gives me a huge difference. For example, in the dataset below, using only one independent variable, the difference in results is huge (estimate, BIC, p.value). Does anyone have experience using weights in `glmer` and confirm if it works as expected or if I am doing it right? A thread on r-sig-mixed models and in github pages, there seems to be an issue with `weights` argument in `glmer`, but since my knowledge of mixed-models is only weeks old, I am not able to follow it. My data:               cv <- structure(list(name = c("AlfF", "AndH", "AntH", "BerG", "BerR","FreZ", "GerB","GerT", "GueS", "GueV", "HanN", "HeiW", "JakW","KarN", "KerG", "KlaS", "ManS", "MarS", "SilN", "TheG", "UweP","WerT", "AlfF", "AndH", "AntH", "BerG", "BerR", "FreZ", "GerB","GerT", "GueS", "GueV", "HanN", "HeiW", "JakW", "KarN", "KerG","KlaS", "ManS", "MarS", "SilN", "TheG", "UweP", "WerT", "AlfF","AndH", "AntH", "BerG", "BerR", "FreZ", "GerB", "GerT", "GueS","GueV", "HanN", "HeiW", "JakW", "KarN", "KerG", "KlaS", "ManS","SilN", "TheG", "UweP", "WerT", "AlfF", "AndH", "AntH", "BerG","BerR","FreZ","GerB", "GerT", "GueS", "GueV", "HanN", "HeiW","JakW", "KarN", "KerG", "KlaS", "ManS", "SilN", "TheG", "UweP","WerT", "AlfF", "AndH", "AntH", "BerK", "BerR", "ChrG", "FraR","FreZ", "GerB", "GerB", "GerT", "GueS", "GueV", "HanN", "HeiW","JakW", "KlaS", "ManS", "MarH", "PetS", "SilN", "TheG", "UweP","WerT", "AlfF", "AndH", "BerK", "BerR", "ChrG", "FraR", "FreZ","GerB", "GerB", "GerT", "GueV", "HanN", "HeiW", "JakW", "KlaS","ManS", "MarH", "PetS", "PetW", "SilN", "SveR", "UweP", "WerT","AlfF", "AndH", "AntH", "BerK", "BerR", "ChrG", "FraR", "FreZ","GerB", "GerB", "GerT", "GueS", "GueV", "HanN", "HeiW", "JakW","KlaS", "ManS", "MarH", "MicH", "PetS", "SilN", "SveR", "UweP","WerT"), prop_yes = c(0, 0.2, 0.6, 0.1, 0, 0, 0.1, 0, 0.3, 0,0, 0, 0, 0.1, 0.8, 0.1, 0.1, 0, 0, 0, 0, 0, 0, 0, 0.3, 0, 0.3,0.3, 0.5, 0.4, 0.778, 0, 0.05, 0.5, 0, 0, 0.4, 0.2, 0, 0, 0, 0.2, 0.2, 0.3, 0.2, 0.6, 0.2, 0.2, 0.1, 0.1, 0.1, 0, 0.1, 0.3, 0.4, 0.1, 0.111, 0, 0.2, 0.1, 0.2, 0.8, 0, 0.111, 0, 0.1, 0, 0.2, 0.3, 0.1, 0.4, 0.333, 0.2, 0.1, 0.2, 0, 0.2, 0.182, 0, 0.1, 0.364, 0.1, 0.3, 0.375, 0, 0, 0, 0.2, 0, 0.1, 0, 0, 0, 0, 0.1, 0.1, 0, 0.3, 0, 0, 0.3, 0, 0.333, 0, 0, 0.667, 0.2, 0.571, 0.2, 0, 0.2, 0.6, 0.2, 0, 0, 0, 0, 0, 0.2, 0, 0, 0, 0, 0.2, 0.3, 0, 0.7, 0.3, 0, 0.2, 0.75, 0.2, 0.1, 0.1, 0.4, 0.1, 0.4, 0.3, 0.222, 0.2, 0.1, 0.1, 0.5, 0.2, 0.6, 0, 0, 0.1, 0.167, 0.333, 0, 0.222, 0.4, 0.5, 0, 0.3, 0.1, 0), size = c(10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 19L, 10L, 6L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 9L, 10L, 6L, 10L, 10L, 10L, 10L, 10L, 10L, 9L, 10L, 20L, 10L, 8L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 9L, 10L, 10L, 10L, 10L, 10L, 11L, 10L, 10L, 11L, 10L, 10L, 8L, 7L, 10L, 10L, 10L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 9L, 10L, 10L, 4L, 6L, 10L, 9L, 9L, 10L, 7L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 8L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 9L, 10L, 10L, 10L, 10L, 5L, 10L, 10L, 10L, 10L, 12L, 12L, 10L, 9L, 10L, 10L, 9L, 10L, 10L, 10L), tmean_winter = c(-3.83, -4.31, -3.97, -5.21, -4.6, -4.09, -4.05, -4.09, -4.85, -4.48, -4.77, -6.66, -4.16, -4.68, -4.48, -5.07, -3.83, -4.28, -4.79,  -4.83, -4.09, -4.43, 2.36, 1.47, 2.13, 1.09, 1.93, 2.26, 2.28, 1.98, 1.66, 1.3, 1.69, -1.01, 2.22, 1.89, 2, 1.23, 2.33, 2.1, 1.68, 1.66, 1.95, 1.38, 1.61, 0.86, 1.82, 0.48, 1.45, 1.74, 1.5, 1.78, 1.14, 0.65, 1.17, -1.59, 1.69, 1.55, 1.44, 0.65, 1.59, 1.16, 1.14, 1.23, 0.81, -1.53, -2.61, -1.52, -2.7, -1.77, -1.54, -1.68, -1.32, -2.16, -2.82, -1.95, -4.56, -1.57, -1.77, -1.76, -2.55, -1.51, -1.98, -2.05, -1.97, -2.62, -4.48, -5.25, -4.04, -4.92, -4.59, -5.34, -5.09, -4.12, -4.36, -5.23, -4.94, -4.7, -5.28, -4.55, -7.07, -4.18, -5.17, -4.56, -4.56, -4.74, -4.58, -4.62, -5.08, -5.25, -1.87, -2.67, -2.84, -2.47, -3.11, -2.3, -2.01, -2.05, -2.96, -2.57, -2.75, -2.54, -4.18, -2.07, -3.04, -1.81, -2.39, -2.24, -2.75, -2.75, -2.79, -2.44, -2.85, -0.35, -1.47, -0.43, -1.02, -0.76, -1.23, -1.57, -0.48, -0.65, -1.18, -0.58, -0.92, -1.58, -1.07, -4.05, -0.52, -2.1, -0.36, -0.75, -1.04, -0.67, -1.05, -1.81, -0.61, -1.64)), .Names = c("name", "prop_yes", "size", "tmean_winter"), row.names = c(NA, -158L), class = "data.frame")          glmer(prop_yes~tmean_winter+(1|name), family='binomial', data=cv)     glmer(prop_yes~tmean_winter+(1|name), family='binomial', data=cv, weights=size)


title: Identifying below-average units in a Poisson model

text: I’m modelling tuberculosis (TB) case rates at the neighborhood level and trying to identify risk factors that are associated with higher rates. I would like to identify communities with below average case rates given their risk factor profile with the idea that these would be the best places for new case- finding. It’s a bit circular I realize to use the same units I’m identifying risk factors with to then identify areas with expected higher rates, but my thinking is that I could identify the communities which are at the lower end of the distribution for each risk factor. My plan then would be to construct a Poisson model with significant risk factors, plug in each neighborhood and calculate the ‘expected rate’. Then calculate the difference by subtracting ‘observed’ from ‘expected’ rates, and then identify the communities with the highest difference. Is this approach legitimate? Or is this using the model in a way it isn’t meant to be used?


title: Identifiability of normal distribution

text: I am working on an exercise problem and am stuck in this problem: Suppose that $X_1,\dots,X_n$ are independent with $X_i\sim\mathrm{N}(\alpha_i + \nu, \sigma^2)$. Let $\theta = (\alpha_1, . . . , \alpha_n, \nu, \sigma^2)$ and consider the family of sampling distributions $P_θ$. Is the above parameterization identifiable?


title: In R, how to sample from the output of combn(a,b) if the "a choose b" is too large?

text: I would like to sample 1 million vectors from the set of vectors $$\mathbf{v} = \\{v_1,v_2,...,v_n \\} $$ subject to   * $\sum_{i=1}^n v_i = 1$   * $\forall i: v_i \ge 0$   * each $v_i$ has at most 2 decimal places As an example I have $n = 10$. Effectively, I am choosing to divide $100$ items into $10$ portion such that there is no minimum number of items that can be in any portion (meaning $0$ item in a portion is possible). This is an interesting combinatorics (albeit elementary) problem. Suppose a rich man wants to divide his \$100 billion wealth amongst his 10 sons. Suppose he wants each son to get at least \$1 billion and the *i*th son should receive $n$ billion where $n$ is an integer greater than 0. Imagine he has stacked his money in \$1 billion piles in a row with some gap between each pile. Then he would have 99 gaps between the piles. To divide his money he has to simply choose 9 cut points from the 99 gaps. So it's 99 choose 9. It general if he has \$c billion and d sons, it would be c-1 choose d ways to divide his wealth. Now suppose he wants to divide his wealth so that there is no minimum and the *i*th son could get $n_i$ billion where $n_i$ is an integer greater than or equal to 0. How many ways are there to do this? He was clever and asked his friend to lend him \$10 billion and he also stacked them into a row together with his other billions. Now there are 110 piles of \$1 billion dollars in a row and there are 109 gaps. He chooses 9 cut points from these 109 gaps. Now he has divided \$110 billion into 10 portions each with at least \$1 billion. Now he retrieves \$1 billion from each portion, and returns the \$10 to his friend. Now the original \$100 billion is divided into 10 piles in such a way that there is no minimum amount! So $109$ choose $9$ would give me the combinations. In R, `combn(a,b)` would give me all the cuts points from all possible combinations "a choose b". However running `combn(109,9)` is not possible on my computer. In R, is there an efficient way to sample from the output of `combn`?


title: Endogeneity in spatially lagged regression model

text: The standard convention in Spatial Statistic is that the spatial lag term in a regression model will be biased due to simultaneity. Looking at the following model, it would be difficult to argue with this: $y_{i} = \rho W y_{i} + X_{i} \beta + \varepsilon_{i}$, as we would expect a feedback effect. In other words, the spatial lag effects $y_{i}$, but $y_{i}$ also influences its own spatial lag. However, I am trying to simulate a spatial regression model and I am not finding this result. I would appreciate it if anybody was able to shed light on why this is the case. This is even the case with a very low number of observations and high spatial correlation. My simulation procedure is as follows. There are 100 observations evenly spaced on a $10 \times 10$ grid. The $W$ weights are given by the rook contiguity, and row standardized. I estimate 4 models. The first two are OLS, with and without a spatial lag term. Unsurprisingly, the coefficients are biased when the spatial lag is excluded. The other models estimate the spatial lag via maximum likelihood and 2SLS. I do not find any substantial difference between the OLS with spatial lag, maximum likelihood, and 2SLS estimators. The following variables are generated: $x_1 = 1$, i.e. the intercept $u_{1i}\sim \mathcal{N}(0,1)$ $u_{2i}\sim \mathcal{N}(0,1)$ $\rho_x = 0.5$, so the independent variables are spatially correlated $x_{2i} = (1-\rho_{x} W)^{-1} u_{1i}$ $x_{3i} = (1-\rho_{x}W)^{-1} u_{2i}$ $\rho = 0.75$, strong spatial correlation $b_1 = 1$ $b_2 = 8$ $b_3 = 2$ $\varepsilon_{i}\sim \mathcal{N}(0,1)$ $M=(1- \rho W)^{-1}$ $y_{i} = M b_{1}x_{1i} + M b_{2}x_{2i} + M b_{3}x_{3i} + M \varepsilon_{i}$, the reduced form model. I also attach the relevant R code.               # some spatial models     library(spdep)          rm(list=ls())     n = 100     data = data.frame(n1=1:n)     # coords     data$lat = rep(1:sqrt(n), sqrt(n))         data$long = sort(rep(1:sqrt(n), sqrt(n)))     # create W matrix     wt1 = as.matrix(dist(cbind(data$long, data$lat), method = "euclidean", upper=TRUE))     wt1 = ifelse(wt1==1, 1, 0)     diag(wt1) = 0     # row standardize     rs = rowSums(wt1)     wt1 = apply(wt1, 2, function(x) x/rs)     lw1 = mat2listw(wt1, style="W")          rx = 0.5     rho = 0.75     b1 = 1     b2 = 8     b3 = 2          inv1 = invIrW(lw1, rho=rx, method="solve", feasible=NULL)     inv2 = invIrW(lw1, rho=rho, method="solve", feasible=NULL)          sims = 100     beta1results = matrix(NA, ncol=4, nrow=sims)     beta2results = matrix(NA, ncol=4, nrow=sims)     beta3results = matrix(NA, ncol=4, nrow=sims)          rhoresults = matrix(NA, ncol=3, nrow=sims)          for(i in 1:sims){       u1 = rnorm(n)       x2 = inv1 %*% u1       u2 = rnorm(n)       x3 = inv1 %*% u2        e1 = rnorm(n)       y1 = b1 + b2*x2 + b3*x3 + e1       y1 = inv2 %*% y1        yl = wt1 %*% y1       data1 = data.frame(y1, x2, x3)            m1 = coef(lm(y1 ~ x2 + x3))       m2 = coef(lm(y1 ~ yl + x2 + x3))       m3 = coef(lagsarlm(y1 ~ x2 + x3, data1, lw1))       m4 = coef(stsls(y1 ~ x2 + x3, data1, lw1))            beta1results[i,] = c(m1[1], m2[1], m3[2], m4[2])       beta2results[i,] = c(m1[2], m2[3], m3[3], m4[3])       beta3results[i,] = c(m1[3], m2[4], m3[4], m4[4])       rhoresults[i,] = c(m2[2], m3[1], m4[1])     }          # assess performance     apply(rhoresults, 2, mean) ; apply(rhoresults, 2, sd)     apply(beta1results, 2, mean) ; apply(beta1results, 2, sd)     apply(beta2results, 2, mean) ; apply(beta2results, 2, sd)     apply(beta3results, 2, mean) ; apply(beta3results, 2, sd)


title: How to prove independence in the following case?

text: Let $f$ be a (possibly randomized) function with finite domain $D$ and finite range $R$, such that for any two $x$ and $x'$ in the domain, $f(x)$ and $f(x')$ are identically distributed; i.e.: $$\forall x,x'\in D, \forall y\in R \quad \Pr[f(x)=y] = \Pr[f(x')=y] \enspace. \qquad (\dagger)$$ I want to prove that the output of $f$ is independent of any input distribution. That is, for any discrete input distribution $X$ with support $D$, we have: $$\forall x\in D, \forall y\in R \quad \Pr[f(x)=y \mid X=x]=\Pr[f(x)=y] \enspace. \qquad (\ddagger)$$ * * * Here's two examples that satisfy $(\dagger)$, and therefore satisfy $(\ddagger)$: **Example 1.** $f(x)$ is a constant function. **Example 2.** Domain $D$ is the set of $n$-bit strings, and $f(x)=x \oplus r$. That is, $f$ picks a random $n$-bit string (internally), and XORs it to its input.


title: Can saved regression factor scores be used to test for normality in spss and used for further analysis?

text: I have already collected my data for my dissertation and as I have many variables I have done a factor analysis to drill down the variables to few factors. Then I have computed the mean of all the relevant variables of each individual factor. However, when I am doing a normality test on the factors, it is being negative, while normality test with the saved regression factor scores are yielding positive results. Which is correct? Also, I wanted to know if I can used the saved regression factor scores to carry out further tests like discriminant analysis, t-test and multiple regression analysis or should I stick to the mean of the factors? I am quite confused on how to proceed. Any help would be much welcomed. Thanks.


title: Estimation and model selectionfor Gaussian AR(p): conditional and unconditional log likelihood

text: Suppose a Gaussian stationary AR(p) process is $$ X_t = \sum_{i=1}^t \phi_i X_{t-i} + a_t $$ where $a_t$ is from iid $N(0, \sigma_a^2)$.   1. For estimating its parameters from a sample path of length $n$, I feel it more convenient to maximize the conditional log likelihood $\log p(X_{p+1}, \dots, X_n \mid X_1, \dots, X_p, \phi_1, \dots, \phi_p, \sigma_a^2)$ instead of maximizing the log likelihood $\log p(X_1, \dots, X_n \mid \phi_1, \dots, \phi_p, \sigma_a^2)$ , since I don't know how to derive the log likelihood for the joint distribution of $X_1, \dots, X_p$. I know $X_1, \dots, X_p$ must have a multivariate normal distribution, but I couldn't find explicit formula for their covariance matrix. So I wonder if you know an explicit formula for the covariance matrix of $X_1, \dots, X_p$ and for the log likelihood?   2. When calculating AIC, since I don't know how to calculate log likelihood, is it acceptable to use conditional log likelihood instead of log likelihood? My concern is that the difference between log likelihood and conditional log likelihood is $\log p(X_1, \dots, X_p \mid \phi_1, \dots, \phi_p, \sigma_a^2)$, which is related to the AR order p. If I use the "conditional" AIC to choose p, I worry I don't compare the conditional AIC's of AR models with different p on the same ground, because, the bigger p is, the more l ignore from the conditional AIC for not including $\log p(X_1, \dots, X_p \mid \phi_1, \dots, \phi_p, \sigma_a^2)$. So is such conditional AIC acceptable? Thanks!


title: Liblinear logisitic regression with L2 regularization for classification

text: I am trying to use the liblinear library for logistic regression with L2 regularization. However, I am finding some issues with it. For eg when choosing the cost parameter, I chose the C parameter to be from 0.5 to 10 with 10 fold cross validation. However, when I tried with higher values of C upto even 32000, the cross validation accuracy is more with C=32000. This is really strange I normally see people using up to 10 but in my case even 32000 seems to be the best. Isn't this weird?


title: Continuity correction error when using normal distribution to estimate Poisson distribution

text: Good morning Everyone, I am doing a self-study exercise which attempts to exemplify a case where the **Normal Distribution** is used to approximate the **Poisson Distribution** , since the population mean is more than 10. I understand that when using **Normal Distribution** to Approximate **Poisson / Binomial Distributions** , there is a need for **Continuity Correction Error** to be managed. From my understanding, this means subtracting 0.5 from the lower bound and adding 0.5 to the upper bound. ![enter image description here](http://i.stack.imgur.com/KJJud.jpg) However, I was faced with the following example, where 0.5 was subtracted from the equation. I will like to clarify: > if this is because it is lower bound, as the value is less than W? and if the following concept is correct > From my understanding, this means subtracting 0.5 from the lower bound and > adding 0.5 to the upper bound. Appreciate any advice please.


title: Parallel lines assumption for ordinal regression model

text: I used the lrm function of the RMS package for a ordinal regression model for prediction. When I use the graphical plot method, the parallel line assumption seems to fail for some variables. Is there any way the parallel lines assumption can be relaxed ? ? Does the lrm package get around the assumption ? Another question would be is whether Hosmer–Lemeshow test would be valuable in interpreting this despite the invalid assumption ? Is there an implementation in R for HS gof test for ordinal regression ? Your response will be much appreciated.


title: Should 'do not know' be included when reporting summaries of survey results?

text: When reporting the results of surveys, we often list the percentages of respondents who gave different answers. However, these percentages could be calculated including people who said 'do not know' and/or 'no basis to judge' (in which case the percentages for the substantive responses would not add up to 1), or without such people. Why would I use percentages excluding the dnk/nbj data?


title: Linear model predictor selection. Which method to use ?

text: From what I understand, there are 3 main types of predictor selection method for linear models, namely, 1 Subset Selection, 2 Shrinkage and 3 Dimension Reduction.   1. The subset selection includes the Best Subset Selection and the Stepwise Selection which could be forward, backward or hybrid. AIC, BIC, Cp or Adjusted R-Square can be used to select the predictors.   2. The Shrinkage includes the Ridge Regression and Lasso. This approach attempts to shrink the coefficients to 0   3. Dimension reduction transforms the predictors and fit the model using the transformed predictor. If forecasting accuracy is my main goal and model interpretability is not important, which method(s) should I use ? What should I do if the methods give inconsistent results ? What are the main advantages and disadvantages of each approach ?


title: How to perform a MANOVA with a random effect in R?

text: I performed an experiment where I raised different families coming from two different source populations, where each family was split up into a different treatments. After the experiment I measured several traits on each individual. Now I would like to have an overall and hence multivariate statistic that tests for an effect of either treatment or source as well as their interaction but correcting for family effects. So, basically I would like to perform a MANOVA with a random effect, which would translate to something like this:                manova(data~source*treatment+Error(family))      However, the standard R MANOVA function does not support random effects and I was unable to find another function that would do the trick. So I would like to ask if anyone has a suggestion that may help. Thanks!


title: In definition of gaussian copula does the marginals also have to be gaussian?

text: I am quite new to this copula idea. In particular I am confused about the definition of a Gaussian copula. For a copula to be a Gaussian copula does the marginals have to Gaussian as well? Or it can be of any distribution? From the wikipedia page it looks like it has to be (http://en.wikipedia.org/wiki/Copula_(probability_theory)#Gaussian_copula) but I thought it needn't be.


title: Why do the normal and log-normal density functions differ by a factor?

text: If a random variable $W$ is Normally distributed, then $\exp(W)$ is Log- Normally distributed. However, the pdfs of these two random variables differ by a factor of $\exp(W)^{-1}$. The Normal pdf for $W$ is $$P(w \in W) = \frac{1}{\sigma\sqrt{2\pi}}\, e^{-\frac{(w - \mu)^2}{2 \sigma^2}}$$ but the Log-Normal pdf for $\exp(W)$ is $$P(\exp(w) \in \exp(W)) = \frac{1}{\exp(w)\sigma\sqrt{2\pi}}\ e^{-\frac{\left(w-\mu\right)^2}{2\sigma^2}}$$ Why is this the case? Shouldn't they be equivalent?


title: Probability Given Datapoint Does Not Appear in a Bootstrap Sample?

text: I received this question for a Statistics / Machine Learning assignment and I'd like to you if any of you know the proper answer. If we have n data points, what is the probability that a given data point does not appear in a bootstrap sample? Sounds simple enough right? I'm reading Introduction to Statistical Learning to try and find the solution but I would definitely appreciate some help


title: Comparing general search engine and a meta search engine. Paired t test or independent t-test?

text: I have built a domain-specific search engine, a meta-search engine. This search engine, of course, take its results from other search engines like Google, Bing, Yahoo and stuff. Then does some classification to filter the results. Now, I'm doing an evaluation on the precision of search results returned by this search engine towards the precision of search results returned from a general search engine. So I queried a number of queries to my search engine, count the precision of each result returned, and compute the average of the precisions. I do the same queries to a general search engine, which in this case is Google. And by so compute the average of the precisions. I use t-test two see if the difference between these precisions statistically significantly differ. What I can't understand now, is it paired t-test or independent one? I heard from my professor that it should be independent t-test, since the sample sets come from two different search engines. But in my mind, since the sample (search results returned) are displayed based on same queries in order, and meta-search engine takes its results from several general search engines, it should be paired t-test. On the other hand I know it as well, that paired t-test should be done on SAME sample, different treatments. Which do I choose? Need help. I enclose the screen shoot of the search results as well, to make my question more clear and understandable. Thanks. ![How it looks like](http://i.stack.imgur.com/o4RsK.png)


title: Using proper scoring rule to determine class membership from logistic regression

text: I am using logistic regression to predict likelihood of an event occurring. Ultimately, these probabilities are put into a production environment, where we focus as much as possible on hitting our "Yes" predictions. It is therefore useful for us to have an idea of what definitive "hits" or "non-hits" might be _a priori_ (before running in production), in addition to other measures we use for informing this determination. My question is, what would be the proper way to predict a definitive class (1,0) based on the predicted probability? Specifically, I use R's `glmnet` package for my modeling. This package arbitrarily picks .5 probability as threshold for a yes or no. I believe that I need to take the results of a proper scoring rule, based on predicted probabilities, to extrapolate to a definitive class. An example of my modeling process is below:               mods <- c('glmnet', 'scoring')     lapply(mods, require, character.only = T)          # run cross-validated LASSO regression     fit <- cv.glmnet(x = df1[, c(2:100)]), y = df1[, 1], family = 'binomial',      type.measure = 'auc')          # generate predicted probabilities across new data     df2$prob <- predict(fit, type="response", newx = df2[, c(2:100)], s = 'lambda.min')          # calculate Brier score for each record     df2$propscore <- brierscore(df2[,1] ~ df2$prob, data = df2)      So I now have a series of Brier scores for each prediction, but then how do I use the Brier score to appropriately weight each likelihood being a yes or no? I understand that there are other methods to make this determination as well, such as Random Forest.


title: Probability of observed data in HMM

text: In a given Gaussian mixture model with observed continues variables $Y$ and latent discrete variables $X$ I want to apply the forward-backward algorithm in order to compute the marginal posteriors $P(x_t|y_{1:T})$. Since this is computed as $$\frac{\alpha_t(x_t) \beta_t(x_t)}{P(Y)}$$ I was wondering how do I obtain the value of $P(Y)$? The only probabilities I have given is a transition probability $P(x'|x)$.


title: Why does bagging use bootstrap samples?

text: Bagging is the process of creating N learners on N different bootstrap samples, then taking the mean of their predictions. My question is: Why not use any other type of sampling? Why use bootstrap samples?