title: Proportionality in Bayesian Models: What Is Absorbed?

text: Considering two Bayesian models:   * Poisson Likelihood & Beta Prior: $p(y|\lambda) \sim \text{Pois}(\lambda)$, $p(\lambda) \sim \text{Be}(a, b)$: $$ p(\lambda|y) \propto \lambda^{a-1}e^{-b\lambda} \times \lambda^{y}e^{-\theta} $$ $$ = \lambda^{a+y-1}e^{-\lambda(1+b)} $$ $$ \sim \text{Ga}(a+y, b+1) $$   * Normal likelihood & Normal Prior: $p(x|\theta, \phi) \sim \text{N}(\theta, \phi)$, $p(\theta) \sim \text{N}(\theta_{0}, \phi_{0})$: $$ p(\theta|x) \propto \text{exp}\left\{-\frac{1}{2}(\theta-\theta_{0})^{2}/\phi_{0}\right\} \times \text{exp}\left\{-\frac{1}{2}(x-\theta)^{2}/\phi\right\} $$ $$ = \text{exp}\left\{-\frac{1}{2}\left(\frac{\theta^{2}-2\theta\theta_{0}+\theta_{0}^{2}}{\phi_{0}}+\frac{x^{2}-2x\theta+\theta^{2}}{\phi}\right)\right\} $$ $$ \propto \text{exp}\left\{-\frac{1}{2}\theta^{2}(\phi_{0}^{-1}+\phi^{-1})+\theta\left(\frac{\theta_{0}}{\phi_{0}}+\frac{x}{\phi}\right)\right\} $$ $$ \sim \text{N}(\theta_{1}, \phi_{1}) $$ where: $\phi_{1} = \frac{1}{\phi_{0}^{-1}+\phi^{-1}}$, $\theta_{1} = \phi_{1}\left(\frac{\theta_{0}}{\phi_{0}}+\frac{x}{\phi}\right)$. Why is it, in the first model, the terms ($a$, $b$, $y$, $n$) in the first model are retained, but in the second model, the terms $\left(\frac{\theta_{0}^{2}}{\phi_{0}}, \frac{x^{2}}{\phi}\right)$ are dropped? Finally, once again in the second model, why is it that the posterior mean is equal to $\theta_{1}$ and the posterior variance is equal to $\phi_{1}$? **EDIT** I'm aware that if some term is a multiple of the parameter of interest (e.g. $c\theta$ or even $\frac{\theta}{c}$), then the scalar ($c$, in this instance) is absorbed by the proportionality sign, leaving just the parameter ($\theta$), but, if possible, could somebody explain the general rules of proportionality with regard to exponents? A few examples: $e^{c+\theta}$, $e^{c\theta}$, $\theta^{c+x}$, $\theta^{cx}$ Assuming $\theta$ is the parameter of interest, what terms would a proportionality sign absorb in the above expressions?


title: Omega squared for measure of effect in R?

text: The statistics book I am reading recommends omega squared to measure the effects of my experiments. I have already proven using a split plot design (mix of within-subjects and between-subjects design) that my within-subjects factors are statistically significant with p<0.001 and F=17. Now I'm looking to see how big is the difference... is there an implementation of omega squared somewhere for R (or python? I know... one can dream ;) Searching on the internet for R-related stuff is a pain the _*_ , I don't know how I manage to find stuff with C. thanks!


title: Combining CostSensitiveClassifer with MultiClassClassifier [RWeka]

text: this is my first attempt at posting here. I looked through CrossValidated and found two similar problems without answers (see How to combine WEKA classifiers and Combine MultiClassClassifier and MetaCost in WEKA and make cost matrices for each binary classifier). Hopefully, third time is the charm here. I'm trying to use the ADTree machine learning algorithm, from the RWeka package, and want to pair it with a cost matrix. I can get an ADTree to work and a cost matrix to work, but not at the same time. If you want to play along, you can use this dummy dataset (100 rows, 6 columns):               # study variables     age <- sample(c(18:65), 100, replace=TRUE)     bmi <- sample(round(seq(from=14, to=40, length.out=100),1), 100, replace=TRUE)     department <- sample(c("Sales", "Marketing", "Finance", "Operations"), 100, replace=TRUE)     location <- sample(c("Georgia", "Florida", "South Carolina"), 100, replace=TRUE)     salary <- sample(c(20000:120000), 100, replace=TRUE)     rank <- as.factor(sample(seq(1:4), 100, replace=TRUE))     d <- data.frame(age, bmi, department, location, salary, rank)      Because the variable I'm predicting (rank) has more than 2 levels, I need to use the MultiClassClassifier w/ ADTree. Since ADTree is no longer part of the base RWeka algorithms, I first need to load it manually:               library(RWeka)     list_Weka_interfaces()     WPM("install-package", "alternatingDecisionTrees")     Sys.setenv(WEKA_HOME = "C:/Users/<...>/wekafiles")  ## my directory to save the interface     WPM("load-package", "alternatingDecisionTrees")     WPM("list-packages", "installed")      With the ADTree interface loaded, I can then build and evaluate an ADTree:               ADTree <- make_Weka_classifier("weka/classifiers/trees/ADTree")     MCC <- make_Weka_classifier("weka/classifiers/meta/MultiClassClassifier")     m.adt <- MCC(rank~., data = d,  control = Weka_control(W = ADTree))     evaluate_Weka_classifier(m.adt, class = TRUE)      However, I want to use the following cost matrix with my ADTree:               c4 <- matrix(c(0,0,0,0, 1,0,0,0, 1,1,0,0, 2,2,1,0), ncol = 4, byrow=TRUE)      I wasn't able to figure out how to call the cost matrix within MCC, but was not able to simply call a cost matrix into the evaluation function as my literature suggested I should.               # attempt at calling cost matrix (in the evaluating function)     evaluate_Weka_classifier(m.adt, class = TRUE, cost=c4)      Well, this code runs, but the results are identical to those if no cost matrix was indicated. I know how to use the CostSensitiveClassifer with something an algorithm not requiring MCC, such as J48:               m.csc.j48 <- CostSensitiveClassifer(rank~., data = d,                   control = Weka_control(`cost-matrix` = c4, W = J48, M = TRUE))     evaluate_Weka_classifier(m.csc.j48, class = TRUE)      But I can't get it to play nicely w/ a MultiClassClassifer. If I run the code w/o MCC:               m.csc.adt <- CostSensitiveClassifer(rank~., data = d, na.action = NULL,                   control = Weka_control(`cost-matrix` = c4, W = ADTree, M = TRUE))     evaluate_Weka_classifier(m.csc.adt, class = TRUE)      I'm told RWeka cannot handle multi-valued nominal class (which is expected as I needed the MCC for ADTree). But, I'm not sure how to get them to speak to each other. Neither of these work:               # attempt #1      m.csc.mcc <- CostSensitiveClassifier(rank~., data = d,                   control = Weka_control(`cost-matrix` = c4,                   W = MCC(control = Weka_control(W = ADTree)), M = TRUE))          # attempt #2     m.csc.mcc <- CostSensitiveClassifier(control = Weka_control(`cost-matrix` = c4,                   W = MCC(rank~., data = d, control = Weka_control(W = ADTree)), M = TRUE))      R tells me _"Error in terms.formula(formula, data = data): argument is not a valid model"_ Can anyone help?? Worst case scenario, I suppose I'll have to use another algorithm, but I'd really like to use ADTree.


title: Probability of MSB=1 for a n-bit discrete random variable while we have a noisy estimate for it

text: **Problem** Let $X$ be a $n$-bit random variable and $s$ be the known estimated value for $X$. Our estimation has an error with normal distribution. In other words, $X|s \sim \mathcal{N}(s,\sigma_z^2)$. We want to calculate $p\left(B_1(X)=1|s\right)$ numerically where $B_1(X)=1$ means the first bit of $X$ equals to $1$. **A Solution** In the page 4 of this paper the authors have calculated this probability as follows: $$ p\left(B_1(X)=1|s\right) = \int_{\{x|B_1(x)=1\}} p(x|s) dx \\ = \frac{1}{2} \sum_{i=-\infty}^{+\infty} \left[ erfc\left( {\frac{(i-1/2)d_{min}-s}{\sqrt{2}\sigma_z}} \right) - erfc\left({\frac{id_{min}-s}{\sqrt{2}\sigma_z} }\right)\right] $$ **Questions**   1. What have they done? I can't understand why there is a $\Sigma$ over $erfc$ functions. Apparently CDF of a gaussian distribution is an error function. Maybe that's because $X$ is from a discrete space. In addition, how they have included the $B_1(X)=1$ condition in their answer?   2. Why they have calculated the integral analytically? As $X \in \{0,1\}^n$ (let say $n=8$) and we just want the numeric value of $p\left(B_1(X)=1|s\right)$, we can enumerate over all $x$'s with $B_1(x)=1$ (e.g. 128 possible $x$ for $n=8$), put them in $p(x|s)$ and sum up the probabilities to get $p\left(B_1(X)=1|s\right)$. Is there anything wrong with this approach? **More details:** While my problem is related to channel coding and LDPC codes, I try to explain it without going into much detail: I am trying to use a LDPC encoder and decoder to transfer $X$. LDPC let me to reduce the bitrate by providing the decoder an estimate of incoming $X$. Let say I have found a way to calcuate $s$ as an estimation of incoming $X$ (we can imagine we observe $X$ from a _virtual noisy channel_ ). The problem is the LDPC decoder cannot benefit from this estimation directly, but I have to feed it with _the probability_ of each bit of $X$ being 0 or 1 given $s$. In other words, the $s$ is not enough by itself to guess incoming $X$ and we should also consider the dependence between $X$ and $s$. In practice, it is supposed that our estimation has a gaussian noise. In other words, we observe $X$ from a virtual channel with gaussian noise. Therefore, we can say that $X|s \sim \mathcal{N}(s,\sigma)$ which means: _knowing $s$, $X$ is distributed around $s$ with variance $\sigma^2$_. Now, all I have to do is to calculate the probability of each bit of $X$ being 1 or 0 given $s$ (instead of feeding the decoder directly with $s$ bits). Therefore, for MSB I should calculate $p(B_1(X)=1|s)$, for the next MSB I should calculate $p(B_2(X)=1,B_1(X)|s)$ in similar way, etc. One obscure point is the $d_{min}$ parameter. As I know, $d_{min}$ is the minimum hamming distance between two coded symbols (i.e. minimum hamming distance of LDPC encoder output for all possible $X$'s). I think what `Dilip Sarwate` said about $d_{min}$ is correct but I am not still sure. In conclusion, I know this problem is a little complicated and needs a good knowledge on both probability theory and channel coding to be fully understood. Anyway, I just want to know why author of mentioned paper have calculated the probability analytically? From probability point of view, can I just compute the summation over all possible $X$'s that meet $B_1(X)=1$ condition (as I asked above in Question #2).


title: When is a repeated measures ANOVA preferred over a mixed-effects model?

text: In response to this question, regarding whether my design where I randomly presented participants with pictures from different categories was an example where I should use a repeated measures ANOVA, I got the answer that I should use a mixed-model instead, with one of the reasons being that I have two forms of dependencies: for subjects and for categories. My question is now: Isn't it always the case that you have two dependencies in this way when doing this type of repeated measures design? That is, under what circumstances would a repeated-measures ANOVA be preferable to a mixed-effects modeling approach and why?


title: Multi task learning

text: I have a dataset where all observations are measured several times and reported outcomes correspond to those measurements. In other words, my set of data points looks like $\{x_i, y_{i_1}, y_{i_2}, ... y_{i_n}\}$ I would like to make a predictive model to predict $y$ from $x$. What should I try to predict? The average of $y_i$? I've heard something about multi task learning, where each observation is considered as $n$ different regression problems and all problems share the same kernel. Can anybody elaborate on this?


title: Using caret train function to perform lm with 10 repeated 10-fold cross validation

text: I'm using the caret package at the moment to perform different forms of analysis on my dataset. I did the following to do 10-fold cross validation, 10 times:-               >library(caret)     >cadets <- read.csv("pathway/CADETS.csv")     >ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 10, savePred = T)     >model <- train(RT..seconds.~., data = cadets, method = "lm", trControl = ctrl)     >model$pred      I took the values obtained from this, and plotted graphs for each repeat, whereby I plot the predicted values against the observed value for each fold. However, the dataset I'm working has been worked on previously by someone else, and the results I have obtained for the linear regression model is completely different to what they've got. I've attached a picture which shows the star difference between our models. I got these values for the error of the model:-               Resampling results            RMSE  Rsquared  RMSE SD  Rsquared SD       2250  0.1       4620     0.118          Whilst they got a RMSE of 179.9618. The same method they had applied to the dataset is the same thing I have done, the only difference being that they worked in matlab (which I have no idea how to use, which is why I've done it in R). I was wondering whether any of you can see maybe what I'm doing wrong in my coding, or what parameters I could tweak maybe that would help make a more accurate model? I'm assuming of course that what he has done is correct. I've also uploaded the file here (the dataset I'm using http://www.sendspace.com/file/wt9vvd) if any of you want to have a look at the dataset I'm using as well. I'd appreciate any help! I know I can help improve my model by doing feature selection, but I just need to know why the multiple linear regression model I have to begin with is way off the reference one. Thank you! ![enter image description here](http://i.stack.imgur.com/33KXX.jpg)


title: Beta distribution fitting in Scipy

text: According to Wikipedia the beta probability distribution has two shape parameters: $\alpha$ and $\beta$. When I call `scipy.stats.beta.fit(x)` in Python, where `x` is a bunch of numbers in the range $[0,1]$, 4 values are returned. This strikes me as odd. After googling I found one of the return values must be 'location', since the third variable is 0 if I call `scipy.stats.beta.fit(x, floc=0)`. Does anyone know what the fourth variable is, and if the first two are $\alpha$ and $\beta$?


title: Finding subsets with Branch-and-bound techniques

text: Hello I am trying to find a Matlab version of a function which performs best subset selection for least squares regression using leaps and bounds (branch and bound) approach. Like in the paper by Furnival and Wilson (1974)*. The R equivalent of this is the `leaps` package. Any suggestions for this function in Matlab? * Furnival, G. M. and Wilson, R. W., Jr. (1974)   _Regressions by leaps and bounds_   Technometrics, **16** , 499–511.   Reprinted in Technometrics, **42** , 69–79.


title: Polynomial regression using scikit-learn

text: I am trying to use scikit-learn for polynomial regression. From what I read polynomial regression is a special case of linear regression. I was hopping that maybe one of scikit's generalized linear models can be parameterised to fit higher order polynomials but I see no option for doing that. I did manage to use a Support Vector Regressor with a poly kernel. That worked well with a subset of my data, but it takes much to long to fit larger data sets so I still need to find something faster (even if trading some precision). Am I missing something obvious here?


title: Sufficient and necessary condition for the integrability of a random variable

text: Could anyone give me a hint on the following problem? Thanks! Let $f$ be a random variable over a probability space $(\Omega,A,P)$. Show that $f$ is integrable $\iff $ $\sum\limits_{k=1}^\infty P(\{|f|>n\}) < \infty.$ Note: Consider the sets $A_n:=\{n-1<|f|\le n\}$, $n\in \mathbb{N}$ and the random variables $L:=\sum\limits_{k=1}^\infty (n-1)1_{A_n}$ and $R:=\sum\limits_{k=1}^\infty n1_{A_n}$,$n\in\mathbb{N}$


title: Cluster centroids visualization

text: I've some data that I've clustered and calculate the similitude within my cluster centroids to spot which class centroids are close to which ones. My goal was to visualize this as a map of points. I thought that the best way was to represent the data as a complete graph where each node represent a centroid and the edges the similitude within the centroids. Plug this into a graph visualization software like gephi let the nodes pull an push each other, get a pretty picture and job done. But i couldn't make it work the graph collapsed all the time. Maybe is just that my data is supposed to collapse, but here's how the similitude within my 36 nodes look like and to me, it seems that the graph visualization I've in mind should go somewhere. The other option is that i'm not familiar with the gephi software and i did not manage to make it work. Do you have any suggestions? software? visualization technique? how would you show a more intuitive visualization of the similitude matrix? thanks![enter image description here](http://i.stack.imgur.com/O8Jl3.png)


title: Exact distribution of sample mean of a continuous distribution

text: Please first note that I do know about the central limit theorem but I wish to derive an exact expression for the sample mean for any continuous distribution with probability density function $f(x)$. I thought I derived it correctly here but there is some mistake I can't spot. Can anyone tell me where I have made a mistake? I take a simple case to derive the sample mean for with only 2 samples. The probability density function $f(x)$ is defined as: $$f(x)=\frac{x}{50}$$ The maximum value of $x$ is $10$ and the minimum value of $x$ is $0$. Note that: $$\int_0^{10} f(x) dx =1$$ The samples are random variables $X_1$ and $X_2$ $$\bar{x}=\frac{x_1+x_2}{2}$$ If the sample mean is to equal $\bar{x}$ then given that the first sample is $x_1$ the second sample $x_2$ must equal $2\bar{x}-x_1$. But if $x_1$ is too large or too small then a value of $x_2$ may not exist in the range [0,10]. We know $$x_2 \leq 10 \quad \therefore 2\bar{x}-x_1 \leq 10$$ $$\implies x_1 \geq 2\bar{x}-10$$ And $$x_2 \geq 0 \quad \therefore x_1 \leq 2\bar{x}$$ $$\implies 2\bar{x}-10 \leq x_1 \leq 2\bar{x}$$ The probability that $X_1$ takes the value $x_1$ and $X_2$ takes the value $2\bar{x}-x_1$ is $f(x_1)f(2\bar{x}-x_1)$ To account for all cases of $X_1$ that could result in $\bar{X}=\bar{x}$ we integrate with respect to $x_1$ $$\int_{2\bar{x}-10}^{2\bar{x}}f(x_1)f(2\bar{x}-x_1)dx_1$$ I did this integration and it ended up as $$\frac{1}{7500}(1000-570\bar{x}-120\bar{x}^2)$$ The sample mean must lie inside $[0,10]$ so the integral of all the eventualities of the sample mean over the range $[0,10]$ should equal 1. But: $$\int_0^{10} \frac{1}{7500}(1000-570\bar{x}-120\bar{x}^2)=-7.8$$ So where have I made an incorrect assumption?


title: Proof of a PD covariance matrix for conditional Gaussian

text: I was looking at the formula for the conditional covariance of a partitioned matrix. I understand the intuition behind the equation for the conditional covariance, but I'm not sure how to show that the covariance matrix is positive definite. $$\text{For }\bf{x} \sim N(\bf{\mu},\Sigma),$$ $$ \bf{x}= \begin{bmatrix} \bf{x}_a \\ \bf{x}_b \end{bmatrix} $$ $$ \bf{\mu}= \begin{bmatrix} \bf{\mu}_a \\ \bf{\mu}_b \end{bmatrix} $$ $$ \Sigma = \begin{bmatrix} \Sigma _{aa} & \Sigma _{ab} \\ \Sigma _{ba} & \Sigma _{bb} \end{bmatrix} $$ $$ \bf{x}_a|\bf{x}_b \sim N(\bf{\mu}_{a|b},\Sigma_{a|b}) $$ where $$ \bf{\mu}_{a|b}=\bf{\mu}_a+\Sigma _{ab}\Sigma^{-1}_{bb}(\bf{x}_b-\bf{x}_a) $$ $$ \Sigma_{a|b}=\Sigma_{aa}-\Sigma_{ab}\Sigma^{-1}_{bb}\Sigma_{ba} $$


title: What methods exist to measure how similar a set of sampled distributions are to each other?

text: I have a set of distributions that are sampled. How do I measure their similarity between themselves or as a group/set? So, given a finite domain for a sampled function$F_j(x)$, with $x_i \in {1,\ldots,N}$ and $j \in {1, \ldots, J}$. I would like a measure $\cdot(F_{1,J})$ for how similar the distributions are, or test of whether whether they are statistically close enough, eg. a test producing some p-value.


title: How to forecast time series with the help of other training time series in R?

text: Consider such a task: I have the figures on last year's bookselling of a certain book store, and the information of the first half of this year too. Now I want to predict the figure at the end of this year. I've checked R's `forecast` package. But it seems to only work on a single time series each time. How can I integrate the information from last year? Thank you.


title: Is $R^2$ value valid for insignificant OLS regression model?

text: I am interested in stating that ___ % of the variance in Y is explained uniquely by $X_1$ and ___ % is explained uniquely by $X_2$.   * Is there some way to obtain this from a multiple regression model, or do I need to obtain adjusted $R^2$ values from a series of residual regressions (sensu Legendre & Legendre et al.)?    * That being asked, is $R^2$ value valid for insignificant OLS regression model?   * Does the value of $R^2$ depend on a statistical test?


title: Propagation of error using 2nd-order Taylor series

text: I am reading a text, "Mathematical Statistics and Data Analysis" by John Rice. We are concerned with approximating the expected value and variance of the random variable $Y$. We are able to calculate the expected value and variance of the random variable $X$ and we know the relation $Y = g(X)$. So, it's possible to approximate the expected value and variance of $Y$ using the Taylor series expansion of $g$ about $\mu_X$. On page 162, he lists 3 equations.   1. The expected value of $Y$ using the 1st-Order Taylor series expansion. It is: $\mu_Y \approx g(\mu_X)$. This is referred to later in my question as $E(Y_1)$.   2. The variance of $Y$ using the 1st-Order Taylor series expansion. It is: $\sigma_Y^2 \approx \sigma_X^2 (g'(\mu_X))^2$. This is referred to later in my question as $Var(Y_1)$.   3. The expected value of $Y$ using the 2nd-Order Taylor series expansion. It is $\mu_Y \approx g(\mu_X) + \frac12 \sigma_X^2 g''(\mu_X)$. This is referred to later in my question as $E(Y_2)$. Note that there are two different expressions for $Y$ because we are using two different orders in the Taylor series expansion. Equations 1 and 2 refer to $Y_1 = g(X) \approx g(\mu_X) + (X-\mu_X)g'(\mu_X)$. Equation 3 refers to $Y_2 = g(X) \approx g(\mu_X) + (X-\mu_X)g'(\mu_X) + \frac12 (X-\mu_X)^2 g''(\mu_X)$. Note that specifically the equation for $Var(Y_2)$ is not given. Later, the author seems to use the equation for the variance of $Y_1$ (Equation 2), when in fact he is referring to expected value of $Y_2$ (Equation 3). This seems to imply $Var(Y_2) = Var(Y_1)$. I have tried to calculate by hand $Var(Y_2)$, and I am getting a somewhat complicated expression. Here is my work (I stopped because at the end I am getting $X^3$ terms in the expectation): $$ \begin{aligned} Var(Y_2) &= E[( g(\mu_X) + (X-\mu_X)a + \frac12 (X-\mu_X)^2 b - g(\mu_X) - \frac12 \sigma_X^2 b )^2] \\ &= E\left[((X-\mu_X)a + \left(\frac12 (X-\mu_X)^2 - \frac12 \sigma_X^2)b\right)^2\right] \\ &= E\left[(ca + \left(\frac12 c^2 - \frac12 \sigma_X^2)b\right)^2\right] \\ & = E[c^2 a^2 + ca(c^2 - \sigma_X^2)b + \frac14(c^2-\sigma_X^2)^2 b^2] \\ & = E[(X^2 - 2X \mu_X + \mu_X^2)a^2 + (X-\mu_X)a((X^2 - 2X\mu_X + \mu_X^2) - \sigma_X^2)b \\ & + \frac14((X^2 - 2X\mu_X + \mu_X^2)-\sigma_X^2)^2 b^2] \end{aligned} $$ Note that in the above equations, $a = g'(\mu_X)$, $b = g''(\mu_X)$, and $c = X-\mu_X$. What is $Var(Y_2)$? Thanks.


title: AdaBoost on a continuum of base classifiers

text: A tutorial on AdaBoost suggests that AdaBoost can be applied to a continuum of classifiers (at the bottom of the first page). Does it mean to simply discretize the classifiers, for example, which are parameterized by one or two real number, or refer to a generalized version of AdaBoost theory, which may be capable of coping with the abstract classifier "space" (though everything should be discretized when put into practice)? I did some googling, but failed to find what I need. Can anyone provide some references? Thank you~


title: Is it valid to assign observations partially to both test and train in n-fold cross validation?

text: I'm looking to use n-fold cross validation for selecting meta-parameters for fitting a model to a dataset. However, dropping observations entirely from the learning-set while fitting the model to each of the folds may create problems for the model fitting. I was wondering whether it would be valid to perform the cross-validation on a weighted basis, e.g. the observations selected for the fold would receive a weight of say 0.95 and the ones excluded from the fold would receive a weight of say 0.05. When assessing predictive performance one then would weight prediction errors according to the complement of these weights, e.g. (1-0.95) for those 'participating heavily' and (1-0.05) for those 'participating lightly'. Is there any merit in this intuition? Such a process is similar to Fay's Balanced Repeated Replication (Survey sampling variance estimation). I have never encountered it in relation to cross-validation. Any literature pointers are appreciated.


title: Random effect nested under fixed effect model in R

text: Looking for help to create a model for my data gathered from a reciprocal transplant experiment: I have 2 populations of fish (A and B) and 2 temperature regimes (warm and cold) that were crossed with each other (= 4 treatment groups- A:warm, A:cold, B:warm, B:cold). Each treatment group consists of two replicate tanks (8 tanks total). There are 100 fish in each tank. The response variable of interest is growth rate and I have recorded growth rate for each of 800 fish.I want to test the effects of population and temperature and their interaction. I also want to include tank effect. I think I should nest tank effect within the interaction effect (Pop:Temp which is equivalent to the treatment group) and include nest effect as a random effect. So I should have Population (fixed effect) and Temp (fixed effect) and tank effect (random effect) nested within the interaction of Pop:Temp. Is it possible to create a random effect nested within a fixed effect interaction? If so, how do I do this in R?


title: The upper bound of the training error of AdaBoost

text: I am reading an overview of AdaBoost written by Schapire, which calculates the upper bound of the training error in Eq. (5), section 3. In fact, it states that $$\prod_{t}Z_t=\prod_{t}\left[2\sqrt{\epsilon_t(1-\epsilon_t)}\right]$$ with $$ \begin{align} Z_t=&\sum_{i}D_t(i)\exp(-\alpha_ty_ih_t(x_i))\\ \alpha_t=&\frac{1}{2}\ln\frac{1-\epsilon_t}{\epsilon_t} \end{align} $$ where $t$ denotes the iteration times, $i$ the index of data points, $h_t$ the classifier selected at the $t$th round, $\alpha_t$ the weight of $h_t$, and $D_t$ satisfies $$D_{t+1}(i)=\frac{D_t(i)\exp(-\alpha_ty_ih_t(x_i))}{Z_t}$$ In order to simplify $\prod_{t}Z_t$, I've tried: $$ \begin{align} Z_t=&\sum_{i}D_t(i)\exp(-\alpha_ty_ih_t(x_i))\\ =&\sum_{i}D_t(i)\left(\sqrt{\frac{1-\epsilon_t}{\epsilon_t}}\right)^{-y_ih_t(x_i)}\\ =&\sum_{y_i\ne h_t(i)}D_t(i)\left(\sqrt{\frac{1-\epsilon_t}{\epsilon_t}}\right)+\sum_{y_i= h_t(i)}D_t(i)\left(\sqrt{\frac{\epsilon_t}{1-\epsilon_t}}\right) \end{align} $$ But I don't what to do next. Can anyone give hints?


title: Modeling a 1D random walk with nonconstant probabilities

text: I have (what I will term, for lack of a better word) a random walk that has a particular property: it tends to be right of the origin some fraction _k_ of the time and left 1- _k_ of the time (and on the origin a negligible fraction). So although the probabilities of going up or down are close to 1/2, there is a central tendency as a result of this propensity. (Note: $k\neq1/2$ else I would use a Poisson process.) It would be best if the process was memoryless in the sense that while the probability can vary (slightly) from 1/2 based on the current position, the count so far is ignored. Thus, it would perform the same at 0 whether it had just spent the last million above 0 as if below. Is such a model possible within my constraints? I want to analyze the probability that the process is, say, left of the origin over an interval of varying lengths, given that the process starts at the origin. I'm looking for a model as natural (parsimonious) as possible in hopes that I can compare my measurements to the model to see if the event is as rare as it (naively) appears to be. The combination of it being a walk (so that it can 'dig itself into a hole' as it were) and its unusual central tendency have confused my intuition as to what is normal... **Edit:** For example, one possible model that seems to fit my requirements: some (N, p) where for position n > N go left with probability p and right with probability 1-p; for position n < N go left with probability 1-p and right with probability p; for n = N use 1/2 either way. (Think of p as just slightly more/less than 50%.) The process doesn't know its score but by choosing N sufficiently far from the origin and p close enough to 50% I get my requirements. It's not quite parsimonious, having one too many degrees of freedom (multiple choices for a given _k_ ) and it really doesn't fit my intuition as to how the process works, but maybe it inspires some better ideas. The underlying problem, for those interested or looking for inspiration on the problem: this models the crossing point where $\pi(n)-\operatorname{li}(n)=0$ after Rubinstein & Sarnak 1994, where _k_ is about $2.6\times10^{-7}$.


title: Analysis of variance with Weibull or Gamma distributions

text: I have been trying to find a method to analyse variance on Weibull and/or Gamma distributions but a Google search for               anovar  Weibull "gamma distribution"      yields nothing helpful. The data I have cannot be fitted to a normal distribution but fits a Weibull or Gamma distribution quite well.


title: Deriving the in-sample error for linear model from the elements of statistical learning

text: From the elements of statistical learning, it was claimed that $$ \frac{1}{N}\sum_{i=1}^N ||h(x_i) ||^2 \sigma^2_\varepsilon= \frac{p}{N}\sigma^2_\varepsilon$$ where $h(x_i) = X(X^TX)^{-1}x_i$. Can someone show me how to prove this ? Thanks This came from the image below ![enter image description here](http://i.stack.imgur.com/3658P.png)


title: How to handle more than one dependent variable (categorical) in logistic regression?

text: My data consists of 6 independent variables (continuous and categorical) and 8 dependent variables on likert scale (categorical). I want to use multinomial logistic regression to find out the relationship between dependent and independent variables. I reduced the 8 dependent variables by principal component analysis (factor analysis) into 2 factors. My question is that could I use the factor results of SPSS or average of factors as dependent variables in logistic regression?


title: Does R have post hoc tests robust to unequal sample sizes/population variances?

text: While reading _Discovering Statistics Using R_ pp. 431-432, Dr. Field says that > "There are a variety of tests designed to deal with these situations > [multiple comparison procedures with unequal group sizes &/or different > population variances], none of which are implemented in R. Hochberg's GT2 is > one such test and is worth mentioning because it is not implemented in R... I didn't find anything in a Google search that indicated otherwise neither have I found an alternative post-hoc test for this situation. So, does anybody know if R has post hoc tests robust to unequal sample sizes/population variances?


title: Are there examples of more informative PCA plots?

text: I am often disappointed with PCA plots in the scientific literature. Typically PCA plots do not provide a breakdown of the variables and their weights, just something like PCA1 (70% variance explained), PCA2 (10% variance explained). How could one tell which variables are strongly loaded into a component? Are there PCA visualizations that can provide more insight into the data?


title: Significance of a 1 state Hidden Markov Model

text: I've been training different observation sequences to obtain different HMMs corresponding to each observed data. Something intriguing is that I get one observation sequence represented by 1 state. What does that represent? Does this mean that it's somewhat repeating a process?


title: Kernel density estimation (CDF) with Epanechnikov kernel in Matlab

text: The code for Kernel density estimation was given in a recent CrossValidated question by Julio Miguel Galvez entitled "Kernel density estimation with Epanechnikov kernel in Matlab", as follows. How can we derive CDF values by using Epanechnikov kernel from below codes?               % Kernel density estimation using mixtures of normally distributed random % variables.     % Set parameters for data generating process.      n=1000; % Number of observations.       data1 = rand(n,1); % Returns data with 1000 observations.     % Initializing the zero matrix of data from the mixture.      data2=zeros(n,1);     % % Generate the data. (mixture of normal random variables)      data2 = (data1 < ones(n,1)/3).*normrnd(2,1,n,1) + ...      (data1 > ones(n,1)/3).*normrnd(-2,1,n,1);     % Generating the density (using an Epanechnikov kernel.)      % Setting the parameters.      s=std(data2); % standard deviation of the data      h=1; % bandwidth parameter     % Evaluate the kernel at all x's in the domain.      x=linspace(-3,3,1000); p=ones(n,1); xi=p*x;      % matrix for the x's where we evaluate the kernel.     data2i=data2*ones(1,n); % matrix for each of the data points in the kernel.      u=(xi-data2i)/h; % matrix of u's.      absu=abs(u); % absolute value of u.      I=(absu<=1);      f=(.75/h)*(1-u.^2).*I;     ff = mean(f,1);      plot(x,ff);


title: Normal-like distribution over a bounded area

text: Is there a distribution that resembles the gaussian (normal) distribution, but such that it's probability density is nonzero only over a defined segment. The question emerged when I tried to model the 'bullet spread' within a circle. Gaussian distribution works fine, but there is always a chance that the bullet would hit outside the circle. So I'd like to find a distribution very similar to Gaussian, but with property that the probability outside the defined segment (or circle) is zero. EDIT: Yes, actually I mean a disk, not a circle. EDIT: And yes, I need only a one-dimensional distribution (along the radius of a disk) which will be circular-symmetrical (not dependent on angle).


title: Can I combine several ordinal questions to create an index and/or composite measure?

text: I am using a survey that contains several questions about various dimensions of performance for policy research institutes. Here, performance in the policy arena is unpacked into things like:   * quality of research,    * overall ability to engage with policy stakeholders,    * quality of recommendations, overall support to, and influence on, policymakers,    * etc.  Each question has a typical 5-level ranking (i.e., ordinal response from "strongly disagree" to "strongly agree", or "very bad" to "very good"). I am thinking of pursing either of two options:   1. Creating a sort of composite measure, where all dimensions of performance are aggregated together so as to have one composite variable of "performance". I could then use this composite variable as the dependent variable (perceptions of performance).   2. Combining 20 of these questions to create a performance index. For each dimension, responses range between 1 and 5. So total scores on the index would thus range from a minimum of 20 up to a maximum total of 100 points. This index could also be used as a dependent variable. Or perhaps only for descriptive statistics.  Does this make sense? Any advice and reference would be greatly appreciated.


title: How to specify uncorrelated random slopes in lmer() syntax?

text: I am currently fitting a mixed effects model to some experimental data as follows:               model <- lmer(Y ~ X + M + (1+X+M|Subject), data=mydata)      The meaning of the variables is not so important here, but $X$ is the predictor of interest while $M$ is a (suspected) mediating variable. All variables are continuous and measured within-subjects. Now the question concerns the random slopes in this model. The above syntax specifies fully correlated random effects. However, I would like to remove the correlation between the two random slopes ($X$ and $M$) without removing the correlation between the random slopes and the random intercept. Initially, I attempted the following code:               model <- lmer(Y ~ X + M + (1+X|Subject) + (1+M|Subject), data=mydata)      This does produce uncorrelated random slopes but lmer() now estimates a random subject intercept _both_ for $X$ and $M$. I am not sure this is correct (or what I require), because I am now forced to introduce an extra variance parameter (simply for removing another one). Is there any way to specify a single subject intercept _and_ uncorrelated random slopes for $X$ and $M$?


title: Proper methods for this case of repeated measures?

text: I am measuring hormone changes in animals. I have three groups: Groups A, B, and C. In each group, blood was collected from the SAME individuals for hormone assay at 0, 15, 30, and 60 min. This was done in males and in females and I would like to keep in mind sex differences. My question: do the changes we see in hormones over time differ between groups A,B, and C? And if so, how? Does the variation of hormone change over time differ between groups A,B, and C? I would like to analyze in one big model, but I'm not exactly confident in the test(s) I'm using. Any advice?


title: How to translate the output from an lm() fit with a cubic spline into a regression equation

text: I have some code and output, and I would like to construct a model. I don't know how to construct a model using this output:                require("splines")      x   <- c(0.2,   0.23,   0.26,   0.29,   0.33,   0.46,    0.53 )      y   <- c(0.211, 0.2026, 0.2034, 0.2167, 0.2177, 0.19225, 0.182)      fit <- lm(y ~ ns(x,3))      summary(fit)      Note that `ns()` generates the B-spline basis matrix for a natural cubic spline. Thus this model regresses `y` against a B-spline for `x` using three degrees of freedom. What would the equation for such a model look like?


title: Likelihood of multiple event times modeled as independent Poisson processes

text: I am modeling three events A, B, and C as Poisson processes with rates $\lambda_A$, $\lambda_B$, and $\lambda_C$ and I would like to calculate the likelihood of observing some data given my model. A distinction that I have come across is the difference between an event happening at a specific time $t$, $P\left(T\left(E\right) = t \right)$, versus by time $t$, $P\left(T\left(E\right) \leq t \right)$, where $T\left(E\right) = E_n - E_{n-1}$ refers to the inter-arrival time between two consecutive occurrences of the event. For a single event with rate $\lambda$, the probability of these scenarios are the exponential likelihood function $P\left(T\left(E\right) = t \right) = \lambda e^{-\lambda t}$ and an integration over it $P\left(T\left(E\right) \leq t \right) = \int_0^t \lambda e^{-\lambda \tau} \, \textrm{d}\tau = 1 - e^{-\lambda t}$. It follows that the odds of the event not occurring by $t$ is $P\left(T\left(E\right)\gt t\right) = 1 - \left(1 - e^{-\lambda t}\right) = e^{-\lambda t}$. Returning to my three-event model, the probability of observing event A at time $t$ but not events B nor C by that time is $P\left(T\left(A\right) = t \cap T\left(B\right) \gt t \cap T\left(C\right) \gt t \right) = \lambda_A e^{-\lambda_A t} e^{-\lambda_B t} e^{-\lambda_C t} = \lambda_A e^{-\lambda t}$ where $\lambda = \lambda_A + \lambda_B + \lambda_C$, the overall rate of events. (It is worth noting that this calculation corresponds with the steps for simulating under multiple Poisson processes.) However, I see two non-equivalent approaches to calculate the probability of observing event A by time $t$ but not events B nor C by that time, $P\left(T\left(A\right) \leq t \cap T\left(B\right) \gt t \cap T\left(C\right) \gt t \right)$. The first method I see is to explicitly evaluate each expression, as I did previously; that is, $P\left(T\left(A\right) \leq t \cap T\left(B\right) \gt t \cap T\left(C\right) \gt t \right) = \left(1 - e^{-\lambda_A t}\right) e^{-\lambda_B t} e^{-\lambda_C t} = e^{-(\lambda-\lambda_A)t} - e^{-\lambda t}$. The second is to integrate over all possible times: $P\left(T\left(A\right) \leq t \cap T\left(B\right) \gt t \cap T\left(C\right) \gt t \right) = \int_0^t \lambda_A e^{-\lambda \tau} \, \textrm{d}\tau = \frac{\lambda_A - \lambda_A e^{-\lambda t}}{\lambda}$. Both approaches seem equally valid to me. Which calculation should I be using, if either, and why? Thanks in advance for your help; I really appreciate it.


title: The vcov function cannot be applied?

text: I originally asked a question about the delta-method in the context of the hyperbolic distribution. I got an answer there, which is useful, except that it says I should apply the `vcov` function to my fitted function what is unfortunately not possible. The `vcov` was there to give me the variance- covariance matrix of the parameter $\zeta,\pi,\delta$. Now I want to calculate this variance-covariance manually. According to this question of me, I know how to get the variance-covariance matrix out of the hessian. Unfortunately according to this this is the hessian for `pi, log(zeta), log(delta), and mu` and not for the desired `pi, zeta, delta, and mu`. So my basic problem is how to get the variance-covariance matrix of the parameter estimates of the hyperbolic distribution (in my original post this was the Sigma solved with the vcov function which does not work)? I know how to get the variance-covariance matrix out of the hessian, but I have the wrong hessian. So what should I do now? Is it possible to transform the hessian from `pi, log(zeta), log(delta), and mu` into the hessian of `pi, zeta, delta, and mu` and how is this done? Or is there an easier way to get vcov running? I am really stuck here and don't know how to continue. **My main question is: How to get the answer of COOLSerdash in thispost working? The vcov function is not working and also the hessian cannot be used, see my description of the error message and my comments to the answer.**


title: Show "Does it affect?" with hypothesis test

text: I am trying to learn ttest and ANOVA for my project (i am a computer science student). I need to know that my subjects answers "effected" by the feedback they get from the first question.   What should i use to show "sth did effect or not sth"? ANOVA or ttest? I am sorry if it sounds very stupid. Thanks in advance. My data is: A person answers a question and receives a feedback as "good/bad" than he/she answers another question. Does the second answer effected by the feedback or not? That is what i want to learn.


title: Filtering a dataframe in R based on multiple Conditions

text: I am new to using R. I am trying to figure out how to create a df from an existing df that excludes specific participants. For example I am looking to exclude Women over 40 with high bp. I have tried several times to use the subset () but I cannot find a way to exclude using multiple criteria. Please Help!


title: Using selected features from a wrapper algorithm to train another model

text: I was wondering if it can be useful to use selected features from a wrapper algorithm (for example SVM-RFE) to train another classification model like k-NN or Linear regression.


title: Statistical Tests with Data Having Measurement Errors

text: For example, can the results of the t test on y1 and y2 be interpreted in the usual way (i.e., like the results of the t test on x1 and x2)? If not, how should I go about testing whether or not y1 and y2 are drawn from the same population?               #Error-free measurements follow a N(10,1)     #I don't know these values     x1 <- rnorm(100,10,1)     x2 <- rnorm(100,10,1)          #Do a t test on the values I don't have     t.test(x1,x2)          #Add 10% relative measurement error     #These are the values I know, and I also know     #that they have a relative uncertainty of 10%      y1 <- x1 + rnorm(100,0,0.1*x1)     y2 <- x2 + rnorm(100,0,0.1*x2)          #Do the t test on the measurements with error     t.test(y1,y2)


title: Should I re-scale the PCA score before classification?

text: This question is about using PCA as a dimension reduction method before feeding the data into a classifier. It's a common procedure to use PCA for a data set which contains a large number of features, and to only use the first several PCA-scores instead of the original features. **My question is** : After the PCA score has been extracted should I need to re-scale them ? (as the scores are in descending order...and can be in different magnitude)


title: List of graph layout algorithms

text: I am a newbie in visualizations. Nowadays, I regularly come across graph-based data that I need to visualize. However, different type of data requires different layouts. And as I do not have enough experience, I am not aware of the different graph layout algorithms. I would like to have a list of graph layout algorithms with pseudo-code or code (in any language) that I can refer to when implementing some graph.


title: How to determine probability of score when selecting from pool of values?

text: The question seems simple but I just can't solve it. I have twelve test scores and three are to be picked at random to determine the overall test grade. How can I calculate the probability of my final grade based on three randomly picked tests from the pool? The test scores are as follows:               100=4     33=2     66=6          100 33 66 66 66 100 66 66 100 66 33 100      Would running a simulation be the simplest solution? Is there a good solution?


title: Strict convexity of Ridge vs Convexity of LASSO

text: Is there any intuition why the ridge regression is strictly convex, while the LASSO is only convex? Does it have to do with the "corners" of the L1 regularization?


title: Dimensionality reduction method for uncorrelated data?

text: In neuroscience it is very common to measure the reaction time (RT) of subjects. Based on the RT conclusions can be made about personal working memory capability, IQ etc. So I have such data, which came from some neuroscience experiment. In this data I have 180 subjects from two groups (90 subjects in each group) for simplicity let say boys and girls, and I have a set of 500 RTs for each subject. The difference between the mean RT of each group is statistically significant (t-test). Now, I would like to construct a classifier which will learn the differences between those two groups. I want the classifier to classify new subjects only based on their RT - and I'm talking specifically about classification of RT. **(A)** I'm searching for dimensionality-reduction method which would fit to this case. ( ~~I have tried PCA and it doesn't work well~~ ). Not sure if PCA is the right thing here as there is no difference between the 100th and the 200th RT. Each RT is independent and there is no specific ordering - so there is not any co-variance here that PCA can capture. (right ?) **(B)** Which classifier is recommended in such case ? Does anyone aware of such kind of work in general and specific for RT ? Any information on this would be appreciated.


title: Guide / tutorial for design & analysis of experiments in retail settings

text: I want to develop my skill set in the design and analysis of experiments in retail, where we conduct experiments (introduction of new package, new product, discount, etc.) in test stores while maintaining the status quo in selected control stores. After the end of the experiment, we conduct analysis of sales to see if there is a lift in the data. Does anyone know of any reading material about these concepts?


title: Choosing a regression model based on missing values

text: I'm trying to predict weight change with an intervention from baseline variables. Literature search yields suggests several predictors. Univariate analyses with weight change as dependent and baseline variables give several more predictors. I've entered both theoretical and GLM univariate predictors into a regression model using forced entry (in blocks according to significance for the GLM univariates). The problem I'm having is that my models vary hugely depending on how I treat missing variables. My study has 329 cases. I thought that listwise was the best method of exclusion (textbook) however this leaves me with insignificant models and only 196 cases. Pairwise and replacing values leave me with much higher cases naturally and give two models with the same significant predictors. My issue is that the predictors remaining are ones which make sense in context with the literature. Is it ok to present these? Could somebody please tell me why pairwise and replacing with means are frowned upon? I know there are other methods of missing value imputation but I have no idea of how to do that in SPSS, so if you're feeling generous and you know, I'd be grateful for any tips!


title: Forecast total for a year given monthly time series

text: I have a monthly time series (for 2009-2012 non-stationary, with seasonality). I can use ARIMA (or ETS) to obtain point and interval forecasts for each month of 2013, but I am interested in forecasting the total for the whole year, including prediction intervals. Is there an easy way in R to obtain interval forecasts for the total for 2013?


title: Strange outcome when performing nonlinear least squares fit to a power law

text: I have a data set (given below in my MATLAB code) y vs. x and my eventual goal is to fit it to a power law $y=ax^b$ to see what exponent $b$ I get. I did some non-linear least squares fitting and since NLLSF requires initial conditions I specified $b_0=10$ since that is what the expected value was. I noticed that after the fitting, $b=b_0=10$ was identical. It doesn't change at all. Then I did this for various values $b_0$ between five and fifteen and each time $b$ doesn't change. It stays identical to what its initial value was which is kind of annoying because we want to know what the exponent is. $a$ on the other hand seems to converge to the same values regardless of what $a_0$ is. Sorry if this gets too long but here I present everything I have done. ![](http://i.stack.imgur.com/tQ1vE.jpg) In this first figure, the top panel is the graph of the data. The second panel is the data (in blue) on a log-log scale to show how well it should fit to a power line. The dataset is definitely linear so I am not trying to fit something nonsensical. Just for giggles, I did a linear least squares fit in log-log space (shown in red) and I get the slope to be 7.568 and $r^2=0.998$ indicating an excellent fit. Then I started doing nonlinear least squares fits on the data. I fixed $a_0=1$ because it doesn't seem to matter. I get the same $a$ values no matter what. And then I varied $b_0$ between five and fifteen and here is the summary of the results. ![](http://i.stack.imgur.com/snC0s.jpg) The top panel shows what $a$ comes out to for various $b_0$. The second panel shows what $b$ converges to for various $b_0$. $b_0$ increments in steps of 0.1 and in every single case, $b$ doesn't change at all. It stays at whatever its initial value was. The third panel shows $r^2$ for various $b_0$ so at least the first half of the fits are good with a peak at around $b_0=8$. So my question is, what is going on here? What is so "special" about this dataset? This is too bizarre for me to have been made up by me. I only have those six data points available, no more, so I can't repeat the experiment nor add more data points. It seems to me that the residual space in $b$ is just very flat. Whatever $b_0$ you choose, the algorithm just stays there and doesn't move around but then for a variety of $b_0$ you get a nice fit as well. How can that be? How can all of the exponents between five and ten let's say ALL fit so well? Is there anyway for me to find out that exponent? Am I doomed to just doing a linear fit in log-log space and taking its slope as my exponent? Does it have something to do with the small number of (six) data points here? Like how trying to fit a seventh degree polynomial to six points will result in many polynomials being a good fit. Is this phenomenon common? Is it well- known/well-studied? Any nice references will be appreciated as well. Thanks everyone in advance. I am also pasting my code in MATLAB in case someone wants to see exactly what am I doing.               close all     clc     clear all          % The data hardcoded     x = 4:0.5:6.5;     y = 1.0e-04 * [0.011929047861902, 0.026053683604530, 0.057223759612162,...         0.117413370572612, 0.242357468772138, 0.462327165034928];          %The data plot     figure('units','normalized','position',[0 0 1 1])          % Shows the data as it is     subplot(2,1,1)     plot(x,y,'-o')     title('The Data y vs. x to be fitted')     xlabel('x');ylabel('y')          % Shows the data on a log log scale to see how well it should fit a power     % law     subplot(2,1,2)     loglog(x,y,'-o')     hold on     title('The Data y vs. x to be fitted on a log-log scale')     xlabel('x');ylabel('y')          % LLSF in log-log space to see how well it should fit a power law     ft = fittype('poly1');     opts = fitoptions(ft);     [fitresult, gof] = fit( log(x'), log(y'), ft)     loglog(x,exp(fitresult.p2)*x.^fitresult.p1,'-r')     hold off          % Will hold the results for various initial b values     fitresultarray=[];     binit = 5:0.1:15;          % NLLSF for various initial b-values are performed and saved     for counter = 1:length(binit)          ft = fittype( 'power1' );     opts = fitoptions( ft );     opts.Display = 'Off';     opts.Lower = [-Inf -Inf];     opts.StartPoint = [1 binit(counter)];     opts.Upper = [Inf Inf];          [fitresult, gof] = fit(x',y',ft,opts);          fitresultarray = [fitresultarray; fitresult.a fitresult.b gof.rsquare];          end          %fitresultarray          % The results of the fits are plotted     figure('units','normalized','position',[0 0 1 1])     subplot(3,1,1)     plot(binit,fitresultarray(:,1),'-o')     title('Nonlinear Least Squares Fitting y vs. x with Fitting Form y=ax^b')      xlabel('b initial');ylabel('a');     subplot(3,1,2)     plot(binit,fitresultarray(:,2),'-o')     xlabel('b initial');ylabel('b');     subplot(3,1,3)     plot(binit,fitresultarray(:,3),'-o')     xlabel('b initial');ylabel('Goodness of fit r^2');


title: Why it is natural to expect equality of variances?

text: When conducting various statistical test why do we expect equality of variances/ homoscedasticity/sphericity etc.?


title: What kind of attribute selection methods exist for high-dimensional regression data?

text: I have a 10,000 dimensioned dataset where all attributes are numeric values. I would like to select the best e.g. 50 attributes out of 10,000 so that I can run regression algorithms on it. I've tried Weka's PCA and CfsSubsetEval algorithms, however they were not capable of handling that much dimension (algorithms never terminated in both case). What kind of attribute selection algorithms (preferably with already implemented tools) exist in machine learning literature for such high dimensions? Summary of what I am doing is: I have a log data where all attributes are numeric (10K attributes in total and 100,000 samples in total) and I want to predict the bandwidth out of it. So I want to select best attributes which can play a role in predicting this bandwidth attribute.


title: Interpreting orthogonality

text: In a multiple partial linear regression setting, the book I'm reading has this sentence: «As a consequence of the fact that residual are orthogonal to explanatory variables, the 'cleaned' variables $M_2Y$ and $M_2X_1$(which are the residuals) are uncorrelated with $X_2$», where $X=(X_1 X_2)$, with $X_2$ being the last g columns of matrix $X$, and $M_2$ is a projection matrix to the orthogonal space of the space generated by the columns of $X_2$. I don't understand this sentence since $M_2Y$ and $M_2X_1$ are vectors of dim $n\times 1 $, and $X_2$ of dim $n\times g $. I understand the first part of the sentence, which means that $X_2^TM_2X_1=0$ and $X_2^TM_2Y=0$. The second part of the sentence I don't get it... Any help would be appreciated.


title: Comparing models

text: Suppose that I have two models, A and B (A nested in B), which are tailored to explain data from a single participant in an experiment. Example: I am modeling response times in a single participant. However, I have 20 participants. I can fit A and B using MLE to each participant independently, and then obtain likelihood ratios for each participant and test the significance with a chi- square test. But how can I make inferences over the population? I read in one (rather obscure and not totally trusted) source that I can simply sum the $\chi^2$ values for each participant. The resulting value will be $\chi^2$ with df=20 in this case.   1. Is this a valid procedure?   2. Can you recommend any sources that go into detail about this sort of method of aggregating model fits?  Bonus: What if my models aren't nested, but all of the above is still true?


title: Estimating demand elasticity econometrically

text: When specifying a production function for regression, it is well known that one of the features of using a log-log model is that the estimated coefficients are the output elasticities w.r.t. their respective independent variables. My question is does it then follow that if one regresses log(production) on log(price), the coefficient on log(price) will be the demand elasticity? That is, if we specify the production model as follows: $$\ln q=\alpha_0 + \alpha_p \ln p$$ (where q is quantity of output and p is output price) and then differentiate w.r.t. $\ln p$ $$\frac{d \ln q}{d \ln p}=\alpha_p$$ then isn't this the demand elasticity? If it is, then does the omission of other important variables in the production function bias the elasticity (if firms are not homogenous in those variables)? Thanks


title: Randomizing variables to get the best combination for high R-squared values in R

text: I've been spending quite some time to figure out how I can get the best R squared value from randomization of some values in a linear regression equation. I have allele frequency data and 14 environmental gradient data. Allele frequency value is fixed, but 2~14 combinations of the 14 environmental variables are used. My aim here is to find a combination of the environmental variables that yield high R squared value. Here is a simple linear regression equation code that returns R squared value.               > summary(lm(allele ~ compositevalues))$r.squared      "compositevalues" is a sum of standardized 14 different environmental values. I want to make 2~14 combinations of variables (with no replacement:i.e. var1+var2, var1+var3, var1+var4, var1+var2+var3, var2+var3+var4, var2+var3, var2+var4, var3+var4....etc. but not var1+var1+var2) as I mentioned above. I would appreciate it if you could instruct me on how to write a code that generate random combination of (sum of ) the variables and returns combinations of variables that are used with R squared value of >0.4. I was looking for permutation and resampling function in R, couldn't find ones that serve my purpose..... Below is a part of my data set.                1.  Location   allele           var1             var2          var3      2.  site1,     0.230271924,    -0.872093023,   -0.696403914,   -0.398671096      3.  site2,     -1.061563963,   0.944767442,    1.104640692,    -0.398671096      4.  site3,     -0.524508594,   0.339147287,    -1.296752116,   0.431893688      5.  site4,     0.027061785,    2.156007752,    -0.096055712,   0.431893688      6.  site5,     0.186726894,    0.944767442,    1.104640692,    -0.398671096      7.  site6,     -0.118088315,   -0.266472868,   -0.696403914,   -0.398671096      8.  site7,     -1.003503923,   0.339147287,    -1.296752116,   0.431893688      9.  site8,     -1.569589312,   0.339147287,    -1.296752116,   0.431893688      10. site9,     -1.119624003,   0.944767442,     0.50429249,    -1.22923588      11. site10,    1.362442702,    -1.477713178,   -0.096055712,   1.262458472      12. site11,    0.215756914,    0.339147287,    -1.897100318,   1.262458472      13. site12,    0.665722223,    -1.477713178,   -0.096055712,   1.262458472      14. site13,    1.086657513,    -1.477713178,   -0.096055712,   1.262458472      15. site14,    -0.001968235,   0.339147287,    1.704988894,    -2.059800664      16. site15,    -1.656679372,   0.339147287,    1.104640692,    -2.059800664      17. site16,    0.433482064,    0.339147287,    1.704988894,    -2.059800664      18. site17,    -0.814808794,   1.550387597,    -1.296752116,   -0.398671096      19. site18,    -0.713203724,   1.550387597,    -0.696403914,   -0.398671096      Many thanks!


title: A Reference for PAC-Bayesian?

text: I've recently came across topic known as PAC-Bayesain, but I cannot find a source to read about it. Any article that I came across are talking about its application in a specific area but there is no introduction to what it exactly is.


title: Explain the meaning of $\alpha$ in a Kolmogorov-Smirnov two-sample test

text: The Wikipedia article on the two-sample Kolmogorov-Smirnov test states that: > The Kolmogorov–Smirnov test may also be used to test whether two underlying > one-dimensional probability distributions differ. In this case, the > Kolmogorov–Smirnov statistic is > > $$D_{n,n'}=\sup_x |F_{1,n}(x)-F_{2,n'}(x)|$$ > > where $F_{1,n}$ and $F_{2,n'}$ are the empirical distribution functions of > the first and the second sample respectively. The null hypothesis is > rejected at level $\alpha$ if > > $$D_{n,n'}>c(\alpha)\sqrt{\frac{n + n'}{n n'}}.$$ It is not clear to me the meaning of the _$\alpha$ level_. Where does it come from and what does it mean statistically?


title: Analytical formula for distribution of partial sum of standardized random variable

text: I would like to know if there is an analytical formula for the distribution of partial sums of standardized random variables. (Of course, if one standardizes a random variable, the sum of _all_ the individual observations will be zero). I have written a Python function that estimates quantiles by bootstrapping. Each bootstrap iteration consists of three steps:   1. n samples of a N[0, 1] r.v. are generated.   2. The resulting sample is then standardized.   3. The sum of the first k elements of the standardized sample is calculated. This procedure is followed niter times, generating a sample of size niter. The quantiles of this sample are then calculated. In case it helps, here is the code:               import numpy as np     import scipy.stats.mstats          def boot_standardized(n, k, niter=1000, prob=None):         if prob is None:             prob = [0.05, 0.1, 0.5, 0.9, 0.95]              x = np.random.randn(n, niter)         x_std = (x - np.mean(x, axis=0))/np.std(x, axis=0)              boot_values = np.sum(x_std[:k,:], axis=0)              return scipy.stats.mstats.mquantiles(boot_values, prob=prob)


title: General location model

text: Let $Z_1$ and $Z_2$ be categorical random variables with $3$ and $2$ categories, respectively. Let $Y_1$ and $Y_2$ be $2$ continuous random variables. Define completely the GLOM (general location model) for the joint distribution of $Z=(Z_1,Z_2)^T$ and $Y=(Y_1,Y_2)^T$ I couldn't solve this problem. Can anyone help me?


title: Interpreting ACF and PACF for ARIMA model estimation for SAS

text: The following output is for a data set for monthly airport arrivals that have been log transformed, with first difference & seasonal difference. My interpretation from the SAS output is p=P=2 d=D=1 q=Q=1 How do you code the initial model estimation function by just interpreting the ACF PACF? diagnostics of the below show poor model selection or is my parameter interpretation completely off from the below chart?               proc arima data=temp;     identify var=logSTD(1,12) nlag=36;     estimate p=(2)(12) q=(1)(12) plot;     forecast lead=60 id=month  out=forecast; run;      ![enter image description here](http://i.stack.imgur.com/6CECB.png)


title: Transfer functions in R (TSA package)

text: In Time Series models’ transfer functions there is a decay parameter in the formula (let’s call it b). In TSA package that decay parameter is not mentioned. When I used other software before (such as SAS) I used to determine b after analyzing ‘prewhitened series’. But in TSA package in R there is no need to specify the decay parameter once you analyze CCF? If not how am I going to know when the decay starts? I understand CCF is used after prewhitening to determine how to filter the outputs but where b comes into the picture?


title: How to model conditional moments using an additive model

text: Let's say that I have the data set $(X,Y)$ where $X$ is a p-dimensional variable and $Y$ is uni dimensional. I'm interested in the following model: $$ \theta_y = E(Y|X) \\ \theta_{y^2} = E(Y^2|X) \\ \theta_{y^3} = E(Y^3|X) \\ \theta_{y^4} = E(Y^4|X) $$ where the components of $X$ enter the model additively. I'm interested on the left hand side for a posterior analysis. I'm barely aware of generalized additive models but I haven't used them before. What I need are pointers to theories and concepts that I should look to model and solve this problem.


title: Face Recognition using NMF

text: I am trying to implement a simple face recognition software using nonnegative matrix factorization. I have tried this on the ORL dataset which has 400 face images (40 people * 10 photographs). The input matrix consists of the column- reshaped-images ordered as V=[p1_1, p1_2, ..., p1_10, p2_1, p2_2, ...]. I expect to see that in the resultant encoding matrix, for the first 10 columns, the first row, and for the second 10 columns the second row must have bigger values, and so forth. But what I got is really bad. What is wrong here? Thanks...


title: Are discrimination parameters in two-parameter IRT models only specific to items?

text: Discrimination parameters in two-parameter model from IRT are usually considered item parameters. But I've come to doubt it. Think about psychophysics; for example, detecting luminance. I don't think anyone would say that discrimination (i.e., slope) parameter applies to luminance itself. It comes from human perception, so I figure it's a person parameter and it could be different for each person. What about test situation? Is it without any doubt that the discrimination parameter is considered as an item parameter? I doubt it is, especially given the fact that there is room for differential item functioning (DIF). So I'm concerned with the fact there's any theoretical basis for saying that discrimination parameters are item-parameters. I think I could model person-specific discrimination parameter model, in which item-discrimination parameter is determined by person-specific discrimination parameters in the group. That would justify DIF. So, my question would sum up like this:   1. Is that if there's any theoretical/empirical basis for saying that discrimination parameters are item-parameters? How can one be sure of it?   2. Is there any related reference for person-discrimination parameter model mentioned above?


title: Group items to perform correlations?

text: I have a question regarding SPSS. I have 7 constructs and each construct has 3 items/variables on a Likert scale. Example: Construct 1 - Mobile Usage, Item a - Liking, Item b - Preferences, Item c - Frequency I want to estimate correlations between these constructs. My questions are:   1. What's the best method to use for testing? Anova, Manova?   2. In order to estimate the correlation, do I need to group/combine the values of items a-c to represent Mobile Usage first and only then estimate the correlation or is there a way to do this in SPSS without grouping them? If so, which method accomplishes this? I really appreciate your help. Thanks very much.


title: Difference between absolute deviation to population median and sample mean

text: I have independent variables $X_i\in[0;1]$ and suppose they are uniformly distributed. If you want to minimize the total absolute deviation to a fixed number, how much can you gain from using the sample mean over the population median? Therefore I am looking for $f(n)=E_{\{x\}}\left(\sum^n|x_i-0.5|-\sum^n\left|x_i-\frac{\sum^n x_i}{n}\right|\right)$ Hope this is the correct notation :) The function depends on the number of variables used and the expectation is over all possible values of these $n$ variables. Apparently $f(1)=0.25$ but it seemed that $f(n>1)\approx 0.17$. It would be interesting to know the analytic expression for that. Any idea? PS: Are there easier or more interesting results if you use the squared deviation, or other distributions, or the mean of population quantity?


title: In an election, how can we tell the certainty that a candidate will be the winner?

text: There was a general election where I live yesterday and the television network started calling out winners long before all ballots were opened. They turned out right on all accounts, and I'm not really surprised they did. I know that statistics are absolutely viable. Still, I'm curious. Assuming:   * we have opened $i$ out of $j$ ballots;   * we have $n$ candidates whose current scores are $c_1, c_2, c_3, ... c_n$; How can we calculate the certainty with which the leading candidate is the winner?


title: Logistic regression, SPSS ignores my reference category and assumes another one

text: I am modelling logistic regressions in SPSS, the same model for different countries (well, with slight differences in the independent variables set due to collinearity diagnosis and stepwise results). The model seems to work fine for most countries. In two countries, I am having some issues with the same variable. The variable has three categories (hierarchical). In the output for those two beautiful countries, seems that SPSS take the reference category out, assumes other as reference category and gives exp(B) for the last category. Is there something wrong going on? What should / can I do about it? * * * Thanks a lot for your answers! The problem cannot be in the syntax since I am using the same for all countries and works perfectly. I am using ENTER. The STEPWISE was used as an exploratory method to identify the strongest predictors. According to the -2LL values my variable is not a 'strong' one for this two countries (0.12%). Is not the first time that I read that STEPWISE is a bad choice... I will that in account! Nevertheless, why do SPSS does that? Is this a real problem, meaning something that I must solve? Or can I use/report the odds for those two countries without that category adding a footnote? Should I run the model without that variable (theorectically the variable is not that important). Thanks ;) * * * Hi, All my variables are recoded to the reference category be 1 and most of my variables have three categories. I am using Enter with the following syntax. LOGISTIC REGRESSION VAR= "DV" /METHOD=ENTER "IV's" /CONTRAST ("IV")=Indicator (1) [...] /PRINT=GOODFIT CI(95) /CRITERIA PIN(.05) POUT(.10) ITERATE(20) CUT(.5). Do you have idea what is going wrong?


title: Are binomial and McNemar's tests right for my data?

text: I have performed an experiment, where I have real and synthesised movie data. My participants were shown either a real movie, or a synthesised movie (of talking heads), in a random order, and had to guess whether the movie was real or synthesised. I got a binomial test result of 0.58 (using MATLAB `binofit`). Of 2040 samples, 1189 were guessed correctly which I think shows significant difference. I also tried McNemar's test. The confusion matrix is:               [506 359]     [449 582]      where top left is 'said synth when synth', top right is 'said synth when real', bottom left is 'said real when synth' and bottom right is 'said real when real'. So anyway, McNemar's test for that gives me 9.91, which again I think shows that there is a significant difference. Are there any other tests I can try?


title: Three Parameter Gamma Distribution

text: Kindly help me with 1) The motivation for the three parameter gamma distribution and the resulting structure of its density. 2) The meaning of the location, scale and shape parameters here


title: Error scores analysis

text: The research paper examines the reaction times on a task and incorrect answers are eliminated as errors. The study does not specify the analysis however report a result of ps > 0.1 . What statistical procedure may they have used?


title: Bayes decision boundary of Figure 2.5 in Elements of Statistical Learning

text: When I read "Elements of Statistical Learning", I met some difficulty in calculating the Bayes decision boundary of Figure 2.5. In the package `ElemStatLearn`, it already calculated the probability at each point and used contours to draw the boundary. Can any one tell me how to calculate the probability? In a traditional Bayes decision problem, the mixture distributions are usually normal distributions, but in this example, it uses two steps to generate the samples, so I have some difficulty in calculating the distribution.


title: Robust correlation in R?

text: I would like to perform a robust correlation on a small sample (n<30). What is the best estimation method to use? I tried to get an overview over the plenty methods for robust statistics provided in R - I would be happy if anyone could give me some recommendations


title: Implication of uniform stochastic boundedness?

text: Let $\theta \in \Theta \subseteq \mathbb{R}^d$ be a parameter vector. Let $Q: \Theta \rightarrow \mathbb{R}$ be a function mapping from the parameter space to the real numbers. Let $Z_T$ be a a random vector representing a random sample, where $T$ is the sample size and $Z_T$ has support $\mathcal{Z}_T$. Let $Q_T(\theta): (\mathcal{Z}_T,\Theta)\rightarrow \mathbb{R}$ be an estimator of $Q(\theta)$. Assume that $$ \sup_{\theta \in \Theta}\sqrt{T}|Q(\theta)-Q_T(\theta)|=O_p(1) $$ i.e. $\forall \epsilon$, $\exists M_\epsilon>0$ s.t. $\mathbb{P}(\sup_{\theta \in \Theta}\sqrt{T}|Q(\theta)-Q_T(\theta)|>M_\epsilon)<\epsilon$ $\forall T$. Does this imply $$ \sqrt{T}|Q_T(\theta)-Q(\theta)|=O_p(1) \text{ } \forall \theta \in \Theta $$ ?


title: Bayesian Forecast with Minimal data

text: The following very simple forecast has been very helpful to me in applying basic methods.. and I have read that Bayesian methods may be superior for small data forecasting but I have not seen any useful comparative examples. Particularly ones that I could use to apply to the following problem. There are 9 years of January Unemployment data and Total unit sales of real estate in the corresponding area and year I have a projection for the 2013 January Unemployment rate From which I want to estimate the 2013 total unit volume. On the whole, in spite of the paucity of data in this model and the fairly wide confidence intervals the existing methods still seem somewhat useful to me,, While I keep reading the texts,, the greek and formulas make my head hurt.. and formulating simple problems seem to be largely overlooked in the texts I keep exploring. Hopefully the example we have plus one that someone :-) can provide showing a reasonable Bayesian approach would be useful to many (it will certainly help me!) Here is the r-arima "solution" arima regression Improving the model and just below is the initial data load for R, (complete code in the above link)                # load the base data as presented in the question      Workbook1 <- structure(list(Year = structure(1:10, .Label = c("31-Dec-04",      "31-Dec-05", "31-Dec-06", "31-Dec-07", "31-Dec-08", "31-Dec-09",      "31-Dec-10", "31-Dec-11", "31-Dec-12", "31-Dec-13"), class = "factor"),      total = c(51439L, 59674L, 58664L, 55698L, 42235L, 37918L,      36234L, 36965L, 39535L, NA), UnemplRt = c(5.7, 4.7, 3.8,      3.7, 4.3, 8.5, 10.9, 10, 8.3, 7.1)), .Names = c("Year", "total",      "UnemplRt"), class = "data.frame", row.names = c(NA, -10L))


title: Hausman test VS Mundlak model; the choice between fixed and random effects

text: I am using panel data and I have to choose between fixed-effects and random- effects models. I run the Hausman test, the H0 (i.e., the difference in the coefficients from the two models is not systematic) is rejected. Thus, I should use the fixed- effects. I also run a second test, which is based on the Mundlak model. So, I run the Mundlak regression (i.e., a random-effects model where I also add the individual means of the time-varying characteristics) and I run an F-test on the extra coefficients (the coefficients of the individual means of the time- varying characteristics). I cannot reject the null hypothesis. Thus, I should use the random-effect model (although also the fixed-effect is unbiased, the random effect is more efficient). I am puzzled. I test the same hypothesis in two different ways, obtaining two opposite results.


title: Sign of coefficient changes when contextual variables are added

text: I have an interpretation question. I am running binary multilevel models on whether or not households have bank accounts. Apart from relevant economic, social and demographical household level variables I also include two contextual variables for the primary sampling units (PSU), which are either villages or urban neighborhoods. These two measure the wealth level (assets) and the educational level of the PSU’s. Because a relevant explanatory variable is whether a bank branch is within the vicinity of the household and because households in urban areas should be expected to be closer to the nearest bank branch I expect the dummy variable for URBAN to be positively associated with having a bank account. The proportion of households with bank accounts are higher in urban (53.1) than in rural (30.9) PSU’s and average wealth and educational levels are higher in urban than in rural PSU’s. However, the dummy for URBAN has a positive sign when the two contextual variables for PSU wealth and education are not included and a negative sign when they are included. How should I interpret this? KML


title: Should I specify Prior or Cost matrix with Tree Bagger in Matlab

text: I'm trying to create Random Forests in Matlab and there are more observations in some classes than there are in others. Do I need to specify this as a cost matrix or as a prior probability or will Matlab figure this out automatically and the fact the data is skewed won't matter. Thanks


title: Random walk on simplex as part of Metropolis-Hastings

text: I would like to perform a random walk on a J-dimensional simplex. However, since this is part of a metropolis-hastings algorithm application, my understanding is that the steps need to be drawn from a symmetric distribution (is this correct?) I was wondering if there is a standard/established way to approach this. Any help/pointers are greatly appreciated! Noushin


title: Testing contrast between low/high and high/low

text: I'm trying to understand the results published in Sagristano, Trope, and Liberman (2002). Specifically, I don't understand what analysis produced the result highlighted in the screenshot below. The basic design is a 2 (time: near vs distant, between subjects) X 5 (probability of winning: .1, .3 ... .9, within subjects) X 4 (payoff: $4, $6, $8, $10, within subjects) ANOVA. They claim to have tested the hypothesis that > participants preferred high-probability, low-payoff gambles A colleague is trying to replicate this study, and we have not been able to figure out what they actually did here. Participants preferred high probability/low payoff to what? Perhaps the (unstated) second half of the hypothesis was "to low probability, high-payoff gambles", but if so I don't understand how they could have tested that. Perhaps the article is simply not clear, or perhaps I'm missing something obvious... The full article is linked above, and the relevant bit is reproduced here: ![enter image description here](http://i.stack.imgur.com/gLnnx.png)


title: Is it possible to combine bayesian SEM with PLS SEM?

text: I have already read some books about both two structural equation models. It seems both SEMs are suitable to the situation with small observations and large variables. I assume to use combine both two SEMs to create a more accurate model to analyse the data and structure.


title: Estimating if a distribution is significantly less conserved than others when one group is involved.

text: I am looking at levels of genes in a dataset and want to identify genes that do not vary much in terms of their expression level. While I can do this using the coefficient of correlation, calculating the covariance or by looking for number of genes within the botton x percent of genes using median absolute deviation those methods appear to be arbitrary. What I am interested in is defining a cutoff based on P.values - is there any way of finding out which genes show significantly less variability than would be expected by chance without having external controls to compare it to? Cheers, Ankur.


title: Negative Binomial Regression and Heteroskedasticity test

text: I am running a negative binomial regression in Stata and would like to know if I need to include the `vce(robust)` option in the model. I know the negative binomial regression is the best for the data I have. I've combed through so many sources and it seems like negative binomials already seem to take care of some heteroskedasticity, but I would really like a test that would let me know if there is still some in the model that needs to be taken care of via the `vce(robust)` option. I know there is the `hettest` but I can only use that in a `regress` model, not a `nbreg` (negative binomial regression) model. I also saw this advice: "You could try plotting the absolute value of Pearson residuals from your NB regression against covariates (at least, covariates you think might affect the overdispersion) or against the fitted means. If there's a discernible trend, this suggests non-constant variance." But what constitutes 'a discernible trend'? How do I know if that discernible trend would be significantly better taken care of by the `vce(robust)` option or it wouldn't matter?


title: Comparison between Two Groups

text: One of my demographic variable is age. Age is measured as continues data not categorical. If I want to test differences between two groups to determine whether there is significant difference between them or not, shall I use independent t-test, or I have to covert it to categorical variable then calculate chi-square. Thank you


title: If original data are available, should replication data be added and everything re-analysed (IPD) or is a meta-analysis better?

text: One common practice which increases the odds of obtaining spurious results is to keep collecting observations after a preliminary analyses are performed[1]. This occurs when the cutoff point for collecting cases is set as the time when significance is reached. Noting that, it seems that meta analyses are a valid way of combining a variety of studies in a way that minimises overall error in the parameter estimate. But I guess that it's also a result of not having access to the original data. **_So my question is_** , what if you _do_ have access to the original data from previous experiments? Is it better to do a meta analysis of the results or to combine the data from the previous studies and run the analyses again, apparently called individual participant data (IPD),? My guess is that the latter is ill-advised although some people do it. For one, it appears to be suspiciously close to the undefined cutoff point discussed in Simmons et al. It also seems to go against other advice I have received on this site regarding dipping into the same data set multiple times. Note that in my particular case, I may have access to _summary_ data of each participant in a study I replicated, not all data points collected for every participant. Essentially I may have access to means, totals, etc for each participant in a previous study. * * * [1]: see, for example, Simmons, J. P., Nelson, L. D., & Simonsohn, U. (2011). False-Positive Psychology: Undisclosed Flexibility in Data Collection and Analysis Allows Presenting Anything as Significant. _Psychological Science_ , _22_ (11), 1359–1366. doi:10.1177/0956797611417632


title: Multiple Imputation Variance

text: I have been reading up on multiple imputation, and I am interested in the between-imputation variance. However, not in the estimation of the parameters, but in the imputed values themselves. (1) Suppose $x$ is vector with some missing values. I then perform multiple imputations an infinite amount of times according to Rubin by drawing from the conditional distributions given the observed data and then take the between- imputation variance of the completed vectors (2) Also suppose I also perform multiple imputations by drawing from the joint probability distribution an infinite amount of times. Then I take the between- imputation variance. Under the MAR assumption, is the variance of (1)<(2)? If the observed values help predict the missing values, it seems like it is intuitively...


title: Plotting a model in R

text: I have used the _vglm()_ function from the VGAM package to obtain a Weibull model for some data. I would like to plot the curve that this model is represented by. How do I do this? Here's my input and output (I am comparing 'p_seen' to 'calculated_logmar', which exist in the dataframe 'dframe1'):               > my_model <- vglm(formula = p_seen ~ calculated_logmar, family = weibull, data = dframe1)          Pearson Residuals:       log(shape) log(scale)     1    0.78928  -0.093136     2   -1.21092  -0.978527     3    0.12519   0.020018          Coefficients:                          Estimate Std. Error   z value     (Intercept):1      1.4797e+01 1.0803e+00  13.69619     (Intercept):2      4.5824e-06 2.3095e-05   0.19841     calculated_logmar -1.5324e+01 1.4448e+00 -10.60582          Number of linear predictors:  2           Names of linear predictors: log(shape), log(scale)          Dispersion Parameter for weibull family:   1          Log-likelihood: 9.46132 on 3 degrees of freedom          Number of iterations: 30       I understand that the coefficients must relate somehow to the model itself, but I have no idea how to turn these coefficients into meaningful plots. PS: I would like to use _ggplot2_ to visualise my data if possible. **Edit** : I have tried using the command:               plotvgam(my_model, se = TRUE, lcol = dframe1$calculated_logmar[1], scol = dframe1$calculated_logmar[1], which.term = 1)      This gives me a plot of 'calculated_logmar' against 'partial for calculated_logmar'. What I'm looking for is the actual curve of the model, showing calculated_logmar plotted against p_seen. _If_ it is possible, I'd like to render this using _ggplot2_. Any help would be appreciated.


title: Calculating confidence intervals for a proportion when there are no 'successes' in the sample

text: I'm looking to analyse 400k replies to a Facebook-equivalent post to determine how many of them are written by bots, and how many of them are written by real people. I don't have the resources to obtain all 400k replies, so I can only work with a sample (say, 4000 replies). I believe it's very likely that virtually all, but definitely not all of the replies are written by bots, but I need to prove it with some numbers. From my reading, I understand that usually what you would do would be to:   * take, say, 40 chunks of 100 replies as your sample,    * for each sample, determine what percentage is written by bots   * calculate the mean and standard deviation of the percentages of all 40 chunks   * calculate the standard error for the sample, based on the standard deviation / sqrt of 40 (chunks)   * calculate 95% significance of the mean percentage via mean of sample chunks +/- 1.96 x standard error To what I understand, this can be done because we expect the data to be roughly normal. With skewed data, you need more samples (chunks, in this case). However, I'm expecting that the % of bots for each chunk will be 100% each time for my 40 samples (let's assume my expectation is right). So what I'm finding is that the sample results are uniform --> Standard deviation = 0 --> SE = 0 --> 99.999% confidence interval = 100%, 100%, 100%. How do you infer on a population when all the sample means are uniform (i.e., all 100%), other than to conclude that the population's mean is also exactly the same (i.e., 100%). But I've only reviewed 400 / 400k replies. I can't really say that I am 99.999% confident that the % of bots for the population is 100% right?


title: Can I use xreg with stl decomposition to handle moving holiday?

text: I am trying to decompose and forecast a weekly time series which is believed to be affected by moving holidays (e.g. Chinese New Year, which happens in different weeks of a year). I would like create a regressor variable to reflect the holiday effect on the series. Is it correct to use the regressor variable as xreg in forecasting stl object / stlf? Also, I would like to know the difference between the following methods, and whether they are doing the job I wanted. 1) decompose using stl, then forecast the decomposed object, i.e.                  model<-stl(tseries,"periodic")        forecast<-forecast(model,h=10,method="arima",xreg=xreg,newxreg=newxreg)      2) use stlf directly, i.e.               forecast<-stlf(tseries,h=10,method="arima",xreg=xreg,newxreg=newxreg)      Thanks in advance


title: PCA + L1 equivalent to elastic net?

text: I am performing a logistic regression on a rather big dataset (700k+ samples and 1k+ features). I suspect that a lot of these features will be highly correlated and multicollinearity can be an issue. I believe that with elastic net regularisation is can perform feature selection and shrink the parameters of correlated features toghether (to prevent overfitting) I understand that after PCA I end up with new (linnearly) uncorrelated (orthonal) features. I'm I correct to think that I now only need an L1 regularisation term in my optimazation function to have a simmular model as in the case with elastic net? I prefer to do the PCA because I can reduce the number of parameters from 1k+ to around 60. (the rest of the PCs have 0 variance) which is more manageble to work with. I am only interesed in a developing a prediction algorithm. I use R. Best regards


title: Differences in AUC calculation between pROC and ROCR

text: Does anyone know the difference in calculation between these two AUC packages? They get different results when I add in positives with predicted value of 0 (simulating a prob model where many outputs will be zero e.g. randomForest). So I would like to know why there is the kink in the pROC results (red) and the AUC 'flips'. The plot below illustrates this behaviour. Also, this stems from the problem of having to compute AUC's on models with many positives assigned probabilities of zero.                library(ROCR)      library(pROC)            start.ds <- data.frame(response=gl(2, 1000, labels = c("1", "0")),pred=ifelse(as.numeric(gl(2, 1000, labels = c("0", "1")))==2,0,1))       add.positives <- data.frame(response=gl(1, 1000, labels = c("1")),pred=0)       add.positives <- rbind(add.positives, add.positives)            nadd <- nrow(add.positives)       theSeq <- seq(1, nadd, 10)       allAUCs <- NULL       for(i in theSeq){            print(i)            temp.ds <- rbind(start.ds,add.positives[1:i,])            temp.ds <- temp.ds[order(temp.ds$pred,decreasing=TRUE),]                 aucPROC <- as.numeric(auc(temp.ds[,1],temp.ds[,2]))            aucAUCPackage <- performance(prediction(temp.ds[,2],temp.ds[,1]),"auc")@y.values[[1]]            allAUCs <- rbind(allAUCs, data.frame(obs=i, aucPROC=aucPROC, aucAUCPackage=aucAUCPackage))              }            plot(allAUCs$aucPROC,col="red",type="l")       lines(allAUCs$aucAUCPackage,col="black",type="l")      ![enter image description here](http://i.stack.imgur.com/IPjU8.png)


title: Logistic regression: Fisher's scoring iterations do not match the selected iterations in glm

text: it happened to me that in a logistic regression in R with `glm` the Fisher scoring iterations in the output are less than the iterations selected with the argument `control=glm.control(maxit=25)` in `glm` itself. I see this as the effect of divergence in the iteratively reweighted least squares algorithm behind `glm`. My question is: under which criteria does `glm` stop the iterations and provides with a partial output? I was thinking about something like "when the new coefficients-old coefficients < epsilon, then STOP". Is this the case? If not, what does make `glm` stop? Thanks, Avitus


title: List of suitable job positions for data analyst

text: While searching for a job for data analyst I found out that many jobs where data analysis is widely used do not include "data" or "statistics" in position title, for example credit specialist, business analyst and so on. What are suitable job positions for people who want to analyze data? Moreover, how they differ?


title: What is the optimal ratio of sample size to the number of parameters in a multiple regression?

text: **What is be the optimal ratio of sample size to the number of parameters in the multiple regression model?** I am wondering as I would like to improve accuracy of prediction. Some sources suggest a ratio of 3 to 1, whereas others suggest 10 to 1. Any other suggestions?


title: Distribution based on rate of distinct values

text: I am investigating the following problem: I have a large set of values, many of which are repeated. Measuring the number of distinct (or unique) values shows that their number grows much slower than the number of total values, so that with 2M values I get just ~260k distinct values, with 20M values it goes up to just ~750k distinct, and after processing 40M values there are just under a million distinct values. ![Graph of distinct values vs values](http://i.stack.imgur.com/54e5V.png) The values are real-world data points created by joining together a number of variables, each of which may take a (potentially large but finite) set of values. These variables are not necessarily independent; on the contrary, they are most likely dependent on each other. Intuitively the values must follow some kind of power law distribution, but I cannot fit the values to a satisfying formula. I can tell you it is not an exponential curve because it isn't straight on a log-log plot; it does not appear to be logarithmic either. The questions that I would like to answer are:   * How many distinct values should I see after processing a given number of values?   * What is the distribution `f(x)` of distinct values? So that `f(1) = 1` (since the first value is certainly distinct) and `f(2e6) ~ 261620`.   * What are the parameters of the distribution? So I can measure them after changing (or adding or removing) some of the independent variables and see how it affects the number of distinct values.   * Related to the above, how does each variable contribute to the number of distinct values? I can experiment numerically, but I would like to be able to characterize the results using numerical parameters. This would tell me how dependent variables are.   * For bonus points: how can I find out how many bits of information each variable carries? This will be related to how dependent it is. I have not found on the Wikipedia any info about distinct values; in the literature there are some articles such as this one but they don't delve on the number of distinct values. Apparently this is called a "birthday problem". Any help is appreciated. Note: initially asked on math.stackexchange.com, where it has been closed. The actual numbers follow (values, distinct values):               1, 1     2000000, 261620     4000000, 367883     6000000, 446604     8000000, 511565     10000000, 565090     12000000, 610213     14000000, 649240     16000000, 684861     18000000, 717790     20000000, 746938     25000000, 810937     30000000, 865702     35000000, 915417     40000000, 966590


title: Least squares robust to outliers

text: I have sparse overdefined system of linear equations. For example I have n variables, m equations(m>n) and k equations from m are "bad" equations that represent outliers. Is there any methods to solve this problem? I already found techniques such as iterative reweighred least squares,LMeds,M-estimator but I'm not sure which is my case? Here is my task in detail: I want to solve computer vision problem which seems to be called "panorama global registration" like described here (they use weighted least squares). Main problem is that when I calculate pairwise relations between images (which is then become equations in my system) even if I prune weak connections by threshold I have some bogus connections, so I need robust algorithm that can handle outliers.


title: Coefficients paths – comparison of ridge, lasso and elastic net regression

text: I would like to compare models selected with ridge, lasso and elastic net. Fig. below shows coefficients paths using all 3 methods: ridge (Fig A, alpha=0), lasso (Fig B; alpha=1) and elastic net (Fig C; alpha=0.5). The optimal solution depends on the selected value of lambda, which is chosen based on cross validation. ![Profiles of coefficients for ridge \(A, alpha=0\), lasso \(B, alpha=1\) and elastic net \(C, alpha=0.5\) regression. Numbers at the top of the plot represent the size of the models.The optimal solution depends on the selected value of lambda. Selection of lambda is based on cross validation. ](http://i.stack.imgur.com/e8dUs.jpg) When looking at these plots, I would expect the elastic net (Fig C) to exhibit a grouping effect. However it is not clear in the presented case. The coefficients path for lasso and elastic net are very similar. What could be the reason for this ? Is it just a coding mistake ? I used the following code in R:               library(glmnet)     X<- as.matrix(mydata[,2:22])     Y<- mydata[,23]     par(mfrow=c(1,3))     ans1<-cv.glmnet(X, Y, alpha=0) # ridge     plot(ans1$glmnet.fit, "lambda", label=FALSE)         text (6, 0.4, "A", cex=1.8, font=1)         ans2<-cv.glmnet(X, Y, alpha=1) # lasso         plot(ans2$glmnet.fit, "lambda", label=FALSE)     text (-0.8, 0.48, "B", cex=1.8, font=1)     ans3<-cv.glmnet(X, Y, alpha=0.5) # elastic net      plot(ans3$glmnet.fit, "lambda", label=FALSE)     text (0, 0.62, "C", cex=1.8, font=1)      The code used to plot elastic net coefficients paths is exactly the same as for ridge and lasso. The only difference is in the value of alpha. Alpha parameter for elastic net regression was selected based on the lowest MSE (mean squared error) for corresponding lambda values. Thank you for your help !