title: Frequentist performance of sequential testing based on predictive power

text: Suppose we perform a hypothesis test from a random sample $(x_i)_{i=1}^n$, assuming for instance $x_i \sim_{\text{iid}} {\cal N}(\theta, 1)$ and $H_0=\{\theta=0\}$ for simplicity. If the test fails to reject $H_0$, it is natural to ask oneself _"if I collect more data, what could happen?"_. This question can be solved by assigning a "predictive distribution" to the new data, and then calculate the probability to reject $H_0$ for a new sample assumed to be distributed according to this predictive distribution. A naive way (suffering from a lack of probabilist interpretation) consists in saying that the new data are $\sim_{\text{iid}} {\cal N}(\hat\theta, 1)$ where $\hat\theta$ is an estimate of $\theta$ obtained from the actual sample $(x_i)_{i=1}^n$. In the Bayesian framework there is a clear notion of predictive distribution. Assuming we consider either the "naive" predictive distribution or the Bayesian predictive distribution derived from a noninformative prior, how to highlight/assess the frequentist performance of the predictive power ? For instance, assume that we adopt the following sequential testing methodology:   * Collect a sample $(x_i)_{i=1}^n$ with $n=10$   * If $H_0$ is rejected then stop.   * Otherwise use the predictive distribution to evaluate the required new sample size to get $80\%$ predictive power to reject $H_0$.   * Then collect the new sample and perform the test. Adopting this methodology, what about the probability to reject $H_0$ in function of $\theta$ ? It would be easy to explore such questions using simulations, but does there exist some theoretical results about such questions ?


title: Tolerance interval for Deming regression

text: I am trying to derive (one-sided) tolerance intervals related to the Deming regression model: $$ x_i=x^*_i + \epsilon_i$$ $$ y_i = (\alpha+\beta x^*_i) + \epsilon'_i$$ where the $x^*_i$'s are nonrandom fixed numbers, $\epsilon_i \sim {\cal N}(0,\sigma_x^2)$, $\epsilon'_i \sim {\cal N}(0,\sigma_y^2)$, and all variables $\epsilon_i, \epsilon'_i$ are mutually independent. The ratio $\sigma_x^2/\sigma_y^2$ is assumed to be known. Put $\tau^2=\sigma_x^2+\sigma_y^2$. Given a new ``true $x$'' value $x^*_\text{new}$ I'm looking for an upper tolerance bound of the difference $$y_\text{new} - x_\text{new} \sim {\cal N}(\alpha+(\beta-1)x^*_\text{new}, \tau^2).$$ Actually this upper tolerance bound is nothing but an upper $100(1-\alpha)\%$-confidence bound about the $100p\%$-quantile of $y_\text{new} - x_\text{new}$. As for classical Gaussian models I'm looking for a tolerance bound having form $$B(x,y)=(\hat{y^*}_\text{new} - x^*_\text{new}) + K s$$ with $\hat{y^*}_\text{new} = \hat\alpha + \hat\beta x^*_\text{new}$ and $s^2= \hat\tau^2$. Denoting by $F_\text{new}$ the cumulative distribution function of the difference $y_\text{new} - x_\text{new}$, the tolerance factor $K$ must be derived in order that $$\Pr\left( F_\text{new}(B(x,y)) > p \right) = 1-\alpha$$ for some given $p\in]0,1[$ (the tolerance level) and confidence level $1-\alpha$ (of course here $\alpha$ is not the intercept parameter $\alpha$ of the model!) Estimates of $\alpha$, $\beta$ and $\tau^2$ as well as their asymptotic variance matrix $V$ and an estimate $\hat V$ of this asymptotic variance matrix are available in some textbooks and/or research papers (see for instance Fuller's book or Iles & Gillard's report 1 and report2). In particular $\hat\tau^2$ is asymptotically independent of $(\hat\alpha,\hat\beta)$. Now I'm trying to derive the tolerance bound by miming the way used for deriving the tolerance bound in a classical regression model. It is based on the following lemma. **Lemma** Let $W \sim {\cal N}(0, \sigma^2)$ independent of $Q \sim \frac{\chi^2_{d}}{d}$. Then the number $K>0$ satisfying $$\Pr\left( \Phi(W+K\sqrt{Q}) > p \right) =1-\alpha$$ is given by a known formula $K=K(\sigma,d,p,\alpha)$ (see Krishnamoorthy & Mathew's book). Now put $X_\text{new}=(1,x_\text{new})'$. One has the three asymptotically independent random variables $$ \frac{(y_\text{new} - x_\text{new} ) - (\alpha+(\beta-1)x^*_\text{new})}{\sqrt{\tau^2}} \sim {\cal N}(0,1),$$ $$ \frac{(\hat{y^*}_\text{new} - x^*_\text{new} ) - (\alpha+(\beta-1)x^*_\text{new})}{\sqrt{\tau^2}} \sim {\cal N}\left(0,\frac{X_\text{new}'VX_\text{new}}{\tau^2}\right),$$ $$\frac{s^2}{\tau^2} \sim \frac{\chi^2_{n-2}}{n-2}$$ The definition of $K$ becomes $$\Pr\left( \Phi(W+K\sqrt{Q}) > p \right) =1-\alpha$$ with $W \sim {\cal N}\left(0,\frac{X_\text{new}'VX_\text{new}}{\tau^2}\right)$ independent of $Q \sim \frac{\chi^2_{n-2}}{n-2}$. Applying the lemma we get $K$ depending on the unknown values of $\tau$ and $V$. Substituting with the estimates $\hat\tau$ and $\hat V$ gives an approximate tolerance factor $\hat K$. Well, I have done it and I have implemented everything in R. But assessing the performance of my tolerance bound with simulation shows that even with large sample sizes _it does not work_ , except for the case $\beta=1$. I'm pretty sure of my formulas for the estimates and the variance matrices since I have also implemented confidence intervals and prediction intervals and they work well. **EDIT:** In fact that works !! There was a silly error in the R code I used for simulations. Should I delete this post ? I think no because the method could interest others.


title: LogLikelihood Parameter Estimation for Linear Gaussian Kalman Filter

text: I have written some code that can do Kalman filtering (using a number of different Kalman-type filters [Information Filter et al.]) for Linear Gaussian State Space Analysis for an n-dimensional state vector. The filters work great and I am getting some nice output. However, the parameter estimation via loglikelihood estimation is confusing me. I am not a statistician but a physicist, so please be kind. Let us consider the linear Gaussian State Space model $$y_t = \mathbf{Z}_{t}\alpha_{t} + \epsilon_{t},$$ $$\alpha_{t + 1} = \mathbf{T}_{t}\alpha_{t} + \mathbf{R}_{t}\eta_{t},$$ where $y_{t}$ is our observation vector, $\alpha_{t}$ our state vector at time step $t$. The quantities in bold are the transformation matrices of the state space model which are set according to the characteristics of the system under consideration. We also have $$\epsilon_{t} \sim NID(0, \mathbf{H}_{t}),$$ $$\eta_{t} \sim NID(0, \mathbf{Q}_{t}),$$ $$\alpha_{1} \sim NID(a_{1}, \mathbf{P}_{1}).$$ where $t = 1,\ldots, n$. Now, I have derived and implemented the recursion for the Kalman Filter for this generic state space model by guessing the initial parameters and variance matrices $\mathbf{H}_{1}$ and $\mathbf{Q}_{1}$ I can produce plots like ![Kalman Filter](http://i.stack.imgur.com/QXS9F.png) where the points are the Nile River water levels for Jan over 100 years, the line is the Kalamn Estimated state, and the dashed lines are the 90% confidence levels. Now, for this 1D data set the matrices $\mathbf{H}_{t}$ and $\mathbf{Q}_{t}$ are just scalars $\sigma_{\epsilon}$ and $\sigma_{\eta}$ respectively. So now I want to get the correct parameters for these scalars using the output from the Kalman Filter and the loglikelihood function $$\log L(Y_{n}) = -\frac{np}{2}\log(2\pi) - \frac{1}{2}\sum^{n}_{t = 1}(log|\mathbf{F}_{t}| + v^{\mathsf{T}}_{t}\mathbf{F}_{t}^{-1}v_{t})$$ Where $v_{t}$ is the state error and $\mathbf{F}_{t}$ is the state error variance. Now, here's where I am confused. From the Kalman filter, I have all the information I need to work out $L$, but this seems to get me no closer to being able to calculate the maximum likelihood of $\sigma_{\epsilon}$ and $\sigma_{\eta}$. **My question is how can I calculate the maximum likelihood of $\sigma_{\epsilon}$ and $\sigma_{\eta}$ using the loglikelihood approach and the equation above?** An algorithmic break down would be like a cold beer to me right now... Thanks for your time. * * * _Note. For the 1D case $\mathbf{H}_{t} = \sigma^{2}_{\epsilon}$ and $\mathbf{H}_{t} = \sigma^{2}_{\eta}$. This is the univariate local level model._


title: Learning what variables have high impact on one other variable

text: I'm using R. My dataset has about 39 different Variables/Vektors and each has about 80 entries.Furthermore I have 1 Variable with 80 entries I want to maximize (lets call it Y). I assume that there are some correlations in my data and I want to find out what values my 39 variables should have in order for my Y to be maximal. I'm also very sure that not all variables have significant impact on my Y. What i already did is calculate the correlations of the 39 Variables to my Y and by knowing which correlations where very high/very low I can assume what values they should have to maximize Y. But im not 100% confident that this is the best way to approach my problem because you can't draw big conclusions just from correlation. However I do not know what else I should do and where I should look to solve this problem.


title: Elasticity of logistic regression model

text: Let the standard binary logit model with utility function: Uni = Vni + eni, where eni obeys iid gumbel, and define Vni as: Vni = a + xni*b. where n denotes individual, and i denotes alternative. Then, the elasticity of xni is given by: Ei = b * xin * (1-Pni) My question is if Utility function is not linear in x but log-linear as: Vni = a + log(xni) * b, then how can we calculate the elasticity of x ? Any suggestions are appreciated. Thanks in advance.


title: Find 1 dimensional sufficient statistic for $Beta(\alpha, 2\alpha)$

text: I am currently getting $\prod x_i$ and $\prod (1-x_i)$ but a 1 dimensional sufficient statistic is asked for and considering that there is only a single parameter ($\alpha$), I think we should get a 1-dimensional sufficient statistic. I am not sure how to go about doing so though.


title: When to use Central Limit Theorem?

text: 1. Is the central-limit-theorem true for all distributions?   2. For what sample size is it true (what is meant with "large")?


title: Long-tailed distribution of time events

text: Suppose you have the logs of a web server. In these logs you have tuples of this kind:               user1, timestamp1     user1, timestamp2     user1, timestamp3     user2, timestamp4     user1, timestamp5     ...      These timestamps represent e.g. users' clicks. Now, `user1` will visit the site multiple times (sessions) during the month, and you'll have bursts of clicks from each user during each session (supposing that when a user visits your site, he'll click on multiple pages). Suppose you want to partition these burst of clicks in the sessions that generated them, but you don't have any additional source of information, only the list of timestamps. If you compute the distribution of intervals between two consequent clicks from the same user, you will obtain a long-tailed distribution. Intuitively, you'd look for a "cut parameter", e.g. N seconds, where if `timestamp_{i+1} - timestamp{i} > N`, then your `timestamp_{i+1}` is the beginning of the new session. The problem is that this distribution in reality is a mixture of two variables: X = "interval between two consequent clicks in the same session" and Y = "interval between the last click of the previous session and the first of the new one". The question is, how to estimate this N, that divides the two distributions (with a bit of overlap, possibly) just by looking at the burst of clicks?


title: How to get Sub-Training and Sub-Test from cross validation in Caret

text: I am using Caret function I have divided my data into training(75%) and test (25%) sets. Now I am running 10-Fold CV on training data. When i fit following model train_control <- trainControl(method="cv", number=10) output <- train(Species~., data=iris, trControl=train_control, method="rpart2") Q1. Is the accuracy shown by output an average Sub-Test set or Sub-Training i.e. (10% of Training data or 90% of Training data) Q2. How can I view the accuracy of both sets for each folds i.e. accuracy of Sub-Test and Sub- Training set for each of 10-Folds


title: Confidence intervals of coefficients in polynomial regression

text: Given a model $a_{n}\cdot x^{n}+ a_{n-1}\cdot x^{n-1}+...+ a_0 = y + e$, where $e$ is some unknown (but assumed normally distributed) error and $x$ and $y$ are observations. Now, how would one to estimate the confidence intervals of the coefficients $a_.$? AFAIU, a simple procedure based on the covariance $(X'X)^{-1}$ is not sufficient, where $X$ is the $n$'th order (zero mean) Vandermonde matrix of $x$. First of all $X'X$ is singular. Also calculation $var(a)= (r'r/m)\cdot diag(pinv(X'X))$, where $r$ is the residual and $m$ number of observations, seems to be vague, because high correlation between columns of $X'X$.


title: How to take cumulative sum of a 2D distribution?

text: How can I take the cumulative sum of a 2D distribution, starting from the highest value in the distribution? Background: my 2D distribution is the probability of different combinations of wave heights and wave periods based on measurements made at an ocean site. Assume the distribution shows that the most common wave (highest probability, 1%) has a height of 3 meters and a period of 8 seconds. This combination represents the "peak" of a "mountain" in the distribution, i.e. adjacent combinations of wave height/period occur often but less frequently (say 0.1%), and combinations further "downhill" occur even less frequently (say 0.01%). Non-zero probability is spread out over many bins of wave height and wave period so it still totals to 1. I want to be able to identify the sea states that occur 10% of the time, i.e. make a contour that encloses the height/period combinations that together add up to a probability of 10%. In other words, I can look at the contour and say "10% of the time, wave heights are between 1 and 4 meters and periods are between 6 and 10 seconds". I believe I have to contour the **cumulative sum of the 2D distribution starting at the "peak" and working "downhill"**. What is confusing is that downhill extends in 2 directions since I am summing over 2 variables, height and period. Any advice is helpful since I am having trouble wrapping my mind around this.


title: Mediation using 'mediate' in R with an interaction

text: I have a question about mediation when there is an interaction between 2 categorical variables, one with two levels and one with three levels. Using the ‘mediation’ package by Tingley and colleagues in R, I’m trying to conduct a bootstrapped mediation analysis for an interaction between a 2-level categorical IV (‘condstorm’) and another 3-level categorical IV (‘condmean’). The mediator (‘trustmed2’) is dichotomous, and the dependent variable (‘AttInd’) is continuous. Categorical 2-level IV ‘condstorm’               table(d$condstorm)              nostorm   storm          604     570       Categorical 3-level IV ‘condmean’               table(d$condmean)          both high mean      413  375  386       Dichotomous mediator ‘trustmed2’               table(d$trustmed2)              0   1          305 866      Continuous DV ‘AttInd’               summary(d$AttInd)     Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's      0.0000  0.3571  0.6200  0.5801  0.8214  1.0000       3       The code and output for the model with the mediator is:               contrasts(d$condstorm)         storm     nostorm     0     storm       1          contrasts(d$condmean)          nonevhigh nonevboth     both         0         1     high         1         0     mean         0         0               med.fit <- glm(trustmed2 ~ condmean*condstorm + gwbelief + trustmed1 + coastdwell + partyID + polor + gender + age + race + educ + income + region, weights=weight1, data = d, family = binomial('probit'))          summary(med.fit)          Call:     glm(formula = trustmed2 ~ condmean * condstorm + gwbelief + trustmed1 +  coastdwell + partyID + polor + gender + age + race + educ + income + region, family = binomial("probit"), data = d, weights = weight1)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.9982  -0.1429   0.2196   0.4594   3.2721            Coefficients:                          Estimate Std. Error z value Pr(>|z|)         (Intercept)          -1.876115   0.284494  -6.595 4.27e-11 ***     condmeannonevhigh     0.170809   0.175454   0.974  0.33029         condmeannonevboth     0.497041   0.181383   2.740  0.00614 **      condstormstorm        0.285569   0.182081   1.568  0.11680         gwbeliefnotvref       0.056512   0.631022   0.090  0.92864         gwbeliefnotvprobs     0.746556   0.129610   5.760 8.41e-09 ***     trustmed1             1.908566   0.121064  15.765  < 2e-16 ***     coastdwellyes         0.320369   0.128445   2.494  0.01262 *       partyIDvdem           0.585477   0.148862   3.933 8.39e-05 ***     partyIDvref          -0.468275   0.654078  -0.716  0.47403         partyIDvrep           0.058146   0.130621   0.445  0.65621         polorvref            -0.316949   0.508007  -0.624  0.53269         polorvlib            -0.058482   0.160609  -0.364  0.71577         polorvcons           -0.369565   0.122246  -3.023  0.00250 **      genderMale           -0.009229   0.105253  -0.088  0.93013         agev1824              0.235127   0.194856   1.207  0.22756         agev2534              0.099723   0.178256   0.559  0.57587         agev3544              0.119399   0.169987   0.702  0.48243         agev4554              0.221921   0.176420   1.258  0.20842         agev5564              0.173343   0.170355   1.018  0.30890        racevhisp             0.033141   0.177162   0.187  0.85161         racevblack           -0.323867   0.188305  -1.720  0.08545 .       racevother            0.475179   0.234313   2.028  0.04256 *       educhsvsome          -0.002255   0.132574  -0.017  0.98643         educhsvcollgrad       0.201980   0.146240   1.381  0.16723         educhsvref            1.304337   0.666427   1.957  0.05032 .       incomev3050          -0.095963   0.159290  -0.602  0.54688         incomev5075           0.184528   0.168946   1.092  0.27473         incomev75100          0.028269   0.191474   0.148  0.88263         incomev100           -0.249873   0.167648  -1.490  0.13610         regionvmw            -0.058581   0.163619  -0.358  0.72032         regionvne             0.315495   0.177227   1.780  0.07505 .       regionvs              0.204790   0.145384   1.409  0.15895         condmeannonevhigh:condstormstorm -0.061962   0.263472  -0.235  0.81407         condmeannonevboth:condstormstorm -0.665467   0.253958  -2.620  0.00878 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 1316.98  on 1170  degrees of freedom     Residual deviance:  767.03  on 1136  degrees of freedom       (3 observations deleted due to missingness)     AIC: 836.49          Number of Fisher Scoring iterations: 6      And the code and output for the model with the mediator included as a predictor on the dependent variable:               out.fit <- lm(AttInd ~ trustmed2 + condmean*condstorm + trustmed1 + gwbelief + coastdwell + partyID + polor + gender + age + race + educ + income + region, weights=weight1, data = d)          summary(out.fit)          Call:     lm(formula = AttInd ~ trustmed2 + condmean * condstorm + trustmed1 +          gwbelief + coastdwell + partyID + polor + gender + age +          race + educ + income + region, data = d, weights = weight1)          Weighted Residuals:          Min       1Q   Median       3Q      Max      -0.79365 -0.12054  0.00739  0.12351  0.72377           Coefficients:                       Estimate Std. Error t value Pr(>|t|)         (Intercept)         1.754e-01  3.411e-02   5.143 3.18e-07 ***     trustmed2           2.777e-01  1.940e-02  14.313  < 2e-16 ***     condmeannonevhigh  -8.754e-03  2.174e-02  -0.403 0.687311         condmeannonevboth   2.548e-02  2.148e-02   1.186 0.235756         condstormstorm      3.977e-02  2.250e-02   1.767 0.077433 .       trustmed1           2.084e-02  2.014e-02   1.035 0.301053         gwbeliefnotvref    -1.063e-01  9.071e-02  -1.172 0.241589         gwbeliefnotvprobs   1.574e-01  1.818e-02   8.656  < 2e-16 ***     coastdwellyes       1.100e-02  1.485e-02   0.740 0.459293         partyIDvdem         6.283e-02  1.618e-02   3.884 0.000109 ***     partyIDvref         -2.558e-01  8.716e-02  -2.934 0.003408 **      partyIDvrep         -4.893e-02  1.688e-02  -2.898 0.003829 **      polorvref            8.509e-02  6.617e-02   1.286 0.198751         polorvlib            6.199e-02  1.742e-02   3.558 0.000389 ***     polorvcons          -7.015e-03  1.571e-02  -0.447 0.655204         genderMale          -5.778e-03  1.267e-02  -0.456 0.648363         agev1824             4.624e-03  2.391e-02   0.193 0.846692         agev2534            -1.773e-02  2.222e-02  -0.798 0.425121         agev3544             2.952e-02  2.107e-02   1.401 0.161448         agev4554             1.590e-02  2.159e-02   0.736 0.461613         agev5564            -9.967e-03  2.073e-02  -0.481 0.630765         racevhisp            4.936e-02  2.030e-02   2.432 0.015174 *       racevblack          -1.585e-02  2.224e-02  -0.713 0.476104         racevother           4.969e-02  2.518e-02   1.973 0.048687 *       educhsvsome         -1.016e-02  1.612e-02  -0.630 0.528528         educhsvcollgrad     -1.209e-02  1.713e-02  -0.706 0.480595         educhsvref           2.383e-01  8.213e-02   2.901 0.003791 **      incomev3050         -9.128e-04  1.964e-02  -0.046 0.962941         incomev5075          1.522e-02  1.996e-02   0.762 0.446132         incomev75100        -3.185e-02  2.280e-02  -1.397 0.162813         incomev100           1.344e-02  2.013e-02   0.668 0.504422         regionvmw            3.467e-02  2.027e-02   1.710 0.087541 .       regionvne            1.861e-02  2.015e-02   0.924 0.355812         regionvs             -7.531e-05  1.716e-02  -0.004 0.996499         condmeannonevhigh:condstormstorm -4.354e-03  3.152e-02  -0.138 0.890175         condmeannonevboth:condstormstorm -5.284e-02  3.090e-02  -1.710 0.087514 .       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.2103 on 1135 degrees of freedom       (3 observations deleted due to missingness)     Multiple R-squared:  0.4546,    Adjusted R-squared:  0.4378      F-statistic: 27.03 on 35 and 1135 DF,  p-value: < 2.2e-16      The model without the mediator is:               out.fit <- lm(AttInd ~ condmean*condstorm + trustmed1 + gwbelief + coastdwell +                      partyID + polor + gender + age + race + educ + income +                       region, weights=weight1, data = d)          summary(out.fit)          Call:     lm(formula = AttInd ~ condmean * condstorm + trustmed1 + gwbelief +          coastdwell + partyID + polor + gender + age + race + educ +          income + region, data = d, weights = weight1)          Weighted Residuals:         Min      1Q  Median      3Q     Max      -1.1727 -0.1371  0.0130  0.1366  0.8077           Coefficients:                                Estimate Std. Error t value Pr(>|t|)         (Intercept)               0.1940313  0.0370161   5.242 1.89e-07 ***     condmeannonevhigh        -0.0019461  0.0236085  -0.082  0.93432         condmeannonevboth         0.0474046  0.0232707   2.037  0.04187 *       condstormstorm            0.0527524  0.0244159   2.161  0.03094 *       trustmed1                 0.1833521  0.0180674  10.148  < 2e-16 ***     gwbeliefnotvref          -0.1093369  0.0985135  -1.110  0.26729         gwbeliefnotvprobs         0.2074035  0.0193781  10.703  < 2e-16 ***     coastdwellyes             0.0241253  0.0161011   1.498  0.13432         partyIDvdem               0.0858665  0.0174835   4.911 1.04e-06 ***     partyIDvref              -0.2999944  0.0945990  -3.171  0.00156 **      partyIDvrep              -0.0468268  0.0183366  -2.554  0.01079 *       polorvref                 0.0560635  0.0718304   0.780  0.43526         polorvlib                 0.0568556  0.0189179   3.005  0.00271 **      polorvcons               -0.0293808  0.0169730  -1.731  0.08372 .       genderMale               -0.0072970  0.0137554  -0.530  0.59588         agev1824                  0.0162995  0.0259507   0.628  0.53007         agev2534                 -0.0134314  0.0241301  -0.557  0.57790         agev3544                  0.0326065  0.0228762   1.425  0.15433         agev4554                  0.0281826  0.0234333   1.203  0.22935         agev5564                 -0.0042263  0.0225102  -0.188  0.85111        racevhisp                 0.0484325  0.0220418   2.197  0.02820 *       racevblack               -0.0292002  0.0241349  -1.210  0.22658         racevother                0.0682464  0.0273067   2.499  0.01259 *       educhsvsome              -0.0130208  0.0175088  -0.744  0.45723         educhsvcollgrad          -0.0052682  0.0185980  -0.283  0.77703         educhsvref                0.3522019  0.0887800   3.967 7.73e-05 ***     incomev3050              -0.0065261  0.0213276  -0.306  0.75967         incomev5075               0.0200734  0.0216792   0.926  0.35468         incomev75100             -0.0336651  0.0247648  -1.359  0.17429         incomev100               -0.0006802  0.0218303  -0.031  0.97515         regionvmw                 0.0293874  0.0220138   1.335  0.18216         regionvne                 0.0317515  0.0218553   1.453  0.14655         regionvs                  0.0085483  0.0186236   0.459  0.64632         condmeannonevhigh:condstormstorm -0.0068725  0.0342330  -0.201  0.84092         condmeannonevboth:condstormstorm -0.0860072  0.0334623  -2.570  0.01029 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.2284 on 1136 degrees of freedom       (3 observations deleted due to missingness)     Multiple R-squared:  0.3562,    Adjusted R-squared:  0.3369      F-statistic: 18.49 on 34 and 1136 DF,  p-value: < 2.2e-16      We want to examine mediation in this interaction effect. We have tried using the code: mediate(med.fit, out.fit, treat = 'condstorm', control.value = 'nostorm', treat.value = 'storm', mediator = 'trustmed2', covariates=list(condmean='both'), boot=TRUE, sims=5000) Which gave the following output:               Causal Mediation Analysis           Confidence Intervals Based on Nonparametric Bootstrap          (Inference Conditional on the Covariate Values Specified in `covariates')          Estimate 95% CI Lower 95% CI Upper p-value          ACME           -0.017705    -0.036626    -0.003284    0.04     ADE            -0.014082    -0.045811     0.023155    0.53     Total Effect   -0.031788    -0.069698     0.012371    0.14     Prop. Mediated  0.516591     1.174065     8.857665    0.04          Sample Size Used: 1171           Simulations: 5000      But this mediation only compares respondents at the level of “both” for the condmean variable and excludes other conditions, so examines the effect of "Both Condmean x Condstorm" MINUS the effect "condstorm", which is -.08 - (.05) = -.03 (the total effect in the output above). The total effect should instead be the same as the coefficient of "condmeannonevboth:condstormstorm" in the second interaction term in the model, which is -.086, p < .01. Is it possible for us to examine this effect using the ‘mediation’ package in R? Any help would be much appreciated!


title: How do I randomize a larger population, from an existing population in R?

text: I have a small 20k (176 variables) test dataset. I want to create a random bigger dataset for testing, using the rows (data) and columns (names) from the small dataset. Currently I have:               real_data <- read.csv("somefilehere.csv")          # load one fake row so assignment of column names won't *ahem* complain     test_data <- data.frame(real_dat[1,])     colnames(test_data) <- colnames(real_data)      So, this sets my test_data dataframe to have the same column names as the real data. How can I run through all the column names and then randomize test_data with existing values from real_data? The following code works perfectly for randomizing a sample from the existing row:               sample(unique(real_data$some_variable_here),random_count,replace=T)      **How do I use iteration of the columns, together with the sample function to populate test_data automatically?** Doing manual assignments for every variable would take forever, maybe there is a function in R that I'm missing as a 3 week beginner? I'm working with R, knitr and LaTeX for processing and would like to use the test_data set for use in my statistical tests and graphs (to see if my document scales well with larger numbers, and my legends still look ok etc). [edit] An example with randomly picked R dataset. Dataset: warpbreaks - load with data(warpbreaks) ![dataset](http://i.stack.imgur.com/4qAck.png) Currently it has 3 already defined columns and 54 observations. How do I create a new dataset called warpbreaks_huge, that has the same column names, but say 25 million observations? Also, from the existing dataset we can see that the tension column has default values L, M and H. The wool column has default values of A and B. I guess my question is, how do you create a randomized _bigger_ populations programmatically? Apologies for my naive phrasing. I'm new at this environment.


Can a confounding variable be correlated with the DV and not the IV? I have heard of the DV being corr. but I can't find IV in any textbooks. I found this def. in Wiki: http://en.wikipedia.org/wiki/Confounding


title: How to assess the fit of a binomial GLMM fitted with lme4 (> 1.0)?

text: I have a GLMM with a binomial distribution and a logit link function and I have the feeling that an important aspect of the data is not well represented in the model. To test this, I would like to know whether or not the data is well described by a linear function on the logit scale. Hence, I would like to know whether the residuals are well-behaved. However, I can not find out at which residuals plot to plot and how to interpret the plot. Note that I am using the new version of lme4 (the development version from GitHub):               packageVersion("lme4")     ## [1] ‘1.1.0’      My question is: **How do I inspect and interpret the residuals of a binomial generalized linear mixed models with a logit link function?** The following data represents only 17% of my real data, but fitting already takes around 30 seconds on my machine, so I leave it like this:               require(lme4)     options(contrasts=c('contr.sum', 'contr.poly'))          dat <- read.table("http://pastebin.com/raw.php?i=vRy66Bif")     dat$V1 <- factor(dat$V1)          m1 <- glmer(true ~ distance*(consequent+direction+dist)^2 + (direction+dist|V1), dat, family = binomial)      The simplest plot (`?plot.merMod`) produces the following:               plot(m1)      ![enter image description here](http://i.stack.imgur.com/6eQ4g.png) Does this already tell me something?


title: How can I model a proportion with BUGS/JAGS/STAN?

text: I am trying to build a model where the response is a proportion (it is actually the share of votes a party gets in constituencies). Its distribution is not normal, so I decided to model it with a beta distribution. I also have several predictors. However, I don't know how to write it in BUGS/JAGS/STAN (JAGS would be my best choice, but it doesn't really matter). My problem is that I make a sum of parameters by predictors, but then what can I do with it? The code would be something like this (in JAGS syntax), but I don' know how to "link" the `y_hat` and `y` parameters.               for (i in 1:n) {      y[i] ~ dbeta(alpha, beta)           y_hat[i] <- a + b * x[i]     }      (`y_hat` is just the cross-product of parameters and predictors, hence the deterministic relationship. `a` and `b` are the coefficients which I try to estimate, `x` being a predictor). Thanks for your suggestions!


title: Formal way of defining a weighted average

text: I am looking for a nice and formal way to define moving averages : $y = \sum_{i=0}^N w_ix_i=<\mathbf{w},\mathbf{x}>$ and especially the fact that weights have to equal one, i.e. $\sum_{i=0}^N w_i=1$. This property seems quite obvious to me, as you somehow want the output to have the same "level" as the input values. However, I don't know how to describe it a more mathematical way .. Can somebody help me out ? Thank you


title: Robust variance estimators

text: Given $N$ data points in $\mathbb{R}^p$ - some of which are outliers (drawn from a different distribution from the inliers) - what sorts of algorithms have been designed to estimate the robust variance of this set? That is we would like to approximate the variance of the inlier space under some assumptions on the outliers perhaps.


title: Is my dataset suitable for a mixed effects model?

text: I've been putting a lot of work over the last few days into bring mixed effects models to bear on some behavioural data I've collected for my thesis, but it's occurred to me that I'm not 100% sure that this kind of model is actually appropriate for my data (I only came across them after starting the experiment). In an experiment, `60 participants` completed `28 trials` of a reasoning task, consisting of `14 problems` (call them "A"-"N"), with participants completing each in `2 conditions`, `x` (conflict) and `y` (control).               Participant  Problem               1  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny               2  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny               3  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny               4  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny               5  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny             ...  ...      ...      ...      ...    ...  ...              60  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny      I'm interested in the difference in a trial-by-trial variable (let's call it `reaction time`) between the conflict and control conditions, and would expect it to be higher for the conflict (`y`) trials. Obviously, I would expect to find a within-subject correlation - some subjects are generally fast, some generally slow. I would also expect to find a within- problem correlation - some problems are answered faster than others, regardless of condition. To account for these, I include random intercepts for these two factors: `(1|participant) + (1|problem)`. The difference between conflict and control conditions may or may not turn out to be the same for each subject, and for each problem. For this reason, I consider including random slopes as well: `(1 + condition|participant) + (1 + condition|problem)`. Putting this together, I'm testing a model that looks either like:               null_model = lmer(reaction_time ~ condition(1|participant) + (1|problem), data=data)     condition_model = lmer(reaction_time ~ condition + (1|participant) + (1|problem), data=data)      or               null_model = lmer(reaction_time ~ (1 + condition|participant) + (1 + condition|problem), data=data)     condition_model = lmer(reaction_time ~ condition(1|participant) + (1|problem), data=data)      . **Please; have I horribly misunderstood how this is supposed to work?** **Edit:** The more traditional approach to analysing this data would be to average across problems within each participant, yielding two data points per participant: `conflict condition mean` and `control condition mean`, and then use a paired-samples t test. Reading this question, I thought that this approach should be largely the same as fitting `lmer(reaction_time ~ condition + (1|participant), data=data)`, but trying this in R, it seems otherwise. **Edit #2:** Bounty added.


title: Multidimensional scaling for big dissimilarity matrix

text: I have a large symmetrical dissimilarity matrix of dimension 300 000. Can you please suggest the multidimensional scaling algorithms that can work with such large data? Input of course can be the original vectors from which the dissimilarity matrix was calculated.


title: Hierarchical Multiple Regression vs Ordinal (logistics) regression

text: This is my first research assignment and I feel overwhelmed and confused. I have been informed by my supervisor to do a hierarchical multiple regression using SPSS. I noticed that my dependent variable (DV) is ordinal. Is a hierarchical multiple regression correct or should I perform an ordinal (logistics) regression. Research question: What proportion of variance in health-promoting behaviour is explained by the linear combination of predictors (independent variables specified below): The (DV) = Health-promoting behaviour (Measured by HPLP 4-point Likert scale; 52 items) Independent variables (IV) 1) Worldview: Measured by OMPI 2-point - 26 items 2) Self-efficacy: Questionnaire 4-point Likert - 13 items 3) Psychological distress: HOP-25 4 point Likert - 25 items 4) Gender 5) Study year (First year, Second year, third year) 6) Race (African, Coloured, White, Indian, Other)


title: How to choose the test statistic in Mann-Whitney test?

text: Reading about the Mann-Whitney test for simple random and independent samples I encountered a small issue. According to the book "Introductory Statistics" by Weiss, the test statistic is obtained using $M = \text{sum of the ranks for sample data from population 1}$ As usual, we use this test statistic to decide whether we reject the null hypothesis or not. But this was a bit confusing because it seems arbitrary to choose a given sample as the first one. However, trying to clarifiy this, I found that there are other so-called test statistic $U$, and sometimes we are supposed to choose $\min(U_{1}, U_{2})$ or the opposite $\max(U_{1}, U_{2})$. For example, in this tutorial, this statistic is used: $U_{1} = R_{1} - \frac{n_{1}(n_{1} + 1)}{2}$ where $R_{1}$ is the sum of ranks in population $1$ as above. It also adds: > Note that it doesn't matter which of the two samples is considered sample 1. > The smaller value of U1 and U2 is the one used when consulting significance > tables. But this procedure doesn't seem to be used in Weiss' book. Which one is the correct procedure? Maybe I'm just confusing different tests with similar names.


title: Distribution of logically constrained parameters in Monte Carlo simulation

text: Papers like Briggs et al. 2002 say that logical constraints on inputs such as probability parameters exclude the the Normal distribution from consideration due to its unboundedness. In this example, values below 0 and above 1. They briefly describe the Bayesian derivation from the proportion data (binomial) to the beta, for use as the sampling distribution for probabilities in the probabilistic sensitivity analysis via Monte Carlo simulation. The beta can take on a normal-like shape, but it can also have very non-normal shapes.   1. Can this Bayesian derivation be reconciled with the frequentist central limit theorem that says with a sufficiently large sample size the sampling distribution of the sampling mean will be approximately normal? For resource use (ie. number of visits to the doctor) Briggs et al 2002 use a Gamma distribution, as shaped in Table 4, which is certainly not Normal. Normality is what would be expected for sufficiently large samples in the frequentist interpretation.   2. In a Monte Carlo simulation, if a Normal distribution was used to draw realized values for probabilities, there would be some iterations where values <0 or >1 are observed. However, for symmetric distributions this should cancel out asymptotically. The distribution is still centered around the mean which is within its logical bounds, and Monte Carlo considers aggregate results not results of a single iteration. So I'm wondering why we care about the logical constraints of parameters such as probabilities, costs, etc. within the context of a Monte Carlo simulation?


title: CLT in a Monte Carlo simulation, small sample

text: A CLT says that asymptotically the sampling distribution of the sampling mean converges to the Normal. I would like to run a Monte Carlo simulation using information on one of the model's variables from the literature which reports some mean and standard error, but n=10. I am not comfortable assuming this n is large enough to justify sampling the sampling mean of this parameter using the Normal in the simulation. This is all the information available on the variable parameters. Is there anything that can be said about the sampling distribution of the sampling mean when n is small? What are the implications of assuming the sampling mean follows a normal distribution with small n? Even n=2 gives an unbiased estimate of the population mean and is enough to calculate the standard error. So isn't it the case that we have valid shape parameters for a normal distribution, and even though the distribution at n=2 might not look normal, we can expect the mean and standard error calculated with small n are representative of the mean and standard error with a larger n?


title: Multicollinearity when adding a confounding variable

text: When you run a regression on ice cream sales with predictor shark attacks, you find a significant coefficient. But that is because there is a confounding variable temperature. But how do you correct for this confounder? In case you also add temperature to the model, you have two correlated predictors and the multicollinearity problem. So in case you found the confounder temperature, is the only option to leave the shark attacks out of the model?


title: Plotting results with many features

text: I need an advice for plotting my results. I got different classification results by including different features and different k values. So, I got different results for each of them: My features: alphaRaw16-1, alphaRaw16-2, alphaRaw16-3, grad-x, grad-y, grad-z.   The range of k value: 1 to 21 And, the results look like:               True_Positive   False_Negative  k  Feature        89%              4%          1   alphaRaw16-1        91%              5%          1   alphaRaw16-2        93%              8%          1   alphaRaw16-3        97%              1%          1   grad-x        66%              3%          1   grad-y        77%              42%         1   grad-z        91%              2%          1   grad-x, alphaRaw16-1             ..               ..          ..  ..        94%              7%          1   alphaRaw16-1,alphaRaw16-2,alphaRaw16-3,grad-x,grad-y,grad-x        92%              9%          2   alphaRaw16-1             ...              ..          ..   ..        90%              12%         21   alphaRaw16-1,alphaRaw16-2,alphaRaw16-3,grad-x,grad-y,grad-x      How to best display these results in a meaningful way? Could you please advise me?


title: C++ library to play with statistics (detecting outliers in time series)

text: I am looking for a C++ library for statistics to play with outliers detection in time series (amongst other). What I need:   * Robust estimators, correlations, hypothesis tests, etc;    * No dependencies with external libraries;   * No GPL; Would be a plus:   * Lightweight;   * Free;   * portable;   * Active and supported;


title: Computing platform for R: System76 Linux laptop with dual core i5 or quad core i7?

text: Even though this is not directly a question related to statistical analyses, I thought this group might have some useful insights. My school windows-xp based workstation takes several hours to run simple logistic regression models with bootstrap corrections and bootstrap procedures for model selection (R code adapted from Dr. Frank Harrell's book on logistic regression). Applying and getting a computing account with the University's computing department has been going extremely slowly (it's been over a month since I started the process). I need to get these analyses done as quickly as possible. As a result, I am considering getting a personal linux laptop from System76 (Pangolin). https://www.system76.com/product_info.php?cPath=28&products_id=100 I am hoping to utilize the R multicore package targeted at linux. However, I am not sure about which intel processor I should opt for. The i5 processors have 2 cores and higher clock speed but the i7 processors have 4 cores and lower clock speeds. The price difference between these two groups is around $500. I am really confused about this. Any insights would be immensely helpful. Thanks very much for your time!


title: Why does entropy increase with dispersion for continuous but not for discrete distributions?

text: For a pdf $f(x)$ (i.e. **continuous** distribution), Entropy (differential entropy) is defined as: $H_C(X) = -\int_\mathbb{X} f(x)\log f(x)\,dx.$ For a **discrete** distribution with p.m.f $F(x)$, Entropy is defined as: $H_D(X) = -\sum_{i=1}^n {F(x_i) \log F(x_i)}.$ The definitions look analogous to each other. However, entropy increases with dispersion for continuous but not for discrete distributions. Why?


title: Difference between two slopes

text: I have a dataset with yearly levels of corruption in a number of countries, as well as whether they changed their government that year.               year, corruption, change of president     2001, 5, 0     2002, 7, 1     2003, 8, 0     etc.      I want to test whether corruption is affected by a change in power (defined as the election of a new president who isn't part of the same political party as the previous one). The idea is to either look at the slope of corruption/year for the two years leading into the election, and the two years after, e.g. $t-2$, $t-1$ and $t$ for the before slope. The other idea is to look at the average level of corruption three years before and after. The rate might make more sense since there are more things that affect corruption and different countries may be on different trajectories. However, there could also be some benefit to just looking at average levels three years before and after. Any thoughts on which one I should look at, and also how to go about measuring this?


title: Rigorous definition of an outlier?

text: People often talk about dealing with outliers in statistics. The thing that bothers me about this is that, as far as I can tell, the definition of an outlier is completely subjective. For example, if the true distribution of some random variable is very heavy-tailed or bimodal, any standard visualization or summary statistic for detecting outliers will incorrectly remove parts of the distribution you want to sample from. What is a rigorous definition of an outlier, if one exists, and how can outliers be dealt with without introducing unreasonable amounts of subjectivity into an analysis?


title: Book on designing clinical trials in oncology

text: I would like to have a good idea on how to design clinical trials in oncology. In that issue, I am looking for a compact book that could give me a good overview, with the emphasis on statistical considerations. Would you have a recommendation for me? Thank you in advance, Marco


title: How should you handle cell values equal to zero in a contingency table?

text: How should you deal with a cell value in a contingency table that is equal to zero in statistical calculations? (Note that such a value can be _structural_ , i.e., it must be zero by definition, or _random_ , i.e., it could have been some other value, but zero was observed.)


title: Confusion related to hidden Markov model

text: I am referring to this tutorial: A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition. I am a bit confused about the forward algorithm. From the tutorial: > If we examine the computation involved in the calculation of $\alpha_t(j)$, > 1<=t<=T and 1<=j<=N, it requires the order of $N^2T$ calculations. I didn't exactly get why it is $N^2T$. Further they have compared it to $2TN^T$ when calculating $P(O|\lambda)$. I am also not sure how they can compare calculation complexity for $P(O|\lambda)$ against $\alpha_t(j)$. They are not even same are they?


title: Maximising the correlation between vectors vs minimising the angle between them

text: In this talk http://videolectures.net/lms08_hardoon_scca/ (4:58) David says that maximizing the correlation between vectors can be viewed as minimizing the angle between them, and gives two references: Breiman & Friedman 1985, and Hastie & Tibshirani 1990. The second of these is just their textbook, and the first I can't find, although they had a paper around that time about Generalised Additive Models. Basically I can't find where they discuss this. Is the claim true? Does anyone have a definitive reference?


title: What does "to center two IVs" mean in the context of an interaction?

text: I am just wondering what does "to center two IVs" mean in the context of understanding an interaction. The IVs are A and gender and the DV is B. Is it necessary for the main effects of a new interaction variable (A*Gender) to be a product (e.g. centered main effect of A * centered main effect of Gender) before the IVs can be used in a (hierarchical) regression framework?


title: Sample size to determine if success proportion is zero

text: I have a lot of $N$ devices, and I'd like to know if any of them are faulty. What sample size do I need in order to make sure with $100(1-\alpha)\%$ confidence that none are faulty?


title: Determining the state of a conditional system

text: I have the following situation: There are two systems, $s_A$ and $s_B$. For $s_A$ there are two states: $s_A = a,b$. For $s_B$ there are three states: $s_B = 1,2,3$. Now, what I want to do is measure $s_B$, but this is not possible. Instead, I have the option to measure $s_A$, and I know the following three relations: (if I write $P(a|1)$ I mean that this is the probability of measuring $s_A = a$ given that $s_B = 1$) $P(a|1) = F_1$ $P(a|0) = F_2$ $P(a|-1) = F_3$ In principle I know the values of $F_i$. Moreover, the probabilities all add up to one of course, so $P(b|1) = 1- F_1$ Now from this it is clear that I cannot decide what $s_B$ is from a single measurement. I can however perform the measurement as often as I want to, to build up statistics. However... Every time I measure, there is a chance $p$ that $s_B$ changes. $s_B$ behaves according to a Markov chain, and there are three possible transitions: 1 to 0, 0 to 1, 0 to -1 and -1 to 0. To begin with it would be easiest to just assume that the transition probabilities are all just equal to 1 $p$, but if possible I'd like to do it for four different probabilities as well. So what I want to do is, given all of the above, construct a method of finding $s_B$ from $N$ measurements of $s_A$. Moreover, I also want to find an expression for how well I can do so, and thus find an optimal $N$. Is there anyone who could help me start with this? What I thought of myself is, if I ignore these transition probabilities, N measurements would be like binomial variables. So for example if $s_B = 1$, then the number of measurements of $s_A = a$ should follow the binomial distribution with parameters $N$ and $F_1$. Similarly for the other two of course. I could then find the intersections of these three distributions, and simply see in which interval the number of $s_A$ measurements is most likely to fall. The probability of being right is then just cumulative distribution functions. But I see no way how to adjust this for the probability of $s_B$ changing.


title: Percentile to Z-score in PHP or JAVA

text: I am searching for a way to calculate the $z$-score for a given percentile $p$. I found  this  site doing such a calculation. They have given a formula for calculating $P(Z\leq z)$: $$ \Phi(z) = P(Z \le z) = \int_{-\infty}^z \frac{1}{\sqrt{2\pi}}e^{\frac{-u^2}{2}}du $$ Now I have the value of $P(Z\leq z)$ and I tried solving it to find out $z$, but it was really complex to solve for $z$. Is there any function in JAVA or PHP which can do this? Or is there any mathematical way in which we can solve for $z$? (I am considering hard-coding the $z$ table as my last choice because I would have to make some approximations in that case.) **EDIT** I forgot to mention this, the main problem I am trying to solve is to draw a random point from a normal distribution. I tried using `inversecdf` function of Normal Distribution. Due to the generalized implementation of that function, it is taking lot of time. Say I want to calculate the point at $0.75$ we can do that in a quick way by finding $z$-value at that point followed by a simple math $\bar{x}+\sigma\times z$. So I wanted to find some way to get $z$-score for a given percentile.


title: What does an ACF graph tell me that a PACF graph doesn't?

text: It seems like the PACF graph is more informative than an ACF graph. What information does one get from an ACF graph that can not be gleaned from a PACF graph. (ACF = Auto-correlation function, and PACF = partial Autocorrelation).


title: Odds for my opponent having the Ace in the card game Pitch if I have the King

text: In the card game Pitch using 1 deck with 52 cards and 4 people playing and six cards are dealt to each person and I have a partner and I have a King and wish to bid in the suit the King is in using the King as high what is the percentage my opponents have the Ace in the same suit as my King out of their 12 cards?


title: Using KNN for prediction, how should I normalize my data?

text: Is it better to constrain the data to a range, say [0,1], or to force a mean of 0 and sd of 1? Why? Does the type of input data matter (I'll be using both continuous and categorical variables)?


title: Variance of parameters

text: I have three estimated parameters, $\hat{\beta_0}$,$\hat{\beta_1}$ and $\hat{\beta_2}$, these follow the ordering $\hat{\beta_1}\leq\hat{\beta_0}\leq\hat{\beta_1}$, where these parameters are estimates from a frequentist likelihood function. I know the variance of $\hat{\beta_1}$ and $\hat{\beta_2}$. Can the variance of $\beta_0$ be given in comparison to $\hat{\beta_1}$ and $\hat{\beta_2}$?. By in comparison, I mean like $var(\hat{\beta_1})\leq var(\hat{\beta_0}) \leq var(\hat{\beta_2})$.


title: Exact derivation for finding k-means from Gaussian Mixtures

text: I am having difficulty in deriving k-means from Mixture of Gaussians. I am following the notation from Bishop (2006), Section 9.3.2: Suppose we have : $$ p(\mathbf{x}| \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k ) = \frac{1}{(2\pi \epsilon )^{1/2}} \exp \left( \frac{ -1 }{ 2 \epsilon } \| \mathbf{x} - \boldsymbol{\mu} \|^2 \right) $$ $$ \mathbb{E} \left[ \ln p \left( \mathbf{X, Z} | \boldsymbol{\mu, \Sigma, \pi} \right) \right] = \sum_{n=1}^{N} \sum_{k=1}^{K} \gamma(z_{nk}) \left \lbrace \ln \pi_k + \ln \mathcal{N} ( \mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k ) \right\rbrace $$ ![enter image description here](http://i.stack.imgur.com/qkE5c.png) ![enter image description here](http://i.stack.imgur.com/AW1HE.png) ![enter image description here](http://i.stack.imgur.com/mBFFl.png) It has been claimed that if we take the limit $ \epsilon \rightarrow 0 $ we will have: $$ \mathbb{E} \left[ \ln p \left( \mathbf{X, Z} | \boldsymbol{\mu, \Sigma, \pi} \right) \right] = -\frac{1}{2} \sum_{n=1}^{N} \sum_{k=1}^{K} r_{nk} \| \mathbf{x}_n - \boldsymbol{\mu}_k \|^2 + \text{const.} $$ It is easy to show that $\gamma(z_{nk}) \rightarrow r_{nk}$. But it is not clear to me how to show the above equation. Can anyone help?


title: Is it possible to interpret the bootstrap from a Bayesian perspective?

text: Ok, this is a question that keeps me up at night. **Can the bootstrap procedure be interpreted as approximating some Bayesian procedure (except for the Bayesian bootstrap)?** I really like the Bayesian "interpretation" of statistics which I find nicely coherent and easy to understand. However, I also have a weakness for the bootstrap procedure which is so simple, yet delivers reasonable inferences in many situations. I would be more happy with bootstrapping, however, if I knew that the bootstrap was approximating a posterior distribution in some sense. I know of the "Bayesian bootstrap" (Rubin, 1981), but from my perspective that version of the bootstrap is as problematic as the standard bootstrap. The problem is the really peculiar model assumption that you make, both when doing the classical and the Bayesian bootstrap, that is, the possible values of the distribution are only the values I've already seen. How can these strange model assumptions still yield the very reasonable inferences that bootstrap procedures yield? I have been looking for articles that have investigated the properties of the bootstrap (e.g. Weng, 1989) but I haven't found any clear explanation that I'm happy with. # References Donald B. Rubin (1981). The Bayesian Bootstrap. _Ann. Statist._ Volume 9, Number 1 , 130-134. Chung-Sing Weng (1989). On a Second-Order Asymptotic Property of the Bayesian Bootstrap Mean. _The Annals of Statistics_ , Vol. 17, No. 2 , pp. 705-710.


title: Cross correlation with R

text: I've done a cross correlation but I'm not sure if what I'm doing is correct... Here are the command lines I've used so far (this one is for the first time series. It represents the numbers of subscriptions for men and women on a website during the hours of the day)               data<-read.table(text="hour men women     00h00 475 295     01h00 321 157     02h00 206 127     03h00 141 61     04h00 73 29     05h00 49 22     06h00 71 30     07h00 100 32     08h00 163 55     09h00 219 126     10h00 300 199     11h00 342 255     12h00 407 247     13h00 480 334     14h00 459 358     15h00 481 281     16h00 490 347     17h00 458 309     18h00 599 475     19h00 618 419     20h00 579 453     21h00 565 530     22h00 659 605     23h00 600 435",header=TRUE)        (this one is for the second time series, which represents the time where ads have been shown on TV. I've add the time of the ads within a dummy time series. see this post for why I've done it in such manner : Time series and cross correlation )               hr<-read.table(text="hours     00h00      01h00      02h00      03h00      04h00      05h00      06h00      07h00      08h00      09h00      10h00      11h00      12h00      12h50      13h00      14h00      15h00      16h00      17h00      17h45      18h00      19h00      20h00      21h00      22h00      23h00      23h10",header=TRUE)      (Then I correlate the time series with the command line below)               ccf(data$men,hr,50)      And I'm having this graph (see below) but it doesn't make any sense for me. If you have any clue, I will be glad to hear them. Thanks. ![enter image description here](http://i.stack.imgur.com/mS8aD.png)


title: Scaling mixed models for PCA using dudi.mix

text: I am trying to do a kselect model from the `adehabitatHS` which uses commands from `ade4` package. I am trying to determine if I need to scale my variables. My surface understanding is that the k-select is basically a fancy PCA. In their example they scaled their variables, but their variables were only continuous measures. I have mixed categorical and continuous variables. In their example, they scale their variables before they use the command `dudi.pca` which from my understanding this is needed to set up the k-select. I know from reading the dudi help in the `ade4` vinette that I should use `dudi.mix` here instead of `.pca`, but what to do about scaling? Do I need to scale my variables? Do I need to scale all variables? Do I need to scale all the variables EXCEPT the categorical variable? I can't find reading material that explains what is happening in the dudi/pca mixed process in sufficient detail. Below is the example code from the k-select help if you would like to see what I am referring to.               data(puechabonsp)     locs <- puechabonsp$relocs         map <- puechabonsp$map     pc <- mcp(locs[,"Name"])     hr <- hr.rast(pc, map)     cp <- count.points(locs[,"Name"], map)          ## prepares the data for the kselect analysis     x <- prepksel(map, hr, cp)     tab <- x$tab          ## Example of analysis with two variables: the slope and the elevation.     tab <- tab[,((names(tab) == "Slope")|(names(tab) == "Elevation"))]     tab <- scale(tab)          ## A K-select analysis     acp <- dudi.pca(tab, scannf = FALSE, nf = 2)     kn <- kselect(acp, x$factor, x$weight, scannf = FALSE, nf = 2)      Again, their example only used continuous variables. Does the dudi.mix scale the variables appropriately in the code that we don't see?


title: Demonstration of sample quantile bias

text: While doing some simulations, I realised that the sample quantile is a biased estimator of the true quantile. And, according to my simulations, a potentially _very_ biased one. I was surprised with that result since the empirical CDF isn't biased, but after some internet research, I discovered it was true. I tried to figure out where that bias comes from, but working with sample quantiles is quite difficult. Does anyone have a demonstration of that bias (and, ideally, a quantification)?


title: In principal components regression, should I weight the regression the same as the PCA? Or at all?

text: I am using PCA on foreign exchange return series to find a market "beta". I am using 10 years of daily data with a 2-year half life weighting in the PCA using the package FactoMineR's PCA function. I extract the first principal component return series (so the product of the first eigenvector and the returns matrix) and I want to regress that against each foreign exchange return vector to find the residuals, that is, to find each currency's returns independently of the market beta. Should I use the same 2 year half-life weighting in the regressions? Will this "double up" the weighting somehow? Conversely if I don't weight the regression, will I implicitly be putting too much weight on PC1 returns that are less "relevant"? For what it's worth market participants tend psychologically to put a higher weight on recent than long past currency behaviour. Thanks for the help.


title: MAD equivalent for standard error

text: As far as I know, one can calculate the relative standard error from the standard deviation of a data sample. I am looking for the Median Absolute Deviation equivalent for standard error. Does one exist? Also which methods exist to calculate the relative standard error, that do not require the standard deviation?


title: Linear combination of Gaussian random fields

text: I am modelling three spatial variables, $v_1, v_2$ and $v_3$ (for example, porosity, water saturation and vshale) that are correlated with each other: $\rho_{12}, \rho_{13}$ and $\rho_{23}$. Let's say I have generated three uncorrelated Gaussian random fields, UG1, UG2 and UG3, each with their own spatial correlation structure. My understanding is that the Gaussian random fields can be linearly combined (using the correlation coefficients) to arrive at correlated Gaussian random fields, CG1, CG2, CG3 that preserve the underlying correlation between the variables. I think I understand how this might be accomplished through a cholesky decomposition approach: http://en.wikipedia.org/wiki/Cholesky_decomposition#The_Cholesky.E2.80.93Banachiewicz_and_Cholesky.E2.80.93Crout_algorithms Is there another way to look at this or calculate the correlated Gaussian random fields?


title: Creating model variogram from experimental variogram

text: I have this data consisting of the variance for different distance, I mean experimental variogram. How can I create a model variogram from this in matlab


title: Strong Predictor in Random Forests

text: My statistical learning book asserts that when decision trees are constructed, m out of p predictors are randomly selected for consideration in a split. The idea here apparently is to decorrelate trees obtained in bagging by weakening the effect of a strong predictor. If we have p predictors to choose from, and at each split we can only use m predictors, they say on average, (p - m)/p of the splits will involve the strong predictor. Why?


title: R: Anova and Linear Regression

text: I am new to statistics and I am trying to understand the difference between ANOVA and linear regression. I am using R to explore this. I read various articles about why ANOVA and regression are different but still the same and how the can be visualised etc. I think I am pretty there but one bit is still missing. I understand that ANOVA compares the variance within groups with the variance between groups to determine whether there is or is not a difference between any of the groups tested. (https://controls.engin.umich.edu/wiki/index.php/Factor_analysis_and_ANOVA) For linear regression, I found a post in this forum which says that the same can be tested when we test whether b (slope) = 0. (Why is ANOVA taught / used as if it is a different research methodology compared to linear regression?) For more than two groups I found a website stating: The null hypothesis is: $\text{H}_0: µ_1 = µ_2 = µ_3$ The linear regression model is: $y = b_0 + b_1X_1 + b_2X_2 + e$ The output of the linear regression is, however, then the intercept for one group and the difference to this intercept for the other two groups. (http://www.real-statistics.com/multiple-regression/anova-using-regression/) for me, this looks like that actually the intercepts are compared and not the slopes? Another example where they compare intercepts rather than the slopes can be found here: (http://www.theanalysisfactor.com/why-anova-and-linear-regression- are-the-same-analysis/) I am now struggling to understand what is actually compared in the linear regression? the slopes, the intercepts or both?


title: Establishing relationship between 2 diseases

text: I have the following data of 2 diseases from 5 areas. I want to see if there is any relationship between the 2 diseases. Incidence Rates of 2 diseases (cases per million per year)                   Areas   Disease 1       Disease 2         1       4.653           0.751         2       6.910           1.121         3       4.957           0.745         4       2.870           0.848         5       2.819           1.166          Actual number of cases are     Areas   Disease 1       Disease 2     1       1152            186     2       2601            422     3       1051            158     4        403            119     5        290            120      I am a beginner in Biostatistics. Kindly advise me in simple terms


title: Model comparison in SPSS

text: I have two models with 24 variables in each. I want to see how adding an extra variable (the variable in study) affects the outcome. I used generalised linear mixed models for both models.(GENLINMIXED comman). How do I save the output of the models and compare them? I would like to do what ANOVA in R does. Compare the two models and find the P value to check if there is significant difference between the two models. Many thanks.


title: Measuring linear correlation of non-normally distributed variables

text: How can I measure linear correlation of non-normally distributed variables? Pearson coefficient is not valid for non-normally distributed data, and Spearman's rho does not capture linear correlation. Thank you


title: Dramatic changes in the results from glmm with slight changes in dataset

text: I am running a GLMM in R using the glmmPQL code from the MASS package. I have been getting a convergence error from the following code.               Calvingglmm2<-glmmPQL(Distance~AllRoads+Natural,random=list(~1|Population,~1|Animal),data=dat,family=Gamma)     summary(Calvingglmm2)      So I was fiddling around with my code to see what happened and I started running the following code on a subset of the data, randomly removing just a single value, and seeing what happens with the results.               Calving.Sample <- dat [sample(1:nrow(dat), 72, replace=FALSE),]     Calvingglmm3<-glmmPQL(Distance~AllRoads+Natural,random=list(~1|Population,~1|Animal),data=Calving.Sample,family=Gamma)      summary(Calvingglmm3)      Much to my surprise the results were drastically different each time I ran it. Sometimes I received the same convergence error, sometimes I didn’t. The coefficient scores varied on orders of magnitude between tests, and the p-values ranged from highly significant to nearly 1. The inner workings of the glmm are a bit of a mystery to me despite my best efforts and many many hours of reading. I was just curious if this is a normal behaviour for a glmm or if there is something I have missed in my code writing. Thanks in advance for any tips or thoughts. Also, if you have any tips on getting the original model to converge I would be extremely grateful (I've tried assigning starting values as best I can, but I haven't found anywhere that explains the theory behind this so I'm just guessing).               Dat<-structure(list(Distance = c(1162.19822264933, 933.7389356964,      12464.4096318671, 2302.50341687867, 2308.50538142165, 3915.18043060571,      17837.96671385, 17473.9826616, 4370.44787757821, 9661.80979659471,      494.88702031796, 6343.92286937286, 10330.586481458, 7624.58478308857,      14897.9088190286, 1243.54364758029, 19007.2191276714, 33453.3717670143,      16970.6123001571, 28644.182422185, 320.560246296333, 246.9661127017,      51813.9245697, 7140.31194765857, 45598.15814312, 41031.9354510571,      268.281266841133, 3260.011811305, 8979.62453351428, 4251.30789989014,      4082.86162842957, 5385.99524586143, 1251.47504851657, 4804.98184018143,      14340.5809291429, 20401.4590098429, 2553.975445222, 4041.81979881614,      1484.13677199571, 1169.66034782414, 13245.1513497714, 1126.19983308043,      8259.06344765571, 1161.20714680757, 4181.80076622643, 40775.1222074,      9293.54510051314, 2776.47225473183, 8576.42940441429, 15609.1535519286,      14714.6743856143, 19679.1616246857, 17846.3540124857, 79828.7094106857,      4493.81152102286, 1601.71799450794, 30602.4498660167, 28884.2104902,      4852.31067817286, 143.320137176271, 9140.66714171429, 15083.7211804857,      519.610817095514, 881.757008642143, 104771.790318714, 32624.0004534714,      20764.3259027857, 793.753925335571, 12702.5984742143, 2203.06065660571,      14479.1747836, 9578.73109709, 8055.94785827714), Animal = structure(c(19L,      19L, 20L, 21L, 21L, 21L, 22L, 22L, 22L, 23L, 24L, 25L, 26L, 26L,      27L, 42L, 43L, 44L, 45L, 28L, 29L, 30L, 31L, 31L, 31L, 31L, 32L,      32L, 33L, 33L, 33L, 34L, 34L, 34L, 34L, 35L, 36L, 36L, 36L, 37L,      38L, 39L, 40L, 41L, 1L, 1L, 2L, 3L, 4L, 5L, 5L, 6L, 6L, 7L, 8L,      8L, 9L, 9L, 10L, 10L, 11L, 11L, 12L, 12L, 13L, 13L, 14L, 15L,      16L, 17L, 17L, 18L, 18L), .Label = c("107", "108", "109", "111",      "112", "170", "172", "173", "175", "176", "177", "179", "180",      "181", "182", "183", "184", "188", "Atiko01", "Atiko03", "Atiko14",      "Atiko19", "Atiko20", "BLDVN01", "BLDVN08", "BLDVN10", "BLDVN13",      "GPS05", "GPS06", "GPS11", "OWL07", "OWL10", "OWL11", "OWL17",      "OWL18", "OWL22", "OWL31", "OWL32", "OWL33", "OWL34", "OWL36",      "ROUND02", "ROUND05", "ROUND06", "ROUND09"), class = "factor"),      Year = c(2000L, 2001L, 2001L, 2006L, 2007L, 2008L, 2006L,      2007L, 2008L, 2007L, 2001L, 2003L, 2004L, 2008L, 2008L, 2001L,      2002L, 2002L, 2002L, 1997L, 1997L, 1998L, 2003L, 2004L, 2005L,      2006L, 2004L, 2005L, 2004L, 2005L, 2006L, 2005L, 2006L, 2007L,      2008L, 2005L, 2006L, 2007L, 2008L, 2008L, 2008L, 2008L, 2008L,      2008L, 2010L, 2011L, 2011L, 2011L, 2010L, 2010L, 2011L, 2010L,      2011L, 2011L, 2010L, 2011L, 2010L, 2011L, 2010L, 2011L, 2010L,      2011L, 2010L, 2011L, 2010L, 2011L, 2010L, 2010L, 2010L, 2010L,      2011L, 2010L, 2011L), AllRoads = c(19018.3333333333, 20042.8,      31917.7142857143, 38803.1666666667, 39712, 38754.5714285714,      873, 10783.1666666667, 1703.42857142857, 27205.5714285714,      3968.6, 1033.28571428571, 4008.2, 1444, 1707.14285714286,      6424, 17395.8571428571, 13087.2857142857, 4435.42857142857,      2236.33333333333, 2513.16666666667, 1119.66666666667, 7531.28571428571,      2055.71428571429, 1183.4, 3636.14285714286, 1019.83333333333,      1035.33333333333, 1240, 2483.28571428571, 1368.57142857143,      3534, 694.857142857143, 814, 3269.14285714286, 191.571428571429,      208.142857142857, 632, 236.428571428571, 1029.57142857143,      2427.71428571429, 910, 638.285714285714, 329.142857142857,      52312.8571428571, 51424.7142857143, 14625.5714285714, 8552.5,      294.428571428571, 5578.14285714286, 11982.2857142857, 135.857142857143,      300.428571428571, 6750.71428571429, 13540.1428571429, 15868.2857142857,      4378.16666666667, 200.166666666667, 208.857142857143, 1015.57142857143,      15205.8571428571, 20078.4285714286, 17717.8571428571, 17995,      48803, 4061.85714285714, 2324.85714285714, 5220.14285714286,      721.285714285714, 1610.14285714286, 1974.42857142857, 244.857142857143,      653.857142857143), Natural = c(778, 445.8, 0, 1917.5, 1445.83333333333,      1850, 514.166666666667, 780.333333333333, 611.428571428571,      847.142857142857, 2163, 52.7142857142857, 807.4, 323.714285714286,      268.142857142857, 803.285714285714, 603, 2488.42857142857,      1528.71428571429, 363, 6.16666666666667, 1056.66666666667,      1617.14285714286, 1659, 676.4, 793.571428571429, 742, 790.333333333333,      382.285714285714, 71.7142857142857, 1120, 942, 842.428571428571,      1000.28571428571, 1255, 312.428571428571, 602.142857142857,      877.142857142857, 107.714285714286, 1499.85714285714, 1015.71428571429,      1431.71428571429, 571.285714285714, 1051.71428571429, 6813.71428571429,      5266.28571428571, 13545.1428571429, 13034.8333333333, 4166.42857142857,      5220.42857142857, 0, 2668.28571428571, 5402.57142857143,      4959.14285714286, 5081, 2020.85714285714, 4490, 6531.5, 5872.42857142857,      8955.85714285714, 3121.57142857143, 3854.28571428571, 12779,      12902.4285714286, 55.1428571428571, 10044.4285714286, 2378.85714285714,      5805.71428571429, 1101, 2650.42857142857, 3269.42857142857,      5046.57142857143, 10149.2857142857), Upland = c(11, 116.6,      2240.28571428571, 641.666666666667, 301, 579.571428571429,      205.333333333333, 42.5, 532.285714285714, 3.14285714285714,      508.6, 116.428571428571, 610, 0, 732.285714285714, 1145.85714285714,      54.7142857142857, 355.142857142857, 369.714285714286, 30.3333333333333,      562.5, 18.3333333333333, 34.7142857142857, 187.285714285714,      201.4, 70.2857142857143, 2.66666666666667, 1, 36.5714285714286,      469, 13.4285714285714, 0, 0, 11.8571428571429, 4.14285714285714,      126, 37.2857142857143, 49.8571428571429, 141.714285714286,      74.7142857142857, 14.5714285714286, 0, 156.571428571429,      10.5714285714286, 16587.7142857143, 12820.8571428571, 7608.14285714286,      5334.83333333333, 1260, 415.285714285714, 2083.85714285714,      204.285714285714, 2937.85714285714, 4947.28571428571, 102.857142857143,      355.285714285714, 5985, 21402.6666666667, 1152, 2780.85714285714,      83.2857142857143, 5.57142857142857, 193.142857142857, 357.571428571429,      106.285714285714, 10172.1428571429, 660.714285714286, 4875.28571428571,      13176, 4738.71428571429, 5525.85714285714, 9168.85714285714,      9033.14285714286)), .Names = c("Distance", "Animal", "Year",      "AllRoads", "Natural", "Upland"), row.names = c(NA, 73L), class = "data.frame")


title: Calculating error in sampling

text: I have a certain population U. The members x of U each either satisfy or don't satisfy a property p. That is, p(x) is either true or false for each x in U. I want to get an estimate of the fraction of members of U that satisfy p. I do this by repeatedly sampling from U (with replacement) for some fixed time such as 60 seconds. (This is a computer program I'm talking about.) After sampling n items, I conclude the proportion is the number of items I found for which p is true divided by n. Let f be this fraction (with 0 <= f <= 1). That's all simple enough. What I want to know is whether I can establish (a reasonably tight) upper-bound, b, on the error of f. I naively assumed that the error b would be the abs(f - f') where f' is the fraction calculated after n - 1 samples. From experimentation, I discovered this assumption was wrong. Is there a way to get a reasonably tight bound b?


title: Identifiability of a random effect in an additive model

text: Can you help me to complete the following argument about the identifiability of the random effect $v$ in an additive model. In an additive model such as $X\beta + v$ the location $v$ is unidentifiable since $X\beta + v = (X\beta + a) + (v-a)$ Thank you! **EDIT** (following mpiktas's comment) Does the following make sense to you? Let us assume $v\sim N(\mu, \sigma^2)$. In the previous model, $\mu$ cannot be identified because we can always add some $a$ to the intercept and subtract it to $\mu$, and this does not change the likelihood. In other words, one needs to impose some constraint, e.g., $E(v)=\mu=0$.


title: How to test for statistical significance of regression coefficients in multivariate multiple regression?

text: How can I test the statistical significance of regression coefficients in multivariate multiple regression?


title: How to deal with differences in subgroup analysis but no significant interaction?

text: In case of a multiple linear regression I found a significant effect in a subgroup analysis (sample restricted to males). In the subgroup analysis restricted to females there is no effect. However in a model of the whole sample the interaction term of gender is not significant. I am aware of the fact that significant result in one group and the lack of significant effect in another group does not necessarily mean that there is a significant difference of the association with the dependent variable between the two groups. However I would like to hear your opinion how to interprete such finding. Could one argue that there is limited evidence for a difference in the sense of that evidence is not as strong as it would be in case of a significant interaction effect of gender?


title: How to use variables derived from factor analysis as predictors in logistic regression?

text: ### Context I have a survey that asks 11 questions about self-efficacy. Each question has 3 response options (disagree, agree, strongly agree). Nine questions ask about self-esteem. I have used a factor analysis of the 11 self-efficacy items and extracted two factors. $x_1$ to $x_{11}$ denote the 11 self-efficacy questions in the survey, and $f_1$ ($x_1$ to $x_6$) , $f_2$ ($x_7$ to $x_{11}$) denote the two factors I got from the factor analysis. $y$ is a Dependent variable. Then I created two new variables:                  f1=mean(x1 to x6);         f2=mean(x7-x11).       So the logistic regression would looks like this:                  y=a+bf1+cf2+....      ### My question:   * Can i use these two factors as predictor variables in my multivariate logistic regression model?    * Should I calculate the mean of each items in each factor and use this mean as a continuous variable in my logistic regression model?   * Is this an appropriate use of factor analysis?


title: Accuracy of advanced parametric methods compared to kNN method

text: If kNN doesn't perform well for classification on a dataset, is there any hope for parametric methods to perform better? Kernel-based methods, SVM, random forests, and neural networks. Could any of these outperform kNN method?


title: What are the problems with doing regression on sample means instead of ANOVA in a time series setting?

text: I have data from an annual survey which seeks to assess the needs of consumers on a 1 (least important) to 5 (most important) scale. Thus, I have data that looks like so: ![enter image description here](http://i.stack.imgur.com/0AKaU.png) My goal is to assess if the needs of consumers have changed across years and thus I believe the 'correct' approach to assess if needs have changed is to do a one-way ANOVA with year being the factor. However, another approach is to consider the sample means across years like shown in the table below: ![enter image description here](http://i.stack.imgur.com/enp1O.png) We can then consider a linear regression of RAM vs Years and check if the slope is significantly different from 0. If it is not significantly different from 0 then we can perhaps conclude that the observed variation in customer needs across years is due to random fluctuation and not due to any underlying shifts. What are the drawbacks of the second approach relative to the first? Are they perhaps addressing different questions? Is the second one a weaker test of changes in customer needs as we lose sample information?


title: Advice for transformation in mixed models

text: I use mixed models (random effect for person; repeated measurements) in lme4 to model a continuous outcome. Predictors are both categorical and continous. My outcome variable (blood pressure; typically between 100 and 170) has a non- linear evolution over time, which is presented in the figure below. ![enter image description here](http://i.stack.imgur.com/kZxfu.png) code used:               model <- lmer(pressure ~ sex + treatment + onset_age + follow_up_time + (1|subject_ID), data=mixed6)      ![residual plot](http://i.stack.imgur.com/aD0tu.png) ![qqplot](http://i.stack.imgur.com/EONPT.png) ![the follow_up_time variable \(actually named "duration"\)](http://i.stack.imgur.com/ul1sR.png) As you can see (figures) my model doesn't satisfy linearity. I have tried the following: * centering the response variable * not centering the time variable, i.e using the time variable as it is * using the logarithm of the time variable. * adding the square of the time variable (I(follow_up_time ^2)) to the model. There are >700k observations in the study. 120k individuals. I've heard that with that many observations, many of the assumptions of linearity become less important. Any advice on transformations? Should i convert to generalized mixed models or perhaps the non-linear (nlme?)? I would really appreciate some opinions on this matter.


title: Layers size in convolution NN

text: I have a **pretty complex function** that I'm trying to make my computer learn using CNNs. It involves 70 X 70 grayscaled images. The final output is the output of the last unit (it's because I want to predict probabilities from 0 to 1). What's the best layer architecture to address complexity but at the same time not resulting in overfitting. Any ideas?


title: Gauss-Newton method for MA parameter estimation

text: Please check my solution below for estimating Moving Average parameter using the Gauss-Newton (Linearization) method. I consider MA(1). MA(1) model: $$z_t=a_t-\theta_1a_{t-1}.$$ **Solution:** The residual of this model is, $$a_t=z_t+\theta_1a_{t-1}.$$ The residual sum of squares of this model is $$S(\theta_1)=\sum_{t=1}^Ta_t^2=\sum_{t=1}^T(z_t+\theta_1a_{t-1})^2.$$ We want to minimize the function of the residual sum of squares of this model. And since Moving Average is nonlinear then we use any nonlinear estimation technique, for now we use Gauss-Newton Method. Gauss-Newton Procedure: Let $\theta_1^{(0)}$ be the initial estimate of the model. Then, the residual sum of squares would be, $$S(\theta_1)=\sum_{t=1}^T(z_t+\theta_1^{(0)}a_{t-1})^2.$$ Since the $a_{t-1}$ are not available, then we need to recursively calculate for the $a_t$ using the initial estimate $\theta_1^{(0)}$. That is, $$a_1=z_1+\theta_1^{(0)}a_{0}$$ $$a_1=z_1+0, since\,\, a_{0}=0$$ Then, $$a_2=z_2+\theta_1^{(0)}a_1=z_2+\theta_1^{(0)}z_1$$ and so on. Then after obtaining this values, we can then compute the residual sum of squares. The $E\left[z_t|a_{t-1}\right]$of MA(1) model is $E\left[a_t\right]-\theta_1E\left[a_{t-1}\right]=-\theta_1a_{t-1}.$ Approximating the $E\left[z_t|a_{t-1}\right]$ with a taylor series expansion about its initial guess value. We have, $$E\left[z_t|a_{t-1}\right]=E\left[z_t|a_{t-1}\right]+\left.\frac{\partial E\left[z_t|a_{t-1}\right]}{\partial \theta_{1}}\right|_{\theta_1=\theta_1^{(0)}}(\theta_1-\theta_1^{(0)})$$ Using the general form of the linear regression model: $$Y_t=f(X_t,\boldsymbol\beta)+e_t$$ The, $E\left[Y_t\right]=f(x_t,\boldsymbol{\beta})$ Then, plugging $E\left[z_t|a_{t-1}\right]$ to the model would become, $$Y_t=E\left[z_t|a_{t-1}\right]+\left.\frac{\partial E\left[z_t|a_{t-1}\right]}{\partial \theta_{1}}\right|_{\theta_1=\theta_1^{(0)}}(\theta_1-\theta_1^{(0)})+e_t$$ Then, $$Y_t-E\left[z_t|a_{t-1}\right]=\left.\frac{\partial E\left[z_t|a_{t-1}\right]}{\partial \theta_{1}}\right|_{\theta_1=\theta_1^{(0)}}(\theta_1-\theta_1^{(0)})+e_t$$ Let $\left.\frac{\partial E\left[z_t|a_{t-1}\right]}{\partial \theta_{1}}\right|_{\theta_1=\theta_1^{(0)}}$ be $D_{t1}^{(0)}$ and $(\theta_1-\theta_1^{(0)})$ be $\delta_{1}^{(0)}$. Then this is equivalent to, $$z_t-(-\theta_1^{(0)}a_{t-1})=D_{t1}^{(0)}\delta_{1}^{(0)}+e_t$$ $$Y^{(0)}=D_{t1}^{(0)}\delta_{1}^{(0)}+e_t$$ Where: $Y^{(0)}=z_t-(-\theta_1^{(0)}a_{t-1})$ Hence by ordinary least square, the estimate of the $\delta_1^{(0)}$ is, $$\widehat{\delta}_1^{(0)}=(D_{t1}^{(0)})^{-2}D_{t1}^{(0)}Y^{(0)}.$$ The true estimate would then be, $$\theta_1^{(1)}=\theta_1^{(0)}+\widehat{\delta}_1^{(0)}$$ The $\theta_1^{(1)}$ will not be the final estimate but a new initial guess value for the parameter. This iterative procedure is done until convergence occur. Is this correct? I really appreciate your help.


title: "Continuing" an ARMA process

text: I have a time series and fitted an AR(I)MA object (I really just want an ARMA process so the order is (p,0,q)) in R. Now I'd like to "continue" that process (I guess what I mean is perform simulations). I found arima.sim() but I'm not sure how to get it to use my existing time series (that I'd like to "continue"). I hope my question makes sense! Thanks!


title: R neural net training and prediction

text: I am trying to form a model on a set of data that I gathered from MT4. The OHLC and some MA slopes. I am trying to get the best guess for price change in the future. I am using               # Packages     require(quantmod) #for Lag()     require(nnet)     require(caret)          #Fit model     model <- train(change ~ Lag(change) + high + low + slope13 + slope50 + slope200, mt4,         method='nnet', linout=TRUE, trace = FALSE,                #Grid of tuning parameters to try:                tuneGrid=expand.grid(.size=c(1,5,10),.decay=c(0,0.001,0.1)))      ps <- predict(model, mt4)          #Examine results     model     plot(CHANGE)     lines(ps, col=2)      The output from the model is not so good from what I understand from the RSquared.               6828 samples        6 predictors          No pre-processing     Resampling: Bootstrap (25 reps)           Summary of sample sizes: 6828, 6828, 6828, 6828, 6828, 6828, ...           Resampling results across tuning parameters:            size  decay  RMSE     Rsquared  RMSE SD   Rsquared SD       1     0      0.00483  0.00021   0.000124  0.000262          1     0.001  0.00483  0.000464  0.000124  0.000669          1     0.1    0.00483  0.00115   0.000125  0.0016            5     0      0.00499  0.000452  0.00074   0.000537          5     0.001  0.00484  0.000271  0.000123  0.000356          5     0.1    0.00485  0.00109   0.000131  0.00183           10    0      0.00485  0.00087   0.00012   0.00129           10    0.001  0.00486  0.000408  0.000121  0.000699          10    0.1    0.005    0.000638  0.000434  0.00121       Questions:   1. Am I training my model incorrectly because of such the bad output?    2. How do I get the best guess for the next bars change that I am training for? Even though this model appears horrible it would be nice to know how to get the output.


title: Variance decomposition in linear regression model

text: Consider the linear model $y = \mathbf{X}\mathbf{\beta} + \epsilon$. The residual variance-covariance matrix is given by $\text{Var}(\epsilon)$. Greene's textbook* states that: $$Var(\epsilon) = E[Var[\epsilon|\mathbf{X}]] + Var[E[\epsilon|\mathbf{X}]]$$ **Why is this the case?** Also, **why is this useful**? * * * *Page 62 in Econometric Analysis by William H. Greene


title: Controlling for variables in pooled OLS estimation in EViews

text: I am working on the Chinese economy and my topic of research is how external political instability can affect Chinese exports. So I want to estimate the Chinese export demand function for 1988-2011 with more than 130 countries. I want to estimate the regression equation given below. $$ \begin{align} \log(\mathrm{export})_{it} &= \beta_0+ \beta_1 \log(\mathrm{real gdp})_{it}+ \beta_2 \log(\mathrm{population})_{it} \\ &\quad+ \beta_3\mathrm{political stability}_{it}+ \beta_4\mathrm{realexchange rate}_{it}+ \varepsilon_{it} \end{align} $$ Where $\log(\mathrm{export})_{it}$ is the (log) of the export from china to other countries and $t=1988, \ldots, 2011$. According to economic theory, the export of a country depends not only domestic GDP and population but also on the GDP and population of other countries. In my research, I want to control for the effect of Chinese GDP and population on Chinese export in EViews in a pooled OLS estimation but I don't know how to do this. If i do not control these two variables I can run a pooled OLS estimation in EViews. Any help is greatly appreciated.


title: How to know dispersion if $\mu$ is close to or below 0 (chance-corrected beta-binomial model)

text: ### Background In sensory science, "replication" means having a panelist in a taste panel do multiple rounds of the same test. You cannot just count those additional rounds as additional panelists, because panelists might show consistent performance. In that case, replications would just confirm the abilities of individual panelists rather than provide more information about the difference between the two stimulants. So, theory suggests using a corrected beta- binomial model in which the abilities of the panelists are assumed to be beta- distributed, not between 0 and 1, but between $C$ (the guessing probability of the taste test used) and 1: $$ P(x|n,a,b,C)=\frac{(1-C)^n}{B(a,b)}\binom{n}{x}\sum_{i=0}^{x}\binom{x}{i}\left(\frac{C}{1-C}\right)^{x-i}B(a+i,n+b-x) $$ $C$ is the guessing probability of the test method used (in my case $C=\frac{1}{3}$),   $n$ is the number of replications (in my case $n=2$ for all panelists),   $x$ is the number of successes. Log likelihood: $$ L=\sum_{i=1}^{k}log[Pr(x_i|n_i,\mu,\gamma)] $$ $k$ is the number of panelists (in my case $k=54$),   $n_i$ and $x_i$ are the number of replications and successes of the $i$th panelist,   $\mu$ is the mean,   $\gamma$ is the over-dispersion parameter (sometimes also referred to as $\rho$). $a$ and $b$ can be calculated from $\mu$ and $\gamma$: $$ a=\mu\left(\frac{1}{\gamma}-1\right)\text{ and }b=(1-\mu)\left(\frac{1}{\gamma}-1\right) $$ (Bi, Sensory discrimination tests and measurements, 2006, as shown in Gacula et al., Statistical Methods in Food and Consumer Research, 2009 § 10.3. See also Brockhoff, The statistical power of replications in difference tests, 2003 and Næs et al., Statistics for Sensory and Consumer Science, 2010 § 7.7.) ### My test results I had 54 panelists perform 2 replications of a sensory test with a guessing probability of 33% (n=2, k=54, C=1/3), all comparing the same two pancake recipes. Depending on the taste difference between the pancake recipes, you would expect the observed proportion of successes $p_c$ to run up from 33% in case of no taste difference to 100% in case of obvious taste difference. However, my test results show panelists correctly grouped the pancakes in only 29 of the 108 rounds (27%). Here is the R-code for my data:               round1 <- c(0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1,          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1,          1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0)     round2 <- c(1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1,          0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,          0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0)     x <- round1 + round2     n <- c(rep(2,length(x)))     dat <- data.frame(id = factor(1:length(x)), x, n)      I have no reason to assume the panelists gave wrong answers on purpose, nor that I mixed up samples, so I guess it was just chance that led to this 'low' performance. ### Maximum Likelihood Estimation As described in Gacula (see link above), I created an R-function of the above- mentioned log-likelihood of the corrected beta-binomial (CBB) model, called it llcbb, then ran nlminb to find MLE-estimates of $\mu$ and $\gamma$. The estimate for $\mu$ is close to 0 and the estimate for $\gamma$ will usually stay near its start value:               fit<-nlminb(start=c(.1, .9), function(x) -1*llcbb(dat,mu=x[1],gamma=x[2],1/3)         ,control=list(trace=1))     fit$par     # [1] 5.986091e-13 8.974503e-01          fit<-nlminb(start=c(.1, .1), function(x) -1*llcbb(dat,mu=x[1],gamma=x[2],1/3)         ,control=list(trace=1))     fit$par     # [1] 3.040535e-12 9.976926e-02      The following picture shows my findings in black and some possible configurations of the CBB-model. The cyan line shows that when $\mu$ approaches 0, it does not really matter anymore what $\gamma$ is. I guess this explains why nlminb leaves $\gamma$ untouched. ![observed results and CBB models](http://i.stack.imgur.com/OPwa4.png) ### Moment estimation Using other formulas provided by Gacula (see link above) I calculated the following moment estimates: $$ \hat{\mu}=-0.097222\text{, }\hat{\gamma}= 0.349006 $$ Wikipedia states the following about negative moment estimates: > Note that these estimates can be non-sensically negative which is evidence > that the data is either undispersed or underdispersed relative to the > binomial distribution. In this case, the binomial distribution and the > hypergeometric distribution are alternative candidates respectively. (Beta-binomial distribution. (2014, February 7). In Wikipedia, The Free Encyclopedia. Retrieved 14:42, February 14, 2014, from http://en.wikipedia.org/w/index.php?title=Beta- binomial_distribution&oldid=594290522) I understand the suggestion to revert to the binomial distribution, but even then I still feel I need to be sure about the dispersion, because the upper limit of the confidence interval on $p_c$ varies quite a bit depending on whether I can use the replication data (in case of no dispersion) or not (in case of over-dispersion):               binom.test(29,108,p=1/3,alternative="less")     # 95 percent confidence interval:     #  0.0000000 0.3475966          binom.test(15,54,p=1/3,alternative="less")     # 95 percent confidence interval:     #  0.0000000 0.3950352      ### Question How to deal with possible dispersion ($\gamma$ or $\rho$) in beta-binomial distribution in case $\mu$ is estimated to approach 0, or to be below 0?


title: Binary logit with time series

text: I intend to run a Binary Logit with Wholesale Price Index of different commodities. I have converted these scale variables into categorical variables with the following criterion: 1=Price has increased when compared to previous period 0=Price has remained stable or decreased in comparison to previous period With the above criterion, all predictors have been converted into categorical variables. Could I then estimate using the usual Binary Logit option in SPSS or should I still make some correction to the data or adopt a different procedure? My aim is to study the inter-dynamics between sectors....Eg: With increase in the price of iron, the odds that there would be an increase in the price of Cement. I would be grateful if anyone could guide me in this regard. Thanks!


title: Do Bayesians interpret the likelihood distribution as subjective as well?

text: One of the main differences between Bayesians and frequentists is that they have a subjective interpretation to probability. However, do Bayesians actually interpret subjectively the probabilities attached to an outcome GIVEN a set of parameters (i.e. for the likelihood), or is it just that they attach a subjective probability to the prior, and also to the posterior as a consequence? (but $p(x | \theta)$ is thought similarly to the way frequentists think about it.)


title: Help understanding an explanation about minimum description length principle

text: I'm reading about MDL principle and my problem is on a book called: "Guide to intelligent data analysis", 2nd edition page 106. I have added here a picture of the page I'm having trouble with. Could someone decode what the author is saying in the highlighted area: ![enter image description here](http://i.stack.imgur.com/sQzQ1.png) It seems abstract and confusing...I get confused when he says that the original data is contained in the decompression rule itself?!...huh?... and because of that I don't understand the rest of it :/ Thank you for any help


title: Can anyone tell me why we always use the Gaussian distribution in Machine learning?

text: For example, we always assumed that the data or signal error is a Gaussian distribution? why? I have asked this question on stackoverflow, the link: http://stackoverflow.com/questions/12616406/anyone-can-tell-me-why-we-always- use-the-gaussian-distribution-in-machine-learni


title: k nearest neighbor with decision tree

text: A dataset has a few attributes. One of the attributes(attribute X) represents a distance with values expressed in meters. I use cross validation to estimate the performance of Decision tree and k nearest neighbor classifiers. My questions 1) Why will the performance of k-nn change if I change the representation of attribute X to cm instead of m ? 2) Why will the performance of k-nn and decision tree not change if I multiply all the attributes with 20?


title: How do I find a 95% percent ML confidence interval of the probability of category 5 strength

text: Let ${\bf X} = (X_1, \dots, X_n)$ denote the log-intensities of all $n = 291$ tropical cyclones from the north Atlantic basin between 1980 and 2006. Log- intensity is defined as the logarithm of its maximum sustained windspeed (measured in knots). A tropical cyclone is a category 5 hurricane if its maximum sustained windspeed is 137 knots or more, i.e. if its log-intensity is 4.92 or more. We are interested in the probability $p \in [0,1]$ of a tropical cyclone attaining category 5 strength. Consider the following model: $X_i \overset{{\rm iid}}{\sim} \mathsf{Normal}(\mu, \sigma^2)$, where $\sigma = 0.33, \mu \in (-\infty, \infty)$. Under this model $p = \Phi\left({\mu - 4.92 \over \sigma}\right)$ where $\Phi(x)$ denotes the standard normal CDF. Find the 95% ML confidence interval for $p$ based on the following summaries of the data $$\bar x = 4.21 \quad \sqrt{n} = 17.06 \quad \#\{x_i \ge 4.92\} = 9$$ * * * **My Work** The formula for an ML interval for $X_i \overset{{\rm iid}}{\sim} \mathsf{Normal}(\mu, \sigma^2)$ with known $\sigma$ is $B_c(x) = [\bar X \mp c{\sigma \over \sqrt{n}}]$. Calculating along these lines, for $c = 1.96$ as the $z$ value that leaves $0.05$ probability in the two tails combined of the normal pdf, I get $B_{1.96}(x) = [4.172, 4.247]$, but when I try to evaluate the endpoints of this interval for $p$ given that $\Phi(x)$ is a monotone transformation, I get nonsense answers. What am I doing wrong? Also, this will never be handed in for any class, I simply came across this problem and wondered how it could be solved.