title: How to choose the best of two highly correlated predictors in cox proportional hazards regression

text: I have built a good model of time-to-stroke under cox ph assumptions using a predictor of stroke risk (Framingham risk score). It incorporates a score according to Age, Gender, controlled / uncontrolled bp, cardiovascular risk factors and can be seen here : http://www.framinghamheartstudy.org/risk/index.html I have hypothesised that it will be a predictor of time to event in my population. However, it incorporates age, in fact it is correlated 85% +/-3% with age. So my question is: How do I effectively assess whether it is a _better_ predictor than age alone. (It's pretty useless clinically if it's not) My intuition has been to include age alone then the framingham score alone in the model with the controlling predictors and then compare the AIC between the two fits, choosing the one which changes the AIC most - compared to the model with only the controlling predictors in it. Is this a good solution to this in general? Or is it just plain wrong.


title: What is the distribution of the euclidean distance between two random points in 2d space?

text: I would like to know the distribution of z as the euclidean distance between 2 points which are not centred in the origin. If I assume 2 points in the 2D plane $A = (X_a,Y_a)$ and $B = (X_b,Y_b)$, where the $X_a \sim \mathcal N(\mu_a,s^2)$, $X_b\sim\mathcal N(\mu_b,s^2)$, $Y_a\sim\mathcal N(\nu_a,s^2)$, $Y_b\sim\mathcal N(\nu_b,s^2)$, then the Euclidean distance between $A$ and $B$, would be $z= \sqrt{(X_a-X_b)^2 + (Y_a-Y_b)^2}$. Now: $X=X_a-X_b$ and $Y=Y_a-Y_b$ are themselves random variables with means $(\mu_b-\mu_a)$ and $(\nu_b-\nu_a)$ and variance $2s^2$, so the problem that I have is determining the pdf of $z=\sqrt{X^2 +Y^2}$, knowing that $X$ and $Y$ are 2 uncorrelated Gaussian random variables with **non-zero means** and the same variance, $2s^2$. The Rician distribution applies when $z$ is the distance from the origin to a bivariate random variable. This has been proven only when the random variables ($A$ and $B$) are **circular** bivariate random variables (A proof can be found in L. C. Andres and R. L. Phillips, _Mathematical Techniques for Engineers and Scientists_ , 2003, Ch. 13, Sec. 13.8.2, p. 680). I would like to know the pdf/cdf of $z$ as a distance between two points (none of them being centred in the origin) when they are not circular. Is there a known parametric distribution for $z$? What would this distribution look like if it is a generalized form of the Rician distribution? Thank you!


title: Kendall's tau and independence

text: Kendall's tau is an association measure between two random variables, say $X$ and $Y$. If $X$ and $Y$ are independent, then $\tau = 0$. However, $\tau = 0$ does not imply independence because kendall's tau only captures monotonic association. Is there some conditions under which one has equivalence? _Edit to account for @Cardinal's comment_ Say that $X > 0$ and $Y > 0$.


title: Negative number of parameters in hierarchical bayesian model

text: I'm using Deviance information criterion to assess the fitness in my Bayesian hierarchical model. The functional form of this criterion is as follows: $$DIC=p_{D}+\bar{D}$$ where $p_{D}=\bar{D(\theta)}+D(\bar{\theta})$ is effective number of parameters and $D(\theta)=-2log(p(y|\theta))$, $p(y|\theta)$ is sampling distribution and $\theta$ are unknown parameters. And the problem is that the effective number of parameters when I use expectation as an estimes of $\theta$ is negative, however if I plug in a medians as estimates of parameters - $p_{D}$ is positive. I know that negative $p_{D}$ indicates prior-data conflict in non-log-concave distributions, but the the visual inspection of the fitted values indicates very good fit. So, is it ok to use median as paramete estimates instead of expectations in DIC?


title: Finding a confidence interval for a parameter

text: > Let $X_1, ... , X_n$ be iid with common density $f(x) = ax^{a-1}$, $ 0 < x < > 1$, where $a$ is an unknown parameter $a > 0$. > > (a) Find the MLE for $a.$ > > (b) Show that $W = - \sum^n_{i=1} \ln X_i$ has a Gamma distribution with > parameters $\alpha = n$, $\beta = 1/a$. Hint: Use the mgf. > > (c) Show that $2aW$ has $\chi^2 (2n)$ distribution. > > (d) Use (c) to find a $100(1-a)%$ confidence interval for $a$. I'm a bit unsure about part (d). Let's say the critical values are $\chi^2_{1 - \alpha/2}$ and $\chi^2_{\alpha/2}$. Then we have $$ (1 - \alpha) = P(\chi^2_{1-\alpha/2} < 2aW < \chi^2_{\alpha/2}) = P \left(\frac{\chi^2_{1-\alpha /2}}{2W} < a < \frac{\chi^2_{\alpha/2}}{2W} \right)$$ But I'm not sure if this works, since I'm dividing by a random variable. So is this correct? If not, can you give me a hint? Thanks in advance


title: Automap or Gstat Question about how to fit the variogram

text: I'm working with autoKrige and Gstat. I was wondering if it's possible to add to two differents models like Exp with Gau. > I know if I put something like this > "kriging_result<-autoKrige(z~1,data,grid,model="Sph",fix.values = > c(0,NA,NA))" the program uses that model to fit the variogram, but how can I > mixed two models? > > on the other hand, in gstat I know if i use fit.variogram the variogram is > fit automaticly but if I use something like this g1<-vgm(60000,"Sph",1500,0, > add.to = vgm(3000, "Exp", 60, nugget = 0)) fit<-fit.variogram(vario,vg1, > fit.sills = TRUE, fit.ranges = TRUE, fit.method = 7, debug.level = 1, > warn.if.neg = FALSE ) it doesn't work So, how can i do this works? Thanks > for the help


title: Conditional heteroskedasticity/variance and uncertainty in estimated residuals

text: Say you've got a simple cross-sectional model $$ y=\alpha + \beta T +\epsilon $$ where $T$ is binary. This is basically a t-test, but lets treat it like a regression. You're interested in the effect of $T$ on the first and second moments of the outcome. (You might have done GLS or IRLS as well, though my question should be orthogonal to that choice.) You might fit a model like $$ g(\hat\epsilon^2) = \alpha +\delta T + \eta $$ If you assume that $\epsilon = \hat\epsilon$, then your estimate of $\delta$ is fine (or as fine as your estimate of the mean model is). But $\hat\epsilon$ is only as good as the estimates of the mean model, which might be highly uncertain. Are there standard ways of dealing with this? I haven't studied time series econometrics in depth, but glancing at definitions of ARCH/GARCH/etc I don't get the impression that uncertainty in the residuals is taken into account when trying to estimate the effect of the independent variables on the variance of the outcome. I can imagine sampling from the multivariate-normally-distributed parameter vector to get many estimates of the mean model, and fitting the variance model many times. But then how would you combine the uncertainty in the residuals with uncertainty in the variance model coefficients at each iteration? Perhaps you'd treat it like a multiple imputation problem, treating the residuals as missing data? But anyway that whole approach would be quite computationally intensive. EDIT: and it would be especially computationally intensive if you were doing IRLS and then did so at every single IRLS iteration. Suggestions for previous work on this would be very much appreciated.


title: Confusion related to a derivation in a paper related to gaussian

text: I was reading this paper where it is mentioned Let F(·) be the CDF of a continuous random variable X and $\phi^{−1}(·)$ be the inverse of the CDF of N(0,1). Consider the transformation from X to Z by Z = $\phi^{−1}(F(X))$. Then it is easy to see that Z is standard normal regardless of F I didn't get how it is true. What does $\phi^{−1}(F(X))$ denote here. I am a bit confused. Explanations guys help me out?


title: GEE: Pairwise comparisons at different levels of a covariate?

text: I'm familiar with basic regression methods, but have no experience using GEEs. I use SPSS, and I'm trying to use a GEE for a dataset that I have, because there is a repeated measures component in my model. Anyway, in my model, I have 3 dichotomous predictors (let's say, $x_1$, $x_2$, and $x_3$), as well as one continuous covariate ($x_4$). My output tells me that I have a 4-way interaction, and just by looking at my plots, it seems that this 4-way interaction is driven by a 3-way interaction at one level of my covariate, but not another. That is, that there is a 3-way interaction between my 3 dichotomous predictors at high levels of $x_4$, but not at low levels. If this were simple OLS regression, I would just center the model at low vs. high levels of the covariate, and test for the interaction at these two levels. However, I'm not sure how I would do this using a GEE. Specifically, when I ask for pairwise comparisons, the model always centers itself at the mean of the covariate. So, I have two questions: (a) is it possible to conduct pairwise comparisons at 2 levels of the covariate (e.g., at the 25th and 75th percentiles)?, if not (b) how should I go about probing this 4-way interaction? Thank you!


title: Basic time series model question

text: If I fit a time series with a AR model, I can tell if this fit is valid by checking the residual correlation (valid if not correlated). But if I fit it with a MA model, how do I tell if this fit is valid or not? In R, when I specify to fit the time series with, say MA(1), it will automatically use white noise for epsilon, right? $y_t = \epsilon_t + \text{coefficient} \cdot \epsilon_{t-1}$ Thanks


title: Relationship between estimates in 2 linear models

text: I am considering two linear models with model matrices related as follows: \begin{equation*} \text{Model} \ 1 \\\ \begin{bmatrix} \mu_1 \\\ \mu_2 \\\ \vdots \\\ \mu_n \end{bmatrix} = \begin{bmatrix} 1 & x_{12} & x_{13} \\\ 1 & x_{22} & x_{23} \\\ \vdots & \vdots & \vdots \\\ 1 & x_{n2} & x_{n3} \end{bmatrix} \begin{bmatrix} \beta_0 \\\ \beta_1 \\\ \beta_2 \\\ \end{bmatrix} . \end{equation*} \begin{equation*} \text{Model} \ 2 \\\ \begin{bmatrix} \mu_1 \\\ \mu_2 \\\ \vdots \\\ \mu_n \end{bmatrix} = \begin{bmatrix} 1 & cx_{12} + x_{13} & x_{13} \\\ 1 & cx_{22} + x_{23}& x_{23} \\\ \vdots & \vdots & \vdots \\\ 1 & cx_{n2} + x_{n3} & x_{n3} \end{bmatrix} \begin{bmatrix} \beta^*_0 \\\ \beta^*_1 \\\ \beta^*_2 \\\ \end{bmatrix} . \end{equation*} $$y\sim N(\mu,I_n \sigma^2)$$ Under what conditions is $\frac{\hat\beta_1}{\hat\sigma_{\beta_1}} = \frac{\hat\beta^*_1}{\hat\sigma_{\beta^*_1}}$, and how can it be proven?


title: Several logistic regressions vs multinomial regression

text: Is it viable to do several binary logistic regressions instead of doing a multinomial regression? From this question: Multinomial logistic regression vs binary logistic regression I see that the multinomial regression might have lower standard errors. However, the package I would like to utilize has not been generalized to multinomial regression (`ncvreg`: http://cran.r-project.org/web/packages/ncvreg/ncvreg.pdf) and so I was wondering if I could simply do several binary logistic regressions instead.


title: Comparing categorized data

text: I am quite rusty on my stats beyond standard deviation and linear regression, so I am not even sure about how to phrase this question. I am looking at a long history of credit card qualification data, so I know how much money has been processed as rewards cards, AMEX, World cards, etc. I have this data for each month. Now I need a way to test if one month of data is unusual, statistically speaking. For example, in January, 21% of the volume was in Rewards, and in February, only 19% was in rewards, an so on through at least a year. I want to know if the volume in Rewards was unprecedented the next January. Can anyone provide a name of a test I can look at to get an answer?


title: Joint distributions for uncorrelated varibles

text: Can someone think of joint distribution of random variables X, Y such that the following three conditions are satisfied: $E[X] = 1$,   $E[Y] = 1$, and   $E[XY] = 0$? A friend of mine asked me this, which seems simple, but I'm not a stat whiz, so This makes my head spin. Thank you for answering.


title: Calculating mean percents

text: Are there any rules to calculating mean percent or is it the same as calculating any other average (i.e., add all of the percents together and divide by the sample size)?


title: Two general questions about log-normal distributions: shift and mixture

text: I have the two following questions:   1. Imagine I have two sets of observations, and both sets have a lognormal distribution. Now, given I look at the union of the two sets - is the distribution still a lognormal?   2. Given now I have again a set of observations with lognormal distribution. Given I add or subtract a scalar to $LogN(\mu, \sigma^2)$, what happens to each observation. Let's say I have a right shift of +2. Does that imply, that the value of each observation is bigger by 2? Thanks


title: How would we get the conditional distribution?

text: Having the marginal distributions, say $f(x)$ and $f(y)$, how would we get the conditional distribution $f(x|y)$? The relation is given by: $$f(x)=\int f(x|y)f(y)dy$$ Do we need to find the derivatives of both sides and then solve a differential equation?!


title: Coda output, joint quantile

text: I am running a hierarchical model in JAGS. It is a piecewise constant survival regression model. So I have 4 sets of regression coefficients and baseline hazard rates. I would like to plot a mean curve and range using the mean and 95% highest density interval (HDI) of these parameters using the covariates from my data. I have output from CODA, and it provides marginal highest density intervals for each parameter. What would be the correct way to get the joint HDI for all parameters?


title: HMM forward algorithm in MATLAB

text: Does anyone know where can I find a pseudo code or MATLAB code of the HMM forward algorithm?


title: What are the pros and cons of learning about a distribution algorithmically (simulations) versus mathematically?

text: What are the pros and cons of learning about a distribution's properties algorithmically (via computer simulations) versus mathematically? It seems like computer simulations can be an alternative learning method, especially for those new students who do not feel strong in calculus. Also it seems that coding simulations can offer an earlier and more intuitive grasp of the concept of a distribution.


title: Obtaining standard error on a data point obtained from linear regression

text: I have data with standard error, included below for clarity,                   X           Y           Error in Y         0.0105574   -28.831027  0.04422         0.0070382   -27.800385  0.04225         0.0052787   -27.314088  0.04209         0.0042229   -27.054207  0.04185         0.0035191   -27.000188  0.04143         0.0030164   -26.891275  0.04108      I have obtained parameters **a** and **b** of the expression **y=a*x*x + b** from a weighted least squares regression using this data ( _fit_ in gnuplot). The regression returned what was called "Asymptotic Standard Error" associated with these parameters. I believe this error was calculated using the deviation from fitted point to actual points (Equation 34/35 here) and is used to assess the quality of a fit. However, this is not the error that I'm interested in. I'm looking to determine the value of the data point at X=0.0 from my fitted function _with standard error_ like my other values. The output of the regression was:                   Final set of parameters            Asymptotic Standard Error         a              = -19389.1         +/- 752          (3.878%)                 b              = -26.7951         +/- 0.03915      (0.1461%)      So, to be quite specific, **how might I determine the standard error at the point (X,Y)=(0.0, -26.7951) using my fitted function**? I expect the error in this calculated point to be much larger than the errors of the values reported in the Y values of my table above. I can see how gnuplot is not the right tool for this, as it only weights my data points using the standard error in my input. What I need to do is propagate the error in my original data points to obtain the error on the regression line. This seems like a pretty basic exercise, sorry for my statistics ignorance. Thanks!


title: Specific estimates about fuel consumption - Simple linear regression?

text: I am working on a story problem for a project. "You work for a small environmental foundation that wants to analyze fuel consumption. Your boss (in the year 2002) has asked you to help her analyze the 2001 data on highway fuel consumption to find out what a change in the tax rate does to fuel consumption. She has provided a datafile (labelled `Fuel2001`) that has data on the following variables: `Drivers` Number of Licensed drivers in the state `FuelC` Gasoline sold for road use (1000s of gal.) `Income` Per capita personal income (year 2000) `Miles` Miles of Federal-aid highway miles in the state `MPC` Estimated miles driven per capita `Pop` Population age 16 and over `Tax` Gasoline state tax rate, cents per gallon She wants you to produce an analysis that does the best job of helping her argue that higher tax rates lead to lower fuel consumption. But her initial analysis of the data (simply including all of these variables in a naive regression model) did not find that tax rates had a statistically significant effect. Is there anything that can be done to find out if there is a relationship? Write her a short memo on that subject including your analysis, including specific estimates of the relationship." I'm supposed to present her with the best case possible to support her cause (raising gas tax to reduce fuel consumption), but I should also include the best model, for her information. I have transformed all the skewed variables by taking their natural log. I ran a bunch of different regressions in Stata until I found one in which tax was statistically significant. When I ran the same regression robust tax was no longer statistically significant. What would you do and show (graphically, including tables or graphs)? What does it mean that the robust regression makes tax insignificant? I know it would be helpful to have the data, but I can't figure out how to upload it. Thanks in advance for any input.


title: Long tailed distributions for generating random numbers with parameters to control tail heaviness

text: I have to generate random numbers for my algorithm based on probability distributions. I want a distribution which has heavy tails and is unskewed, which can produce numbers far away from location parameter. There should be a parameter to control the tail heaviness (e.g., like levy distribution where alpha determines tail heaviness). I have identified the t-distribution (for smaller degrees of freedom) and the laplace distribution as two possibilities.   * Are there any reasons to prefer t or laplace for my purpose?   * Apart from t-distribution and laplace distribution, is there any distribution except cauchy or levy that would be useful for my purpose?


title: Weka: Test set not getting processed in Knowledge Flow? ClassAssigner fails?

text: I am trying to predict labels for a test set (around 15000 instances of around 400 features) based on a model trained with a training set (around 86 instances of around 400 features). I am using Weka 3.6.11 on a humble Windows laptop. Both training and test set are loaded from a .csv file (one file for each). The test set did not have labels (this is for prediction, not validation), so I added a colum of NaNs for this attribute in the .csv file using Matlab. (At first I tried without this step and encountered errors, I read on forums that this was the proper setup). This is my first time using Weka for predicion. I tried two different setup in the KnowledgeFlow GUI, and failed both times. **SETUP 1: Use ClassAssigner within FilteredClassifier**                   filter: weka.filters.MultiFilter           weka.filters.unsupervised.attribute.ClassAssigner <- Class is assigned here (last         columns)           weka.filters.unsupervised.attribute.NumericToNominal           weka.filters.unsupervised.attribute.Remove           weka.filters.unsupervised.attribute.Remove           weka.filters.unsupervised.attribute.Normalize         classifier: weka.classifiers.bayes.NaiveBayes         execution slots: 2 (this is the default, not sure what it means)      Then what I did was manually click on "start loading" for the training set, and then "start loading" for the test set. This resulted in the Classifier getting interrupted during training, and then during test an error is thrown: No attribute set in test data. I don't understand this error since "ClassAssigner" is the first module to run after "TestSetMaker". Here is the error log if that helps:                   22:08:11: [Loader]  CSVLoader$22564265|-M ? -E ",'|loaded   filename_training             22:08:11: [Classifier] FilteredClassifier$24309584|-F "weka.filters.MultiFilter -F                         \"weka.filters.unsupervised.attribute.ClassAssigner -C last\" -F         \"weka.filters.unsupervised.attribute.NumericToNominal -R last\" -F         \"weka.filters.unsupervised.attribute.Remove -R first\" -F         \"weka.filters.unsupervised.attribute.Remove -V -R numbers \" -F         \"weka.filters.unsupervised.attribute.Normalize -S 1.0 -T 0.0\"" -W          weka.classifiers.bayes.NaiveBayes --|starting executor pool ( 2slots)...          22:08:11: [Classifier]  [Same information as above] --|setup output queues.           22:08:11: [Classifier]   [Same information as above]  --|scheduling run  1fold   1for execution...         22:08:11: [Classifier] [Same information as above] --|ERROR: no class attribute set in test data!         22:08:11: [Classifier] [Same information as above] --|( run  1fold1) interrupted!         22:08:21: [Loader]  CSVLoader$11661342|-M ? -E ",'|loaded   filename_test             22:08:21: [Classifier]  FilteredClassifier$24309584|-F "weka.filters.MultiFilter -F \"weka.filters.unsupervised.attribute.ClassAssigner -C last\" -F \"weka.filters.unsupervised.attribute.NumericToNominal -R last\" -F \"weka.filters.unsupervised.attribute.Remove -R first\" -F \"weka.filters.unsupervised.attribute.Remove -V -R numbers \" -F \"weka.filters.unsupervised.attribute.Normalize -S 1.0 -T 0.0\"" -W weka.classifiers.bayes.NaiveBayes --|ERROR: no class attribute set in test data!      and below is a screenshot of the KnowledgeFlow: ![setup1](http://i.stack.imgur.com/wYnpR.png) **SETUP 2: Use ClassAssigner before FilteredClassifier** Same as SETUP1, but the Class is assigned for the training and test set separately before using the FilteredClassifier (screenshot below). Then I don't get any error, the training part seems to work fine, but nothing happens in the test set, it is just loaded. For example if I try to visualize it after assigning the class (workflow below), I don't get any results to see. ![enter image description here](http://i.stack.imgur.com/FDluX.png) Here is the log:                   22:33:31: [Loader]  CSVLoader$13145810|-M ? -E ",'|loaded   finename_train             22:33:31: [Classifier] FilteredClassifier$13512466|-F "weka.filters.MultiFilter -F \"weka.filters.unsupervised.attribute.NumericToNominal -R last\" -F \"weka.filters.unsupervised.attribute.Remove -R first\" -F \"weka.filters.unsupervised.attribute.Remove -V -R numbers -F \"weka.filters.unsupervised.attribute.Normalize -S 1.0 -T 0.0\"" -W         weka.classifiers.bayes.NaiveBayes --|starting executor pool ( 2slots)...          22:33:31: [Classifier]  [Same as above] --|setup output queues.           22:33:31: [Classifier]  [Same as above] --|scheduling run  1fold   1for execution...         22:33:31: [Classifier] [Same as above] --|storing model for run 1fold 1         22:33:40: [Loader]  CSVLoader$22945919|-M ? -E ",'|loaded   test_FNC_SBM_1      I am at a complete loss and have been trying to make this work for several hours. I checked the test .csv and the column of NaNs is there. I visualized the test set in Weka by using DataVisualizer directly after TestSetMaker and the Class seems to be read as missing values: if i set the data points to be colored by the class they are shown as gray "M"s which i guess stands for "missing". If i set the points to be colored by another attribute they are shown as blue-to-orange crosses. So i suppose the data is here, but the workflow is incorrect. Any idea ? Thank you. It seems there is a similar question on StackExchange: http://stackoverflow.com/questions/8112041/error-while-using-weka-api-in-java- code-class-attribute-not-set but this is not using the GUI and it seems to suggest a setup similar to setup 2 which did not work in my case.


title: Is there a name for this process/ distribution?

text: Does the equation below have a name, or is it similar to some other well-known process/ equation? Equation of interest: $$S_c = S_{c-1} + S_{c-1}\omega_c\delta_c$$ $\delta\sim\mathcal{N}(0,1)$ is a standard normal shock $\omega\sim\mathcal{B}(1,~0.5)$ is a coin flip with a 50% chance of being 1 (0 otherwise) $c=(1,~2,~\dots,~C^{\ast})$ is sequence of coin flips ($C^{\ast}$ total flips) $S_0\sim\mathcal{N}(0,1)$ initializes the process $S_{C^{\ast}}$ is the value of $S_c$ at the end of the iterations, and is the value of interest. Note that the $c^{th}$ coin flip dictates whether or not the $c^{th}$ standard normal shock occurs, and thus whether or not $S_c$ is different from $S_{c-1}$. Below is some R code that will generate one value of $S_{C^{\ast}}$ when $C^{\ast}$ is 4. It's set up to make the omegas and deltas the same each time, but to allow the order in which they accumulate to be random. Note that the order of accumulation for a given set of omegas and deltas does not affect the value of $S_{C^{\ast}}$.               set.seed(2) # set the RNG seed so that all the shocks and coin flips are the same          Cstar <- 4 # number of coin flips/ number of iterations     S0 <- rnorm(1) # initial standard normal shock     Sc <- c(S0, rep(NA, Cstar)) # vector to hold the shocks as they accumulate          deltas <- rnorm(Cstar) # standard normal shocks     omegas <- rbinom(Cstar, 1, 0.5) # coin tosses          rm(.Random.seed, envir=globalenv()) # remove the non-random seed, so the order in which the shocks/ coin flips are accumulated is random     c.2.Cstar <- sample(1:Cstar, Cstar) # a scrambled sequence from c=1 to C*     ctr <- 1 # the counter used for indexing Sc (but not the omegas and deltas)          for(i in c.2.Cstar){         ctr <- ctr+1         Sc[ctr] <- Sc[ctr-1] + Sc[ctr-1]*omegas[i]*deltas[i]      }          SCstar <- Sc[ctr]     print(c.2.Cstar)     print(SCstar)


title: Little help with prediction

text: I have energy consumption data from a device for the past two weeks. I have to predict energy consumption of this device for another week. I am not sure how to forecast this. Does anybody know how to do that? If there a way to do that in R? I've just started learn the R language.


title: Creating index from three variables in panel data

text: I have three variables which I want to combine in one. The obvious choice is factor analysis or principal component analysis but it's not so simple with panel data because it needs means and standard deviations and in absence of information about panel structure, software calculates the mean and sd for the whole pool. So here is the list of my variables (all continuous variables): 1) human capital variable index based on average years of schooling and returns to education 2) proxies for quality of human capital 2(a) journal publications per worker 2(b) patent applications per worker I'd like to combine these three into one index which will be called quality adjusted human capital. I can't find a way to do it other than taking simple average of ratios in 2(a) and 2(b) and multiply it with (1) but this doesn't sound right. Any ideas, please?


title: Maxent Markov models in R

text: Is there a package that implements Maxent Markov Models in R? (http://en.wikipedia.org/wiki/Maximum-entropy_Markov_model). I understand that package crf implemets conditional random fields which are related, and there are packages implementing hidden markov models. But I am specifically looking for packages that implement maxent markov models.


title: What is the most appropriate method to recover words transcription from the phoneme sequence with errors?

text: I conduct some experiments on continuous speech recognition. After the initial application of the recognizer I have the phoneme sequence that includes some errors (three types of errors: substitution, deletion and inclusion). For instance: I obtain _h e l' l o n y k f r a n d s t_ instead of _h e l l o m y f r i e n d s_ So, the task is use the information about language (correct transcription of every word, valid word order etc.) in order to recover the correct transcription from the phoneme sequence corrupted by noise. I wondered about Hidden Markov Models, but I am not shore how to use it here. Any idea is appreciated. Thanks in advance!


title: Does PCA have any advantages or usages in the frequency domain?

text: My question is about analysis of signals with PCA in the frequency domain.   * As frequency analysis offers a powerful tool for signal processing, does Principal Component Analysis (PCA) have any "well-known" advantages or usages in the frequency domain?    * Can it be stated that PCA is doing a kind of "frequency analysis" as the principal components represent the most important frequencies in the signal?


title: Are these descriptions of batch gradient descent algorithm conflicting each other?

text: The first one is from Andrew Ng ![enter image description here](http://i.stack.imgur.com/JNfEv.png) The second one is from Francis Bach ![enter image description here](http://i.stack.imgur.com/hvFdv.png) I might be a little confused, but why is there a summation of partial derivatives in the second description and none in the first one? Are they talking about the same thing after all?


title: How to interpret the conditional expected value in a regression?

text: If you want to estimate $E[y| x = \text{some value}]$ is this just a matter of plugging $x$ into the regression equation? Because you estimate the regression coefficients, therefore the value you get is an estimate also?


title: How are epidemic models simulated in case of mobility?

text: I don't know whether this is a pure statistics, math or programming question so kindly let me know if there is a better place to post this question. I am trying to implement the SIS epidemic model when the nodes have mobility. I understand how to perform this simulation in an analytical fashion. However, things get rather confusing when nodes are mobile. The model assumes that each node can infect **any** node and hence the equations are valid. But when nodes are mobile, each node is **not able** to infect every other node (the other node might not be within the range) and has to explicity send a message to a node that is susceptible in order to infect it. In that case, give an infection rate `B`, how do I simulate this when the nodes are mobile? Currently, the way I am doing this is in the following way:               def Controller():         for i in range(1,100):             randNum = getRand()             if (randNum <= InfectionRate):                 neighbors = getNeighbors(i)                 ScheduleTransmission(getCurrentTime(), i, neighbors)         Schedule(getCurrentTime() + 1, Controller)      My problem is that I am not understanding if the infection rate can now be captured through a single value (which was previously `B`). If not, how does one analyze this scenario? Do I set the `InfectionRate` as `B/numNodes` so that the overall probability will be `B`? Any suggestions? **UPDATE: Back calculating Beta and improved infection strategy**               def Controller():         for i in range(1,100):             neighbors = getNeighbors(i)             k = len(neighbors)             for j in neighbors:                beta = -k log(1-c)                if (beta <= InfectionThreshold):                      ScheduleTransmission(getCurrentTime(), i, j)         Schedule(getCurrentTime() + 1, Controller)


title: Include specific interaction terms

text: I have a question regarding interactions in GLM. I run a Poisson regression with the purpose of predicting claims in insurance. I have the following problem :               model 1 : CLAIMS ~ X1 + X2     model 2 : CLAIMS ~ X1 + X2 + X1:X2      X1 has 10 levels, X2 has 4 levels When I do a deviance test (Chi Square), my result is that model 2 is significantly better compared to model 1 (at the 5% level). On the other side the AIC of model 2 is worse. (higher) When I look at the individual parameter estimates and their corresponding standard error, I see that only 3 out of 27 levels from X1:X2 differ significantly from the base level. I am not sure how to interpret these observations. Can someone help me and say something meaningful about it? And more important, how to deal with it? My attempt would be: "There are interactions in the model that are meaningful, however most of them are not. This explains the increasing AIC, although the deviance test says that model 2 is the better one." Should I try to construct a few dummy variables for these three interactions that seem significant? (However, I am not sure this is a correct approach, because it is only significant relative to the base level...)


title: Which observation has the largest variance of the residual?

text: > a) For which of the following observations (obs1, obs2, obs3) is the > variance of the residual the largest? Which observation has the highest > leverage? And which one the smallest? Explain why. > > b) Which of the three observations has the largest DFBETAS? > > c) About which observation would you worry the most? > > ![plot](http://i.stack.imgur.com/uk9lO.jpg) I thought that the answer to a) would be obs 1, because I think that will have the highest leverage because the distance from the cloud is the greatest (not quite sure whether that is true though). The rest of the questions I had no clue about. Could anyone please help?


title: Does including both raw  and per capita measures as predictors reduce significance of either predictor?

text: I'm running a regression on independent variables, some of which are measured in different units, for example:   * The amount of broadband connections in a country   * The amount of broadband connections in a country per 100 people Since the second is not perfectly correlated with the first (it is divided by the population divided by 100), I assume there is no issue including both measurements in an OLS regression. The question is, will it diminish the significance of either coefficient?


title: Liblinear types of solver

text: There is many variants of type of solver in liblinear but I don't understand their differences.Which one I must choose? Also why data must be scaled? duo to some numerical issues? > >     -s type : set type of solver (default 1) >       for multi-class classification >        0 -- L2-regularized logistic regression (primal) >        1 -- L2-regularized L2-loss support vector classification (dual) >        2 -- L2-regularized L2-loss support vector classification (primal) >        3 -- L2-regularized L1-loss support vector classification (dual) >        4 -- support vector classification by Crammer and Singer >        5 -- L1-regularized L2-loss support vector classification >        6 -- L1-regularized logistic regression >        7 -- L2-regularized logistic regression (dual) >


title: Is it inappropriate to do a One-way independent ANOVA when one of the levels has only one participant?

text: I sense this is not a good idea, but am struggling to think of an reason why it is wrong, other than that the assumption of homogeneity of variances would surely be violated. Is there some more fundamental reason why this is wrong? Does it depend on how many participants are in the other levels of the factor, or on how many other levels there are?


title: Island Hopping with Metropolis Algorithm

text: From John Kruschke's book, Chapter 7, Pg. 120 (summarised for succinctness): > A politician is constantly travelling from island to island on a chain of > islands... His goal is to visit all the islands proportionally to their > relative population, so that he spends proportionally more time on the more > populated islands...He doesn't even know exactly how many islands there > are!.. His advisers can ask the mayor of the island they are on how many > people are on the island. And when the politician proposes to visit an > adjacent island (to the east or west, by flipping a coin), they can ask the > mayor of that adjacent island how many people there are on that island. I wrote some Python code, where `A` is the island chain and `run_Metropolis(A, N, start, burn_in)` does the computations. The code is intentionally dumbed- down to make it more human-readable.               import random     def run_Metropolis(A, N, start, burn_in):              ''' A: the distribution to be predicted (a list).             N: the number of trials / moves             start: starting index in A             burn_in: the number of  initial moves to be discarded.         '''              Posterior = [0] * len(A) # to be returned         cur_pos = start # starting index              for i in range(N):                  # proposal distribution is (-1, +1) cordinates, takng into account the list boundaries             if cur_pos<len(A)-1 and cur_pos>0:                 proposed_positions = [cur_pos-1, cur_pos+1]             elif cur_pos == len(A)-1:                 #proposed_positions =[cur_pos-1] # if this is enforced, the posterior does not look like A                 proposed_positions = [0,cur_pos-1]              elif cur_pos == 0 :                 #proposed_positions = [cur_pos+1] # if this is enforced, the posterior does not look like A                 proposed_positions = [len(A)-1, cur_pos+1]                  proposed_pos = random.choice(proposed_positions)                  # decide whether you will move             if A[proposed_pos] > A[cur_pos]: # definitely move                 cur_pos = proposed_pos             else: # move with a prob proportional to how close A[cur_pos] is to A[proposed_pos]                       u = random.random() # sample from uniform distribution from 0 to 1                      if float(A[proposed_pos])/A[cur_pos] > u: # move to proposed                      cur_pos = proposed_pos                  else: # stay at cur_pos                     pass                   # if past brun_in point, increment the location by 1 in the posterior distribution.             if i > burn_in:                  Posterior[cur_pos] +=  1              return Posterior           A = [20,30,10,50,80,20,90] # populations of the islands      posterior_A = run_Metropolis(A, 4000, 3, 1000)      # Validate the sampling works:      [n/float(sum(A)) for n in A]      > [0.0666, 0.1, 0.0333, 0.166, 0.266, 0.066, 0.3]       [n/float(sum(posterior_A)) for n in posterior_A]      > [0.0593, 0.0910, 0.0393, 0.177, 0.263, 0.0666, 0.302]          So, `posterior_A` does reflect the population distribution of `A`. What Kruschke does not make clear is that **in order for this to work, the island chain needs to be cyclic** , i.e. you need to be able to travel from the 7th island to the 1st one and vice versa. Otherwise, `posterior_A` does not look like `A`. You can test this by changnig the boundary proposal distributions with the commented out ones, which don't assume cyclicity. **Questions:** 1.Is this because the ergodicity assumption is being violated? 2.Can we modify the acceptance criteria (move decision) to account for the lack of symmetry for the easter-most and western-most islands? (I guess this is Metropolis-Hastings) 3.If so, can you please modify the code to illustrate?


title: How to visualize this data?

text: ![enter image description here](http://i.stack.imgur.com/sX1yH.png) Here I have shown the first two lines of my data frame. It goes on for a number of lines. I am wanting to make plots to show what is happening over time . e.g a plot showing the sequence of grip changes over time and who is making the changes P1 or P2 . I am using R, so things I can implement in R would be relevant - and guidance on how to do so would be especially useful. I have tried using code such as `plot(mydata$Start,mydata$Name)` and `stripchart(mydata$Name)` but I have not been successful. In some cases I get an error such as               Error in stripchart.default(mydata$Name) : invalid first argument      I am unsure how to approach this. Should it not be a data frame maybe ?


title: How to apply multiple testing correction for gene list overlap using R

text: I have 2 studies looking at the patient response to the same drug. Study 1 found 10,000 genes expressed above the background and 500 of them are differentially expressed and referred to as the drug response signature. Study 2 found 1,000 genes representing the drug response signature. The overlap between the two signatures is 100 genes. I want to calculate the statistical significance of the overlap between the signatures. If I understand correctly, one way of doing it (based on the posts here: Calculating the probability of gene list overlap between an RNA seq and a ChLP-chip dataset and here: Using R's phyper to get the probability of list overlap) is via `phyper()`:               > overlap  <- 100     > list1    <- 500     > totalPop <- 10000     > list2    <- 1000     >      > 1-phyper(overlap-1, list1, totalPop-list1, list2)     [1] 4.103051e-12        1. Does that sound reasonable?   2. If I wanted to apply the Bonferroni correction, I would need to multiply this p-value by the number of comparisons. What would the number of comparisons correspond to in this case? List2? Alternatively, what would be a quick way of doing less conservative correction (e.g., Benjamini-Hochberg)?


title: How to use Ridge Regression for classification? (Or other suggestions)

text: I have a sparse regression problem (Sparse because a few inputs are factors so we have a lot of columns of 1s and 0s). I am thinking of Ridge Regression because of the sparsity, but also because a lot of the terms will have interaction effects. I also want an interpretable model. Is there a way to use the Ridge Penalty for a Linear Regression Classifier? If not, is there any base learner which allows me to include interaction effects and still result in a sparse solution.


title: Smoothing algorithm for saturating function

text: I have a noisy readout of a curve that is monotonically increasing or decreasing for a narrow range of points and then quickly saturates. I don't know exactly where the saturation point is, but from the underlying (biological) process I know that it exists and once it's reached, strong deviations from a flat line are very likely to be outliers. I am looking for the best way to smooth these data. Something like this (the data is actually not for time series, but it's a good analogy): ![Example of data and a curve I'd like to fit](http://i.stack.imgur.com/MhPhA.png) So far I have tried using:   * the Hampel filter   * loess regression (including its robust version, family="symmetric" in R)   * a combination of both The problem with loess is that since the non-saturated bit is so small, loess tends to "underfit" it, "dragging" the trendline up or down towards the saturation point. The problem with the Hampel filter is that while it does a great job in filtering out strong outliers, it's not that great for smoothing as such, i.e., polishing out small deviations from the trendline. It seems like the best thing would be to find the saturation point and then fit two separate models "before" and "after" it. This would be straightforward without the noise, but I'm not sure how to do it in a noisy setting. I'm also open to any other suggestions. This is part of a large project in R, so if a suggested algorithm has an R implementation this would be particularly great.


title: Number of areas in conditional autoregressive models

text: This is a simple question on Bayesian spatial modelling via conditional autoregressive modelling. What is, according to your judgement (and possibly some methodological source), the minimum number of regions for which CAR is suitable? Would you say that with a low number of regions, say 9, it is highly advisable to use a model with as little parameters as possible, that is, the intrinsic model instead of the better performing Leroux specification? Reference: Lee 2011: A Comparison of Conditional Autoregressive Models Used in Bayesian Disease Mapping, Spatial and Spatio-temporal Epidemiology 2/2, p 79–89


title: Comparing multiple datasets from a controlled experiment

text: I have a system to which I make certain changes. I want to compare the impact of those changes on the system performance. For this, I record the performance as a function of a set of inputs multiple times, before and after the changes have been made. How do I quantify the impact of these changes given that random environmental effects apply on the system at all times? I believe this is a generic problem in statistics, and have read ANOVA, Student's t-test and K-S test. However, what considerations have to be made while deciding on a method? There is no assumption of normal distribution, so it seems that K-S test is the best solution. Can you point me to a resource that discusses this class of problems? I am pretty confused by the fragmented explanations I have found online so far.


title: How should I interpret the result of my OLS (involving H-statistic)?

text: Please consider the following excerpt from "Measuring and Explaining Competition in the Financial Sector" by Bikker and Spierdijk (2008): > "In many economic theories, competition is related to the (relative) size of > a mark-up on the cost price as a component of the output price. However, > data on the price-cost margin (PCM) are generally not available in the > financial markets. Therefore, we have to measure competition indirectly. > [...] This paper gives a survey of measures of competition of over 100 > countries based on the **Panzar-Rosse (P-R) model**. [...] This approach > measures how total interest revenues of banks in a country or market react > to changes in input prices. A firm’s competitive behaviour in the market is > reflected by the degree by which input price changes are passed through to > output prices and to changes in output volume. The P-R model produces a > certain **H-statistic** which under certain conditions reflects the degree > of competition with H=1 pointing to perfect competition and H≤0 indicating > monopoly or a perfect cartel. The range where H is 0 to 1 denotes > monopolistic competition or oligopoly of some sort." I have used this H-statistic in my (OLS) regression model where companies' social behavior is the dependent variable (on a scale from 0 to 100). There are about 2,500 companies from 30 countries in the dataset (each country having an `H-statistic`). My hypothesis was that competition has a curvelinear effect on social behavior:   * too little competition: low social scores because consumers have no alternative anyway so companies have no incentive to act socially responsible;   * intermediate competition: high social scores because social responsible behavior becomes a unique selling point;   * too much competition: low social scores because companies cannot be bothered with social programmes since the intense competition drives them to save money where they can in order to turn profit. To see whether this hypothesis is right, I first included the `H-statistic` in the OLS and second the `(H-statistic)^2`. I would have expected a positive coefficient in the first case, and a negative one in the latter. However:               H-statistic       0.150 (p<0.01)     (H-statistic)^2   0.159 (p<0.01)  >>> instead of -0.159 what I would have expected      My question is: **how should I interpret this result?** In other words: based on this outcome, how would you describe the relationship between competition in the banking industry and the social scores of those banks?


title: Implementing Kalman filter in R

text: I am attempting to implement a Kalman filter in R. Working from "Time Series Analysis" by Hamilton the state space equations are defined as: State: $z(t) = F_z(t-1) + v(t)$ Observation: $y(t) = A'x(t) + H'z(t-1) + w(t)$ Where $x(t)$ are exogenous variables. The model I am attempting to implement is defined this way. There are multiple packages that implement state space models in R, but through reading the documentation it appears that the equations are defined differently, with $A'x(t)$ including in the state equation. Short of writing the maximum likelihood estimation code myself, does anyone know a way to fit this model in R? Update: Deleted my past update. I found another paper that implements this in a way that is much more clear, will update after I have digested it. Update 2: So I found another paper that does something similar to what I am attempting, but it is much more clear to me how to set up the state space equations. The problem I am running into now is running the filter using package dlm. If anyone has any experience that would help I would be greatful. I am able to set up the model, and then run the filter, but only 2 filtered values are returned (initial state and t=1) as opposed to t-values. I am not really clear why that is. From the dlm documentation: Observation: $y= \mathbf{F}\theta + v_t$ State: $\theta = \mathbf{G}\theta_{t-1} +w_t$ Where m is the number of observations and p is the number of parameters in $\mathbf{F}$   * $\mathbf{V}$ is m x m   * $\mathbf{G}$ is p x p   * $\mathbf{W}$ is p x p   * $c_0$ is p x p  Sample code:               require(dlm)          FMat<-cbind(rnorm(123),rnorm(123))       GMat<-diag(1,2)     WMat<-diag(0.02,2)     Vmat<-diag(0.02,123)          m0Vec<-c(1,1)     c0Vec<-diag(0.05,2)     y<-rnorm(123)          ssMod<-dlm(FF=FMat,V=Vmat,GG=GMat,W=WMat,m0=m0Vec,C0=c0Vec)     kFilter<-dlmFilter(y,ssMod)          kFilter$m      > kFilter$m [,1] [,2] [1,] 1.00000000 1.00000000 [2,] -0.02417179 0.05107227


title: Models for learning acquistion process

text: Imagine such a memory test for a mice. The mice performs an experiment $E$ with two possible issues $0$ (failure) and $1$ (success). If the mice gets $1$ it is "rewarded", if it gets $0$ it is "punished". The experiment $E$ is repeated many times with the same mice, and we consider the stochastic process $(X_n)$ modelling the successive experiments. I wonder about the way to analyze such data. The mice is supposed to learn over time so we cannot assume independence and we cannot even assume a homogeneous Markovian transition. We can assume inhomogenous Markov for simplicity. We are interested in the transition probabilities $p^n_{ij}=\Pr(X_{n+1}=j\mid X_n=i)$ ($i,j \in\\{0,1\\}$). Actually the two probabilities $p^n_{01}$ and $p^n_{11}$ determine the Markov step at time $n$. With several mices it should be possible to estimate the $p^n_{ij}$. But with a moderate sample size we could assume a parametric constraint to have "power". Even with a large sample size it is interesting to assume a parametric constraint in order to see a "trend" over time. For instance we could assume a logit relation $\text{logit } p^n_{ij} = \alpha + \beta n$. At first glance it should be easy to estimate this models with maximum- likelihood. My requests are:   * Are there some statistical models of the same spirit in the literature ? I'm interested to get some references (not necessarily restricted to the Markov assumption).   * I am also interested in getting some references --- if there are some --- about the more complex case of a $\\{0,1\\}$-valued continuous-time process $(X_t)$ modeling in addition the duration of the experiment $E$. More precisely, the mice can achieve success in a more a less rapid time, and the experiment is stopped after one minute if the mice has not yet succeed. For instance, we could assume for simplicity that the mouse always succed experiment $E$ after a time following an exponential distribution.


title: Detect outliers in very small data set

text: I have a data set that includes the different response times of a user that is visiting a web application. For example, a visitor enters www.test.com in the browser and navigates through this domain watching child pages like www.test.com/news, www.test.com/overview, www.test.com/overview/current, etc. If a user watches a web site, it is called a user action. Let's say a user has performed 5 user actions with the response times 200ms, 500ms, 350ms, 1200ms, 154ms. Now I want to find the outliers that express either fast page loads or slow page loads. Is that somehow possible? Thanks EDIT: I want to detect outliers because I want to determine the user experience depending on the response time. Let's say I have three ux-states, namely happy, ok and unhappy. All user actions are ok except the outliers. They are either unhappy if the response time is too high or happy if the response time is very low.


title: Interpret scale parameter from GEEGLM output in R

text: In GEE, the variance of each outcome is equal to the variance function multiplied by the scale parameter. For a continuous outcome the variance function =1. I understand that the scale parameter for a continuous outcome is always estimated in a GEE analysis, as opposed to set at 1 as for a poisson or binary outcome. Therefore this means that the variance is always equal to the scale parameter (when estimating a continuous outcome with GEE). Is this correct? If so, does this mean that a higher scale parameter is better than a lower one? e.g., a scale parameter of 0.9 means that 90% of the variance is explained by the GEE? Thank you in advance!


title: Identifying outliers for within-patient numerical data

text: I have a simple dataset of people's heights, many people with measurements on multiple days (once a year for 10 years, say). I have the date of each measurement. Some of the height values are absurd. I already drop values that are 'impossible' (e.g. height values above 3m). However, I'd like to identify unusual values within a patient. If a patient has 5 records at around 1.8m, I'd like a red flag if there's a value for 2.1m I am considering trying the ESD test for outliers - it seems straightforward to implement - but I thought I'd ask if there are better ideas before I get started. Thanks. EDIT: The ESD test I refer to is the [[Generalized] 'Extreme Studentized Deviate'][1] test; similar to the Grubbs Test, it 'is used to detect one or more outliers in a univariate data set that follows an approximately normal distribution'. I'll have measurements for over 100,000 people - however, each person will only have between 0 and around 20 measurements. I'm hoping to generalise this to other measurements such as BMI, lab tests (Eosinophils, White Blood Counts), respiratory tests (FEV1), etc. I suspect the latter ones will be difficult as they're more varied, rising and falling with little dependency on time - hence my starting with height!


title: Imputation and Distributions

text: Suppose you impute a variable using a normal distribution with mean 10 and sd 5. Is it better to draw 1000 random samples from this normal distribution, take the average, and then use this to impute this missing value? Or can I just take one draw from this distribution and use this as the value to impute the missing data?


title: Propensity score matching: using alternative methods to create a distance measure

text: I would like to use a greedy nearest neighbour method to do propensity score matching. Though I've little experience here, it seems that the distance measure used is generally a propensity score generated from a logistic regression. My question is: why logistic regression? Why not a random forest, SVM or another method? Is there some logic to suggest this would be unfruitful? My plan is currently to use the `MatchIt` package in R and input my own distance measure calculated off the back of a random forest (you can input your own propensity score into the argument `distance` of the `matchit` function). I'm deterred by the fact that modern propensity scoring packages such as PSAboot don't have built in facility to do this _for nearest neighbour methods_. They do use `party` and `rpart` but only to match using _strata_ (additionally they only use single trees). I'm intrigued that methods with (typically) greater predictive power than logistic regression do not appear to be harnessed to create a 'better' propensity score for one-to-one matching. Is there someone out there that can shed some light on this? Linked question: Propensity Score


title: Predicting team performance

text: Are there any Machine Learning/ Datamining research paper which deals with predicting the performance of two team containing each $x$ players?


title: Pearson's $\chi^{2}$ value?

text: Pearson's $\chi^{2}$ value is 115.778 and 4 cells (33.33%) have expected count has less than 5. I am confused how to interpret this output. I have a 2 by 6 contingency table. Can someone tell me what has gone wrong?


title: What are problems of many hidden layers?

text: Is there any problem if we use too many hidden layers in Neural Network? Can anyone simply describe what problems can occur if we have too many hidden layers.


title: Bounded multiple linear regression

text: I am attempting to build a model a for which the only valid output in the range [0,100]. I was wondering if would be possible to reduce the penalty on values under 0 and over 100 as they will be constrained anyway in order to get a better fit.


title: How to define the rejection region for a permutation test?

text: I am trying to find out the rejection region for a permutation test about the difference between medians where:                H0: Difference between medians equals 0;       Ha: Difference between medians is greater than 0.       Considering that I have the distribution of my test statistics on resampled datasets, how can I find the critical value for the test so that I can infer the rejection region?               > sum(permstats>=obs.median)/length(permstats) #pvalue      It is not completely clear to me how can I find the rejection region according to the permutations distribution( and maybe the p-value??)


title: Significant logistic regression model but poor goodness of fit

text: What should one do, if the parameters for a logistic regression model are significant as well as the model itself, but the test for the goodness of fit show that the model is bad (i.e. the P-values of deviance etc. are low)?


title: How is  the SAS REPLICATE statement in PROC NLMIXED different from a FREQ statement?

text: Note: I do not have access to `SAS` to work this out for myself. I'm trying to work out the difference between the `FREQ` statement in `PROC GENMOD` (for example) and the `REPLICATE` statement in `PROC NLMIXED`. The `SAS` explanation is not illuminating the practical mathematical difference for me. The explanatory test for `REPLICATE` is (Documentation link here): > The REPLICATE statement provides a way to accommodate models in which > different subjects have identical data. This occurs most commonly when the > dependent variable is binary. When you specify a REPLICATE variable, PROC > NLMIXED assumes that its value indicates the number of subjects having data > identical to those for the current value of the SUBJECT= variable (specified > in the RANDOM statement). Only the last observation of the REPLICATE > variable for each subject is used, and the replicate variable must have only > positive integer values. Note that the REPLICATE mechanism is different from > using a FREQ statement in other statistical modeling procedures, such as > PROC GLM, GENMOD, GLIMMIX, and LOGISTIC. A FREQ variable is used to identify > grouped values for observations, essentially multiplying the log likelihood > or sum of squares contribution for the observation. A REPLICATE variable is > used to multiply the contribution of a subject that comprises one or more > observations. The explanatory text for `FREQ` in `PROC GENMOD` is (Documentation link here): > The variable in the FREQ statement identifies a variable in the input data > set containing the frequency of occurrence of each observation. PROC GENMOD > treats each observation as if it appears _n_ times, where _n_ is the value > of the FREQ variable for the observation. If it is not an integer, the > frequency value is truncated to an integer. If it is less than 1 or missing, > the observation is not used. In the case of models fit with generalized > estimating equations (GEEs), the frequencies apply to the subject/cluster > and therefore must be the same for all observations within each subject. Assuming the replicate variable contains only integer values, that sounds like the `REPLICATE` statement in `PROC NLMIXED` is just a `FREQ` statement. But the `SAS` documentation says otherwise. What have I missed with how the commands affect the number of observations used, and the flow on effects to estimates? In particular, how does this affect repeated measures studies? Update: could the answer lie with differences in treatment of the degrees of freedom? I don't understand how `SAS` uses these commands to figure this out for myself.


title: Analysis of infection data between two years using parametric or nonparametric tests?

text: I have three years of infection data from 6 different hospitals.                     2010     2011    2012     1      1.3       2.5    1.03     2     3.35       2.5    2.6     3     4.1        1.8    2.5     4     2.0        3.65    .7     5     3.4        3.75    2.4     6     1.2        3.4     1.3      These numbers are rates (infections/10,000ptd). I need to determine significance between years. What I plan to do is run either an unpaired t-test for 2010/2011 and 2011/2012 or a Mann-Whitney U test. Since I only have three years of data, I feel that to compare all the data with ANOVA there wouldn't be high enough power to find significance. Does this sound reasonable or are there alternate testing I can do to find out if there is statistical significance between the years?


title: Is expectation of univariate variable dependent on the multivariate

text: Suppose I have a mutivariate Gaussian distributed variable $u\sim\mathcal{N}(\mu,\Sigma)$, where $\Sigma$ is a dense matrix. I wish to calculate the expectation of $f(u_i)$. Is $$E(f(u_i))=\int f(u_i)\mathcal{N}(u_i|\mu_i,\Sigma_{ii})\,du_i$$ If so provide a simple proof? Also does this apply to other multivariate distributions? eg. Multivariate T. Otherwise do I have to consider that $$E(f(u_i))=\int f(u_i)P(u_i|u_{\backslash i})\,du_i \,P(u_{\backslash i})\,du_{\backslash i}$$ where, $u_{\backslash i}$ indicates variable $u$ without the i-th dimension. Aside: in my particular case $f(u_i)=\exp(-u_i)$, however I'm more concerned of the general case.


title: How does one fit conditional Poisson regression in SAS?

text: I'm interested in fitting a conditional Poisson regression model using PROC GENMOD in SAS to analyze a matched cohort study. However, it's not quite clear to me how I should exactly go about it. My impression is that a REPEATED statement should be used along with the events/trials syntax, but if so then how does one account for continuous covariates in the model, such as age? I would be most grateful if anyone had any pointers. Thus far, Google and the internet have disappointed me. Cheers.


title: Count data forecasting/prediction

text: I would like to know if the normal forecasting methods apply for count data, in specific a dataset that contains several zeros? I have data set that counts the usage of a service on an hourly basis, and I would like to make a forecast on the hour in which the service is going to be used based on the historical data. Also, I would like to know if the correct term for what I want to do is forecast or prediction?


title: Figuring out an effect size with a P value for a t-test, sample means, and sample sizes

text: An author did not report an effect size so I am hoping you might be able to help me calculate it. I know the Means are:               Condition 1: 326     Condition 2: 558      Sample size:               Cond1: 11     Cond2: 12                    T-test p value: .04


title: Modifying stochastic meta-descent for variable weight

text: In [1], Schraudolph presents an algorithm for updating a weight $w$ and some auxiliary variables $p$ and $v$ given a vector gradient $g$ and another vector "$Cv$". The update rules attempt to make it so that v goes to $C^{-1}g$ where $C$ could be, for example, the Hessian. $Cv$ is a vector like $g$ that is passed into the algorithm. The update rules are: \begin{align}\DeclareMathOperator{\diag}{diag} \DeclareMathOperator{\max}{max} w_{t+1} &= w_t - \diag(p_t)g \\\ p_t &= \diag(p_{t-1})\max\left(\frac12, 1 + \mu\diag(v)g\right) \\\ v_{t+1} &= \lambda v_t + \diag(p_t) (g - \lambda Cv_t) \end{align} where $\mu$ and $\lambda$ are constants. How would you modify these update rules so that the application of training can be weighted by $0 \le u$ such that it's similar to training $u$ times with the unmodified algorithm. My first attempt was to replace $g$ with $ug$ and replace $\lambda$ with $\lambda^u$. I'm looking for insight. [1]: Schraudolph, N. N. (2002). Fast curvature matrix-vector products for second-order gradient descent. Neural computation, 14(7), 1723–38.


title: A question about paired sample t test?

text: Your help on the following question is appreciated. Thanks. I know paired-sample t test is for the same group of people to test their pre and post tests and see if they have a significant gain between their pre and post tests. If the treatment is progressive, can I use the same people to collect data over time and run paired-sample t test together because the sample size is very small (only 6 persons)? Year 1: Treatment 1 with 3 persons' (ABC) pre and post tests Year 2: Treatment 2 (improvement of treatment 1) with the same 3 person's (ABC) pre and post tests Year 3: Treatment 2 (improvement of treatment 1) with the other 3 persons' (DEF) pre ad post tests In this case, can I run paired-sample t test for 9 data points (3 data points [ABC] from Year 1, 3 data points [ABC] from Year 2, and 3 data points [DEF] from Year 3) and then report the effect of treatment 1 in general (since treatment 2 is based on the improvement of treatment 1)? I know some people might argue that I use the same people twice (ABC) in year 1 and year 2. But can I argue that paired-sample t test is concerned with the "gain" between pre and post tests. Thus, although they might have a higher score for the pre test in the second year, we expect them to have a higher score for the post test too. That is, we expect ABC to have a "gain" in both Year 1 and Year 2. Your insights and advices are much appreciated. References are much appreciated too.


title: Estimating growth over sets of clustered data

text: The title may not be entirely accurate for what I'm looking for. I am measuring the growth of a database table over time in order to estimate how much disk I will need in the future. With just a single table I guess I could apply linear regression to make my estimate, getting more accurate as more measurments are made. The issue, though, is that the database in question is multi-tenant. Customers range from very frequent users to very infrequent users, and each user has their own table. More frequent usage would equal greater growth over time. Given this, how can I make a reasonably accurate guess of the growth? In my understanding, the infrequent users would skew the data, and vice versa for frequent users. I dont see how I could cluster the users, since they might be considered a frequent user at measurement 1 but an infrequent user at measurement 2.


title: What if sample size is larger than population size?

text: I read something about standard error, which tells that sample mean is not accurate estimation because we do not sample full population of size N. But, what if sample size n = N or exceeds N, i.e. n > N? Can standard error can be more accurate than the standard deviation?


title: Fisher's exact test 2 x 4 contingency table possible to test one-sided?

text: Is it possible to test a one-tailed hypothesis with Fisher's exact test for more than a 2x2 contingency table? I am asking because the online calculators I found for this case only give a 2-sided value (e.g., like this one: http://vassarstats.net/fisher2x4.html). So would it nonetheless be legitimate to divide the given p-result by two and say I conducted a one-sided test?


title: Cronbach's Alpha and Factor Analysis for Repeated Measures Design

text: I have a repeated measures (time 1, time 2) experimental design and would like to run a confirmatory factor analysis on the scales I've used and to follow up with a cronbach's alpha. However I am unsure if I should run the analyses on time 1 or time 2 data. Any advise?


title: BMI at baseline & followup with exposure at baseline; model change or BMI at FUP? Control for BMI baseline?

text: For a prospective occupational cohort where everyone is exposed to one or more chemical agents, examining BMI at follow-up compared to a specific chemical exposure at baseline, is it necessary to control for baseline BMI? Is it better to model change in BMI or BMI at followup? There is no unexposed group -- just cohort members unexposed to some agents versus others. All analyses are within-cohort.


title: Entropy calculations and Process ranking

text: Still working on my financial datasets, I am now trying to observe the effects of the pre-processing on the results of my study. To sum it up: I have thirty time series of prices (the main indices, forex & bonds over ten years) and I want to observe markets trends. The main idea is to observe the correlation matrix of normalized log-returns and its highest eigenvalue. This should give the overall market trend. I obtained interesting results wich match the reality, and the dynamic of correlation we are supposed to have. I encountered a lot of problem with daily market prices. And there is many more solutions. For exemple:   * Market synchronization: due to different opening hours we could not really take the correlation calculated day by day. The solutions range from looking at lag one correlation, or at weekly moving average price, or at weekly price, to using GARCH estimation models.   * Data normalization: To use PCA data have to be normalized, but due to heteroskedasticity, it is tempting to use local normalization. Local normalization process use a lot of parameters, like time window and weight.   * Data filtering: We can use a moving average to remove HF noise, a kalman filter, or random matrix theory (see here: http://arxiv.org/pdf/0910.1205v1.pdf). All these methods use parameters, some are choosen intuitively to fit data properties, but others are only choosen with the wet finger technique.   * data window: I'm trying to look at the evolutions of correlation so I had to choose a moving windows and its width. The width is a compromise between the precision of the result (shape), and the number of time windows obtained. Not to forget about it, each time I had to think about a moving average or a windows, there is also a question of weight, for example the exponentially weighted moving average is commonly used in finance. Each time I had to choose a parameter or a method I had one question in mind: how to choose the better solution ? My technique was to use the correlation. I think real correlation is always underestimated in finance, so a transformation wich increase correlation is always a good transformation. To rank process I used the idea that an higher growth in correlation means a better transformation. I recently learned about entropy in statistics and its use in information theory. I now want to take a look at what it would bring to my study. I was thinking about using it's value to know if a pre-pocessing method keeps enough information to be used. So I have one or two questions: The entropy is defined by $ H(X)= - \sum_{i=1}^np(x_i)\log_b p(x_i)$. I don't have probability distribution, should I estimate it ? Or should I estimate entropy directly ? I am quite surpised that the dynamic is not taken into account. I mean that only the distribution of values is important. I have seen some study where the "shuffled" values are gaussian distributed but where the dynamic plays an important role. I think that returns are not independants from one day to another and we have to take it into account. Am I right ? So, how entropy is used to ranks methods in practice ? Thinking about entropy, I also wondered in wich way it could reveals something direclty of the market information itself. Have you heard of any study in this domain ? what was their principle ?


title: Simulating responses to a test for item-response theory

text: I am developing an online assessment system and I need to calibrate the bank of questions but I do not have enough people to implement a pilot test. That is why I decided to simulate the responses of the question bank. I will use the Rasch model for the development of a computer adaptive test (CAT); I have only questions from the bank and the number of options for each. How I can do a simulation using only this data? I had thought to make it generate a random number $X$, if $X < 1 / N$, where $N$ is the number of options in a particular question, then the person has gotten the answer correct, otherwise it is incorrect. **Is this a correct approach for calibrating the item bank to a Rasch model via simulation?** PS: What I want to do is to calculate the level of difficulty of each question. It is my understanding that to calibrate a question bank, that bank must be applied to a sample of subjects, but I do not have enough people to do this, so I decided to try a simulation approach. I tested the function `sim.rasch` in the `psych` package, and I see that I need to take into account the number of options for each question, as this may influence their difficulty.


title: Energy estimation through machine learning

text: Greedings to everybody. I have the dataset which you can find here, containing many different characteristics of different houses, including their types of heating, or the number of adults and children living in the house. In total there are about 500 records. I want to use an algorithm, that can be trained using the dataset above, in order to be able to predict the electricity consumption of a house that is not in the set. I have tried every possible machine learning algorithm (using weka) (linear regression, SVM etc) . However I had about 350 mean absolute error, which is not good. I tried to make my data to take values from 0 to 1, or to delete some characteristics. I did not managed to find some good results. I also tried to use R tool, and I did not have good results either... I would be very grateful, if someone could give me some advice, or if you could examine a little the dataset and run some algorithms on it. What type of preprocessing should I use, and what type of algorithm?


title: Rank correlation and demonstrating reasonable proficiency in rank ordering

text: **Background** I'm working on a new effort estimation model (dissertation work) which involves software engineers organizing tasks (rank ordering) in order of increasing effort based solely on their understanding of the requirements. For example, writing a function to _calculate an average_ should require less effort than writing a function to _calculate an integral_. A software engineer's ability to rank order tasks in a _reasonably proficient_ manner is essential to successfully constructing an estimate with my new model. Put another way, a software engineer should be able to get enough of the items in the correct rank order that the estimate becomes useful for project planning purposes. Some errors are expected and can be accounted for by the model (estimation is not an exact science, after all). I believe I can use one of several rank correlation measurements to compare the software engineer's rankings to the "correct" rankings (as defined by historical data of actual effort). **Question** How can I define _reasonably proficient_ in terms of Spearman's ρ, Kendall's τ, etc.? One possibility that has occurred to me is "better than random." Assume that anything better than the average of all possible arrangements indicates some level of intelligence in the ranking (ρ > 0 for Spearman, maybe?). This doesn't inspire a lot of confidence, though. Several references talk about the statistical significance of 10%, 5%, 1%, etc., but give no guidance on why these values are selected. I can't just make up a required level of significance. I do have a reference indicating that estimators should be be 90% confident in their estimates before using them for project planning purposes. Perhaps this can be applied in some manner. **Additional Information** The Pearson product-moment correlation coefficient article on Wikipedia describes an interpretation approach from a 1988 book (page 78) of Small, Medium, Large, and None based on the "size" of the correlation. Any thoughts or ideas?


title: Finding windows of high-correlation across coordinates in R

text: I've got a dataset that looks like this:               id      start   end     score1  score2     [...]         LmjF.31 280000  290000  .       .     LmjF.31 290000  300000  .       2686     LmjF.31 300000  310000  .       .     LmjF.31 310000  320000  .       74     LmjF.31 320000  330000  .       1897     LmjF.31 330000  340000  .       4116     LmjF.31 340000  350000  705     3156     LmjF.31 350000  360000  663     4145     LmjF.31 360000  370000  170     4864     [...]     LmjF.34 1140000 1150000 2284    2751     LmjF.34 1150000 1160000 .       .     LmjF.34 1160000 1170000 .       .     LmjF.34 1170000 1180000 572     .     [...]      (dot=`.` means score is 0) I would like to rank the regions for each `id` where the correlation between score1 and score2 is highest . The second and third columns are `start` and `end` positions for sliding windows that at their minimum can be of size 1, but with smaller windows, fewer of them have non-zero on both scores, so I am already using bigger windows. Since these windows in different chromosomes for two different species, correlations in `score1` and `score2` that don't overlap but that are close to each other would still be relevant, but if they partially or fully overlap, then even better. That's why the sliding window size is relative, since I would like to use a method that takes into account cases like the ones for LmjF.31 start=340000 to end=370000, where there is three consecutive rows with a score bigger than 0. Any suggestions?


title: Statistics, war stories, data intuition

text: I think it is fair to say statistics is an applied science so when averages and standard deviations are calculated it is because someone is looking to make some decisions based on those numbers. Part of being a good statistician then I would hope is being able to "sense" when the sample data can be trusted and when some statistical test is completely misrepresenting the true data we're interested in. Being a programmer that is interested in analysis of big data sets I'm relearning some statistics and probability theory but I can't shake this nagging feeling that all the books I've looked at are kind of like politicians that get up on stage and say a whole bunch of things and then append the following disclaimer at the end of their speech > Now, I'm not saying that this is good or bad but the numbers say it's good > so you should vote for me anyway. Maybe you get that but maybe you don't so here's a question. Where do I go to find war stories by statisticians where some decisions were based on some statistical information that later turned out to be completely wrong?


title: Meaning of partial correlation

text: From Wikipedia > Formally, the partial correlation between $X$ and $Y$ given a set of $n$ > controlling variables $Z = \\{Z_1, Z_2, …, Z_n\\}$, written $ρ_{XY·Z}$, is > the correlation between the residuals $RX$ and $RY$ resulting from the > linear regression of $X$ with $Z$ and of $Y$ with $Z$, respectively.   1. It says earlier that > partial correlation measures the degree of association between two random > variables, with the effect of a set of controlling random variables removed. I was wondering how the partial correlation $ρ_{XY·Z}$ is related to the correlation between $X$ and $Y$ conditional on $Z$?   2. There is a special case for $n=1$. > In fact, the first-order partial correlation (i.e. when $n=1$) is nothing > else than a difference between a correlation and the product of the > removable correlations divided by the product of the coefficients of > alienation of the removable correlations. The coefficient of alienation, and > its relation with joint variance through correlation are available in > Guilford (1973, pp. 344–345). I was wondering how to write the above down mathematically? Thanks and regards!


title: Cox proportional hazard model and non-randomly selected sample

text: Are there any methods to correct bias in Cox proportional hazard model caused by non-randomly selected sample (something like Heckman's correction)? **Background** :   Lets say the situation looks as follows:   \- During first two years all clients are accepted.   \- After those two years a Cox PH model is build. Model predicts how long clients will use our service.   \- Due to the policy of the company from now on only clients with probability of surviving 3 month greater than 0.5 are accepted, the others are rejected.   \- After another two years a new model needs to be built. The problem is that we have target only for accepted clients and using only these clients might cause some serious bias.


title: Fitting multilevel models to complex survey data in R

text: I'm looking for advice on how to analyze complex survey data with multilevel models in R. I've used the `survey` package to weight for unequal probabilities of selection in one-level models, but this package does not have functions for multilevel modeling. The `lme4` package is great for multilevel modeling, but there is not a way that I know to include weights at different levels of clustering. Asparouhov (2006) sets up the problem: > Multilevel models are frequently used to analyze data from cluster sampling > designs. Such sampling designs however often use unequal probability of > selection at the cluster level and at the individual level. Sampling weights > are assigned at one or both levels to reflect these probabilities. If the > sampling weights are ignored at either level the parameter estimates can be > substantially biased. One approach for two-level models is the multilevel pseudo maximum likelihood (MPML) estimator that is implemented in MPLUS (Asparouhov et al, ?). Carle (2009) reviews major software packages and makes a few recommendations about how to proceed: > To properly conduct MLM with complex survey data and design weights, > analysts need software that can include weights scaled outside of the > program and include the "new" scaled weights without automatic program > modification. Currently, three of the major MLM software programs allow > this: Mplus (5.2), MLwiN (2.02), and GLLAMM. Unfortunately, neither HLM nor > SAS can do this. West and Galecki (2013) give a more updated review, and I'll quote the relevant passage at length: > Occasionally, analysts wish to fit LMMs to survey data sets collected from > samples with complex designs (see Heeringa et al, 2010, Chapter 12). Complex > sample designs are generally characterized by division of the population > into strata, multi-stage selection of clusters of individuals from within > the strata, and unequal probabilities of selection for both clusters and the > ultimate individuals sampled. These unequal probabilities of selection > generally lead to the construction of sampling weights for individuals, > which ensure unbiased estimation of descriptive parameters when incorporated > into an analysis. These weights might be further adjusted for survey > nonresponse and calibrated to known population totals. Traditionally, > analysts might consider a design-based approach to incorporating these > complex sampling features when estimating regression models (Heeringa et > al., 2010). More recently, statisticians have started to explore model-based > approaches to analyzing these data, using LMMs to incorporate fixed effects > of sampling strata and random effects of sampled clusters. > > The primary difficulty with the development of model-based approaches to > analyzing these data has been choosing appropriate methods for incorporating > the sampling weights (see Gelman, 2007 for a summary of the issues). > Pfeffermann et al. (1998), Asparouhov and Muthen (2006), and Rabe-Hesketh > and Skrondal (2006) have developed theory for estimating multilevel models > in a way that incorporates the survey weights, and Rabe-Hesketh and Skrondal > (2006), Carle (2009) and Heeringa et al. (2010, Chapter 12) have presented > applications using current software procedures, but this continues to be an > active area of statistical research. Software procedures capable of fitting > LMMs are at various stages of implementing the approaches that have been > proposed in the literature thus far for incorporating complex design > features, and analysts need to consider this when fitting LMMs to complex > sample survey data. Analysts interested in fitting LMMs to data collected > from complex sample surveys will be attracted to procedures that are capable > of correctly incorporating the survey weights into the estimation procedures > (HLM, MLwiN, Mplus, xtmixed, and gllamm), consistent with the present > literature in this area. This brings me to my question: does anyone have best practice recommendations for fitting LMMs to complex survey data in R?


title: Estimating correlation between random slopes and random intercepts using the lme4 package in R

text: For answering my research question I am interested in the correlation between the random slopes and random intercepts in a multilevel model, estimated using the R library lme4. The data I have is: Y (test-scores from students), SES (socio-economic status for each student) and schoolid (ID for each school). I am using the following syntax to estimate random intercepts and slopes for the schools:               library(lme4)     model3 <- lmer(Y ~ SES + (1 + SES | schoolid))      The reference I used for this syntax is this pdf: http://www.bristol.ac.uk/cmm/learning/module-samples/5-concepts-sample.pdf On page 19, a similar analysis is described. It is said that by defining the random intercepts and slopes toghether, it is indirectly specified that we want the random intercepts and slopes to covary. Therefore, also the correlation between random slopes and random intercepts is estimated. Basically, exactly what I need for answering my research hyptohesis. However, when I look at the results:                summary(model3)      I am getting the following output:               Linear mixed model fit by REML ['lmerMod']     Formula: Y ~ SES + (1 + SES | schoolid)          REML criterion at convergence: 8256.4          Scaled residuals:          Min      1Q  Median      3Q     Max      -3.1054 -0.6633 -0.0028  0.6810  3.5606           Random effects:      Groups   Name        Variance  Std.Dev. Corr      schoolid (Intercept) 0.6427924 0.80174                 SES             0.0009143 0.03024  1.00      Residual             0.3290902 0.57366           Number of obs: 4376, groups: schoolid, 179          Fixed effects:                  Estimate Std. Error t value     (Intercept) -0.036532   0.060582  -0.603     SES          0.062491   0.009984   6.259          Correlation of Fixed Effects:         (Intr)     SES 0.226       As stated in the output, the correlation between the random slopes and random intercepts equals 1.00. I find this hard to believe. When I call in R:               VarCorr(model3)$schoolid      I am getting the following output which gives the correlations and covariance matrix:                               (Intercept)          SES     (Intercept)  0.64279243 0.0242429680     SES          0.02424297 0.0009143255          attr(,"stddev")     (Intercept)         SES       0.80174337  0.03023782           attr(,"correlation")             (Intercept) SES     (Intercept)           1   1     SES                   1   1      It seems as if the correlation between the slopes and intercepts was set to 1.00 by R. I did not see this in the output of anyone else when I was searching the internet on references on multilevel modelling. Does anybody know what can be the cause of this correlation? Can it be that the correlation is set to 1.00 because otherwise the model would not be identified? Or is it because the variance of the random slopes is so small (0.0009) that the correlation can not be estimated? I have tried to simulate data in order provide the code for a small reproducible dataset. I was however not yet able to reproduce this output by means of simulated data. As far as I have the code I will eidt my post and add the code. Edit: In response to a comment by Roman Luštrik, the following plot:               ggplot(data[1:261,], aes(x = SES, y = Y)) + geom_point() + facet_wrap(~ schoolid) +        geom_smooth(method=lm)      As there are in total 179 schools the plot becomes quite chaotic, therefore I included the first 10 schools only to make it readable: ![enter image description here](http://i.stack.imgur.com/gSPI5.jpg)


title: Is a large control sample better than a balanced sample size when the treatment group is small?

text: I am running an experiment looking at brain volume changes in a rare disorder. We have a small number of patients (n = 8) but a large control group (n = 100). Some colleagues have suggested that a balanced group should be used, ie. selecting 8 control brains, but I find this counterintuitive. Doesn't it make more sense to use the entire control group to obtain a better estimate of population brain volume?


title: Anova: 'Base case' for interactions when regressing upon dummies

text: Say I want to perform a two-way ANOVA, and I want to use the sources A, B . For A I test at levels A1 and A2, for B at B1,B2 and B3. So when I do the analysis I regress on the dummies with, for instance, A1, B1 as base case. So I get the model $$y = \beta_0 + x_{A2}\alpha_2 + x_{B2}\beta_2 + x_{B3}\beta_3$$ where the $x$'s are the dummies. Now, I also want to include interactions, and I read in my statistics book that I should add the dummies $x_{A2}x_{B2}$ and $x_{A2}x_{B3}$ with corresponding coefficients. But what is the base case here? In the first model I have that the "mean $y$" for A1 x B1 is $\beta_0$. And then "mean $y$" for A1 x B2 is $\beta_0 + \beta_2$ etc, because all other dummies will be 0. But when I have the interactions. I will sort of get all interactions $x_{A1}x_{B2}$, $x_{A1}x_{B1}$, $x_{A1}x_{B3}$ as the base case, but does not seem correct.


title: Calculate BIC to determine the optimal number of clusters (k-means clustering)

text: I have a set of data and want to know whether they fall in 1, 2 or 3 groups. I started exploring the question by using k-means in MATLAB. By just looking at the distance from the centroid of each cluster by eye it looks like I might have two clusters, but I would like to be able to back this up with some statistics. I have read that **using BIC** might be a good way. I am new to clustering, which means that when I read all the answers to questions similar to this one that refer to research papers on the topic, I just cannot follow them, and I was hoping that someone can give me some basic information on whether BIC or something else might be suitable, and some MATLAB code (computing BIC or something similar) with explanations for dummies (I am also new to MATLAB), so that I can use it.


title: Hard thresholding a covariance matrix

text: I am new to the concept of thresholding a variance-covariance matrix and am having trouble understanding the exact process. I am following Bickel and Levina (2008) in choosing a hard threshold. What troubles me is their equation number (3) for the threshold operator: $$ T_{s}(M) = [m_{ij}1(|m_{ij}| \ge s)] $$ My interpretation of that equation is that the thresholding operation applies to the diagonal elements of matrix $M$. This doesn't make much sense to me. In a variance-covariance matrix I am not sure why you would want to set any of the variances equal to zero. To be explicit, my questions are:   * Does the thresholding operator apply to the diagonal elements?   * If I only apply the thresholding operator to off-diagonal elements will that result in a bad estimate of the variance-covariance matrix? The context which my problem comes up is I am estimating a probit model with endogenous regressors via generalized method of moments following Wilde (2008). I have a large number of regressors and a number of them are indicator variables. With some specifications of the model the variance-covariance matrix is singular which presents a problem. I am open to any and all solutions but one solution I read about is this thresholding operation. I want to mention that I am going to bundle the estimation of an endogenous probit model via GMM into an R package. I would really appreciate any help on making it robust and useful to the statistical/econometric community.


title: Orthogonalizing predictors for least squares estimation

text: I know that orthogonalization in LS is to avoid inverting _X'X_. The idea behind it is to find variables _Z_ that are orthogonal to each other. Although the process to find those is clear to me, I don't get the way to find the coefficients. The algorithm leaves $Z_1=X_1$ (I leave the first transformed variable like the original), and then calculates the coefficients in the new space of $z,\ \alpha_1=\frac{<Y,Z_1>}{{||Z_i||}^2}$ (in point 4 in the pseudo- code). Next, $Z_i=Z_i-\frac{<Z_j,Z_i>}{{||Z_i||}^2}$ (point 5 in the pseudo code). That is OK. I understand that. The thing is that to get the $\beta$ (coefficients in the original space), (point 10 in the pseudo code): $\beta_m=\alpha_m-\frac{\sum<\beta\cdot X,Z>}{{||Z_m||}^2}$. And it starts with the last $Z$ , so when m=M, the $\beta_m=\alpha_m$ (because $i=m+1$ and the sum ends at $m$). Finally when it reaches $\beta_1$, it's different from $\alpha_1$. But we were doing $Z_1=X_1$, **so I don't understand why $\beta_1\ne\alpha_1$.** Can anyone give me an idea of why? Thanks! The algorithm would be (in Spanish but you can realize what it's doing): > ![enter image description here](http://i.stack.imgur.com/Aj94C.png)


title: Similarities and differences between correlation and regression

text: If I want to investigate how two continuous variables are linked, what is the difference between calculating the _correlation_ coefficient (Pearson's r) versus calculating the (simple linear) _regression_ coefficient? I see people who, if the regression coefficient is significantly different from zero, talk about the two variables as if they are correlated, which is confusing as it suggests that the two coefficients (correlation, regression) are the same thing. Having said that, isn't r a measure of the (regression line) slope anyway? I'm confused!!