title: Producing class labels from sliding window

text: I have got a dataset: > The attribute of interest: Count. I work out the linear model of Count, then > the root mean square error (`gmdl.rmse`). > > `buffer1` and `buffer2` of length `w`. In this case `w=3`. Therefore > `buffer1` is the the first three elements [1:3] of Count, and `buffer2` is > [4:6]. I work out the root mean square error of the buffers, then `abs(b2.rmse-b1.rmse) >= gmdl.rmse`   * If FALSE, buffer1 and buffer2 slide by 1 position (`b1[2:4], b2[5:7]`) and then I work out the rmse, and repeat the above code. If FALSE again slide by 1 again, repeat process until TRUE.   * If TRUE, say for `b1[2:4],b2[5:7]`, then a change point is recored (C1), then `b1=b1+w,b2=b2+w` (i.e. `b1[5:7],b2[10:12]`). For the new b1 and b2, if TRUE again, then another change point is recored (C2), and b1 and b2 slide by w (`b1[10:12],b2[13:15])`; say TRUE again, then C3 recored, then `b1[13:15],b2[16:18]`, and TRUE again then C4 recorded and b1, b2 slide by w. (100 change points in total). Define, > A1=the average of elements in the interval C1:C2, > > A2=the average of elements in the interval C3:C4, If `A1 > A2`, elements in the interval C1:C2 are labelled as "HIGH", in an array called `class`. I've got all the change points, but am unsure how to code it (in R) for the `A1 >A 2` and labeling stuff. For labeling a while loop with an if statement is all I need, apparently. so I assume it goes something like this:               class<-array("Low",c(700))     A<-array(0,c(50))     index1<-1     Low<-1     High<-2     while ((High)<length(Count))     {            A[index1]<-mean(Count[CP[Low]:CP[High]])         Low<-Low+2         High<-High+2         index1<-index1+1         i<-2         ii<-4         if (A[i]>A[(i-1)])         {             class[CP[(ii-1)]:CP[ii]]<-"High"             i<-i+2             ii<-ii+4         }     }      where `CP` is a vector of all the change point positions. After a lot of trial and error, the above code comes close to what I want. Any ideas as to how to get the desired results? Thanks.


title: How to replace the top 1% of the data with a new value in R?

text: In `R`, given a set of numerical data, how do I decide the top 1% of these data, and replace them with a new value, say, e.g. $10^{10}$


title: Cohen d for 2x2 interaction

text: There is a nice answer to this question, but it assumes that you have the ANOVA table available. My problem is different. Say I'm reading a paper describing a male vs female, disease vs control experiment, and I know the n's, means, and standard deviations for all four groups (healthy females, ... , diseased males), but I don't have the original data. I estimate Cohen d for the sex:health interaction. Numerator: (male_disease- female_disease)-(male_healthy-female_healthy). Denominator: I pool the standard deviations. In simulations with perfect data (normally distributed, etc), my estimated Cohen d does not match the Cohen d calculated from eta-squared: $$2\sqrt\frac{\eta^2}{1-\eta^2}$$ I am sure there must be a closed-form estimator for Cohen's d for this case, but having looked all over google and online stats books I can't seem to find the answer. Apologies if I've missed something obvious.


title: Issues in pattern recognition and plots

text: I am trying to compare the performance of a classification result with Bayes classifier, K-NN, principal component analysis (PCA). I have doubts regarding the following (please excuse my lack of programming skills since I am a biologist and not a programmer thus finding hard to follow the matlab documentation)   1. In the Matlab code               Class = knnclassify(Sample, Training, Group, k)     Group =  [1;2;3]   //where 1,2,3 represents Class A,B,C respectively.      What goes in the sample because my data is a 100 row 1 column for each of the classes? So Group 1 contains data like $[0.9;0.1;......n]$ where $n=100$. Would the sample be a vector containing random mixtures of the data points from the three classes? Same query/doubt for Training matrix.   2. I cannot follow `crossval()` & `cvpartition()` function given in MATLAB documentation crossval(). Would be obliged if a simpler version of it is provided here.   3. Is there a code for multi-class ROC since Internet resources provide information about binary class ROC. How to proceed for ROC plot when there are three or more classes like in my case?   4. Is it possible to obtain a surface plot from the confusion mmatrix from the command `confusemat()`?


title: Fixed point iterations for expectation propagation using energy minimization

text: **EP Primal** In 1, it is finding the EP iterations by solving a saddle-point problem on the energy function. First, the primal is claimed to be $$ \min_{\hat{p}_i} \max_{q} \left[ \sum_i \int_{\mathbf{y}} \hat{p}_i(\mathbf{y}) \log \frac{ \hat{p}_i(\mathbf{y}) }{ t_i(\mathbf{y}) p(\mathbf{y}) } d\mathbf{y} - (n-1) \int_\mathbf{y} q_\theta(\mathbf{y}) \log \frac{q_\theta(\mathbf{y})}{p(\mathbf{y})} d\mathbf{y} \right] $$ with the local moment matching constraints $$ \mathbb{E}_{ q_\theta(\mathbf{y}) }\left[ \phi(\mathbf{y}) \right] = \mathbb{E}_{ \hat{p}_i(\mathbf{y} ) }\left[ \phi(\mathbf{y}) \right], \forall i \quad \quad $$ **EP Dual** The dual energy function is the following; $$ \min_{\nu} \max_{\lambda} \left[ (n-1) \log \int_\mathbf{y} p(\mathbf{y}) \exp \left( \nu^\top \phi(\mathbf{y}) \right) d\mathbf{y} - \sum_{i=1}^{n} \log \int_\mathbf{y} \hat{t}_i(\mathbf{y}) p(\mathbf{y}) \exp \left( {\lambda_i}^\top \phi (\mathbf{y}) \right) d\mathbf{y} \right], $$ $$ (n-1) \nu = \sum_i \lambda_i. $$ **EP fixed point iterations** And using the dual energy function, we should be able to find the fixed point iterations: _Message elimination_ : Choose a $\tilde{t}_i$ to do approximation with.  Remove the factor $\tilde{t}_i$ from approximation, $\; q_\theta^{-i} = \displaystyle \frac{q_\theta}{ \tilde{t}_i }$  _Belief projection:_ Project the approximate posterior, with $\tilde{t}_i$ replaced with $t_i$, on the approximating family, $$ q^{new}_\theta(\mathbf{y}) = \text{proj}\left( \hat{p}_i(\mathbf{y}) \rightarrow q_\theta(\mathbf{y}) \right), $$ where, $$ \hat{p}_i(\mathbf{y}) = \frac{1}{Z} q_\theta^{-i}(\mathbf{y}) t_i(\mathbf{y}), \; \; Z = \int q_\theta^{-i}(\mathbf{y}) \times t_i(\mathbf{y}) d\mathbf{y} $$ $ \tilde{t}_i = \arg \min_{\tilde{t}_i} \text{KL} \left( \displaystyle \frac{ t_i \prod_{j\neq i} \tilde{t}_j }{ \int t_i \prod_{j \neq i} \tilde{t}_j } \parallel q_\theta(\mathbf{y}) \right) $.  Message update: Compute the new approximating factor, $$ \tilde{t}_i = Z \frac{ q^{new}_\theta(\mathbf{y}) }{ q_\theta^{-i}(\mathbf{y}) } $$ Here are the questions:   1. I know how how to derive dual from primal, but it it not clear to me where the primal is coming from.    2. I don't see how can I find the EP iterations from Dual. Any idea?


title: Does Julia have any hope of sticking in the statistical community?

text: I recently read a post from R-Bloggers, that linked to this blog post from John Myles White about a new language called Julia. Julia takes advantage of a just-in-time compiler that gives it wicked fast run times and puts it on the same order of magnitude of speed as C/C++ (the same _order_ , not equally fast). Furthermore, it uses the orthodox looping mechanisms that those of us who started programming on traditional languages are familiar with, instead of R's apply statements and vector operations. R is not going away by any means, even with such awesome timings from Julia. It has extensive support in industry, and numerous wonderful packages to do just about anything. My interests are Bayesian in nature, where vectorizing is often not possible. Certainly serial tasks must be done using loops and involve heavy computation at each iteration. R can be very slow at these serial looping tasks, and C/++ is not a walk in the park to write. Julia seems like a great alternative to writing in C/++, but it's in its infancy, and lacks a lot of the functionality I love about R. It would only make sense to learn Julia as a computational statistics workbench if it garners enough support from the statistics community and people start writing useful packages for it. My questions follow:   1. What features does Julia need to have in order to have the allure that made R the de facto language of statistics?   2. What are the advantages and disadvantages of learning Julia to do computationally-heavy tasks, versus learning a low-level language like C/++?


title: Software packages that can construct feature representation for a given text file using N-gram model

text: Are there any open source software packages, (Java, Matlab, R) that can generate a feature representation for a text document using N-gram model?


title: How to generate a group summary categorical variable for survival analysis in Stata?

text: I have registry data for which I have data for individuals for each hospital attendance over a 10 year period. I wish to determine the effect of chronic antibiotic use upon acquisition of specific infections (outcome). However, often the choice of antibiotic may differ over time and so I want to be able to summarise individuals' antibiotic use so that eventually I can group individuals with similar exposures. In order to do this I have generated an indicator variable that records each hospital visit as first, second visit, etc. However I cannot work out how I may then use this to summarise individuals' exposure. Any ideas?


title: PLS in R with the pls package

text: I am very new in PLS and I try to understand the output of the R function plsr(). Let us simulate data and run the PLS:               library(pls)     n <- 50     x1 <- rnorm(n); xx1 <- scale(x1)      x2 <- rnorm(n); xx2 <- scale(x2)     y <- x1 + x2 + rnorm(n,0,0.1); yy <- scale(y)     p <- plsr(yy ~ xx1+xx2, ncomp=1)      I was expecting that the following numbers $a$ and $b$               > ( w <- loading.weights(p) )          Loadings:         Comp 1     xx1 0.723      xx2 0.690                          Comp 1     SS loadings       1.0     Proportion Var    0.5     > a <- w["xx1",]     > b <- w["xx2",]     > a^2+b^2     [1] 1      are calculated in order to maximize               > cor(y, a*xx1+b*xx2)               [,1]     [1,] 0.9981291      but this is not exactly the case:               > f <- function(ab){     + a <- ab[1]; b <- ab[2]     + cor(y, a*xx1+b*xx2)     + }     > optim(c(0.7,0.6), f, control=list(fnscale=-1))     $par     [1] 0.7128259 0.6672870          $value     [1] 0.9981618      Is it a numerical error, or do I misunderstand the nature of $a$ and $b$ ? I would also like to know what are these coefficients:               > p$coef     , , 1 comps                     yy     xx1 0.6672848     xx2 0.6368604       **EDIT** : Now I see what is p$coef :               > x <- a*xx1+b*xx2     > coef(lm(yy~0+x))             x      0.9224208      > coef(lm(yy~0+x))*a             x      0.6672848      > coef(lm(yy~0+x))*b             x      0.6368604       So I think I'm right about the nature of $a$ and $b$. EDIT: In view of the comments given by chl I feel my question is not clear enough, so let me provide more details. In my example there is a vector $Y$ of responses and a two-columns matrix $X$ of predictors and I use the normalized version $\tilde Y$ of $Y$ and the normalized version $\tilde X$ of $X$ (centered and divided by standard deviations). The definition of the first PLS component $t_1$ is $t_1 = a \tilde X_1 + b \tilde X_2$ with $a$ and $b$ chosen in order to have a maximal value of the inner product $\langle t_1, \tilde Y \rangle$. Hence it is equivalent to maximize the correlation between $t_1$ and $Y$, isn't it ?


title: Testing strict exogeneity in time series

text: One of the important OLS assumptions is a strict exogeneity assumption, i.e. $E(\epsilon_i | X) = 0, \forall i$. I'm interested in testing empirically this hypothesis, notably in the context of time series. It is known that exogeneity is rarely true in time series but assuming that the model is well specified, e.g. contains no lagged dependent variables, is there a formal test (in R) that could indicate there's a problem? What is the usual approach to testing this assumption?


title: k-means in R generates same number of clusters but different cluster label

text: Hi there I am running a k-means code in R with the same data and with the same number of clusters, in this case 3, but each time that I run the code, the cluster label changes, for example. In the first run I got this results Col A Cluster.K-means FirstOBS Cluster1, SecondOBS Cluster1, ThirdOBS Cluster2, FourthOBS Cluster3 then I run the code again, specifying the same number of clusters=3, and the output column changes the cluster number, for example FirstOBS Cluster2, SecondOBS Cluster2, ThirdOBS Cluster3, FourthOBS Cluster1 In this case the second run (with the same data, same code and same number of clusters) changed the (1 by 2), the (2 by 3) and the (3 by 1) How can I fix the result? because in the next step I have to compare the cluster number between two different data, for instance the cluster results from January against March results. This is the code `cv1.km3 <\- kmeans(cv1.f, 3) cv1.km3$size` So each time I run the cv1.km3$size I got different combinations Thank you


title: How can conditional be substituted with a function of same variable?

text: I was studying on a PAMI article and I have seen an equation like this: $$ \begin{align*} P(B|X) &= \prod_c{P(B^c|X)} \\ &= \prod_c{P(B^c|A^c)} \end{align*} $$ where $A^c$ is a function of $X$. I know to fully understand the derivation details are needed, but what I am wondering is in under what circumstances we can substitute the conditional with a function of the variable? It does not seem nice to me, because the result of the function may not carry all the information. The details, motivation, explanation are not given regarding to this substitution in the paper, so I am stuck. The rest of the paper is well written except a few unimportant typos, so I am assuming this is something that I don't know. Edit: $c$ is not for complement. $B$ is a vector where every element of is independent of each other. $A$ is a vector with same dimensions as $B$. $c$ is an index for the vector's element. My main area of work is not statistics, sorry for any unclear points.


title: PLS regression is not working on weighted data

text: I was running PLS regression on the data which is weighted and gives the following error message:               Error: The dataset is weighted, but this procedure cannot be used with weights      I do not understand the reason for this as earlier with weights on, there was an outcome.


title: Dividing an estimate with confidence limits

text: I have an estimated population size for which individual covariates are known for only a smaller sample of the population. So for example, my total population size is estimated to be 26.9 individuals/ha, SE = 4.4, lcl = 19.6, ucl = 36.9. I know that the population is 30% male, and 70% female (from a smaller sample of the population). Is there a way to estimate the total number of males and females with confidence limits? Any suggestions would be greatly appreciated.


title: Why use a likelihood ratio (and its relation with p-value)?

text: Lets' say I have a Poisson distribution for which I use a _maximum likelihood_ defined as: $P_i=\frac{m_i^{n_i}}{e^{m_i}n_i!}$ where $m_i$ represents the **model** value of bin $i$ (real, $m_i$> 0) and $n_i$ its **observed** value (integer, number of counts in bin $i$). The cumulative likelihood for the whole data set is then: $L=\prod\limits^{N}_{i=1} \frac{m_i^{n_i}}{e^{m_i}n_i!}$ whwre $N$ is the total number of bins. Now, as I understand it, $P_i$ denotes the probability that the observation $n$ is drawn from model $m$. (1) So a _bigger likelihood_ value $L$ means a _better fit_ between the model and the data (Is this correct?) What I don't get is the necessity for a _maximum likelihood **ratio**_ , defined as (2) the ratio of probability of drawing $n_i$ points from model $m_i$, from that of drawing $n_i$ points from model $n_i$ (Is this correct?) whose cumulative form is: $LR= \prod\limits^{N}_{i=1} \frac{\frac{m_i^{n_i}}{e^{m_i}n_i!}}{\frac{n_i^{n_i}}{e^{n_i}n_i!}} = \prod\limits^{N}_{i=1} \left(\frac{m_i}{n_i}\right)^{n_i}e^{n_i-m_i}$ (3) Won't maximizing this value $LR$ give me the same results as maximizing $L$? Finally: (4) How can I calculate the **p-value** associated with a _Poisson likelihood_ $L$ or a _Poisson likelihood ratio_ $LR$ for a model($m$)-observation($n$) analysis?


title: Covariance of linear regression coefficients in weighted least squares method

text: I'm a Physics student and we had a lab where we had to process some data in the form $x_i$, $y_i$, $\sigma_i$, where the $\sigma_i$'s are the uncertainties of the $y_i$'s. With the weighted least squares method we obtained coefficients $A$ and $B$, the estimated parameters for a straight line in the form $y = A+Bx$ that would fit the data. In particular, $A=(s_x s_{xx} - s_x s_{xy})/D$, $B=(s_1 s_{xy} - s_x s_y)/D$, where $D=s_1 s_{xx} - s_x s_x$, $s_1 = \sum_i \frac{1}{\sigma_i^2}$, $s_{xx} = \sum_i \frac{x_i^2}{\sigma_i^2}$, $s_{xy} = \sum_i \frac{x_i y_i}{\sigma_i^2}$, $s_x = \sum_i \frac{x_i}{\sigma_i^2}$, $s_y = \sum_i \frac{y_i}{\sigma_i^2}$. In addition we have named $\delta_{AB}$ the variable $-s_x/D$. In light of the use of $\delta_{AB}$ in the subsequent data analysis, it seems that this variable represents the covariance between $A$ and $B$. In my statistics courses we never covered the covariance between linear regression coefficients, so my question is what is the general theory behind $\delta_{AB}$? Is it possible to prove simply, in the simple case of two variables, that this is indeed the covariance between $A$ and $B$?


title: Testing for trends in partial proportional odds models

text: I am trying to build a model to explain an ordinal response variable $y$ with 4 levels: $y_0$, $y_1$, $y_2$ and $y_3$. The independent variable in this model is $v$. $v$ is a categorical variable with three categories $v_a$, $v_b$ and $v_c$. The assumption of proportional odds fails for $v$. Therefore, I am using the partial proportional odds (PPO) model. The model gives me 3 ORs for $y_1$, $y_2$ and $y_3$ that correspond to the odds for events $\{y>=y_1\}$, $\{y>=y_2\}$, and $\{y>=y_3\}$. I would like to perform a test for trend for each of these 3 events such that each test gives a p-value. For example, the test for trend for $y_1$ would explain if there is a trend between the change in level of v and in the odds for $\{y>=y_1\}$. Is there a way to use the partial PO model directly to compute a trend test for the 3 events? Here’s a solution I am thinking of:   1. I first convert the categorical variable, v into a continuous variable with values 1,2, and 3 and build a PPO model using this continuous variable - let’s call this model PPOC.   2. For each outcome level $i$, I divide the samples in the data into 2 sets: $\{y>=y_i\}$ and $\{y < y_i\}$ and compute the likelihood of this split using PPO models with a) only an intercept term, and b) both intercept and exposure (i.e. PPOC). The difference in likelihoods of a) and b) gives me a test for this level.  Is this the right approach? Will this test be equivalent to a test for trend? Are there any other approaches that I can adopt to determine if a trend exists between $v$ and $y$ for each of the 3 levels of $y$? PS: I posted this question on the MedStats forum last week. However, there have been no responses thus far. I am hoping someone might have an opinion about this problem here.


title: Feature construction for text mining

text: In the text mining, besides N-gram model, what are the state-of-art models for building feature space while capturing the dependence among the different words, or capturing the semantic meaning in the article?


title: Firth logistic regression in R

text: I'm using the `logistf` package in R to perform Firth logistic regression on an unbalanced dataset. I have a logistf object:               fit = logistf(a~b)      Is there a `predict()` function like on that's used in the `lm` class to predict probabilities for future data points? Or do I have to manually input the estimated parameters from the Firth regression. Thanks!


title: Is there such a thing as a fair die?

text: Is there such a thing as a fair die? On dice where the number is represented by a scooped out dot, surely that makes a difference? Has anyone done any research? In fact thinking about it, why would a coin flip be fair? the physics on each side is completely different.


title: Statistical test for a series of data over time

text: I have a set of data that looks at the number of "hits" a specific program makes over the course of time. The data goes back to September 2010, and includes data up to March 2011, so the data points are monthly. What I want to see if the most recent data (March 2011) shows a statistically significant decrease in the number of "hits" this program makes. I have a feeling there might not be a test that would fit this perfectly, as the data is a bit limited. I can also pull data weekly for the same time frame, which would build 31 points (at which point I would still want to look at the most recent unit for comparison). There hasn't been a population mean built for this data as of yet, as the data can only be pulled as far back as Jan 2010 (but the data from then is not reliable). For reference, just 9 weeks data (as i pulled that first) reveals mean= 1013.67 n=9 st.dev= 53.57 Most recent week= 991 Just eyeballing it does not appear statistically significant as a drop in "hits", however I'll need to perform this analysis every few weeks, and wondering if there's something reliable I can use. Thanks ahead of time for the input!


title: Application of machine learning methods in StackExchange websites

text: I have a **Machine Learning** course this semester and the professor asked us to find a **real-world problem** and solve it by one of machine learning methods introduced in the class, as:   * Decision Trees   * Artificial Neural Networks   * Support Vector Machines   * Instance-based Learning (kNN, LWL)   * Bayesian Networks   * Reinforcement learning I am one of fans of **stackoverflow** and **stackexchange** and know database dumps of these websites are provided to the public because they are awesome! I hope I could find a good machine learning challenge about these databases and solve it. **My idea** One idea came to my mind is predicting tags for questions based on the entered words in question body. I think the Bayesian network is the right tool for learning tags for a question but need more research. Anyway, after learning phase when user finishes entering the question some tags should be suggested to him. **Please tell me** : I want to ask the stats community as experienced people about ML two questions:   1. Do you think tag suggestion is at least a problem which has any chance to solve? Do you have any advice about it? I am a little worried because stackexchange does not implement such feature yet.    2. Do you have any other/better idea for the ML project that is based on stackexchange database? I find it really hard to find something to **learn** from stackexchange databases.  * * * **Consideration about database errors:** I would like to point that although the databases are huge and have many instances, they are not perfect and are prune to error. The obvious one is the age of users that is unreliable. Even selected tags for the question are not 100% correct. Anyway, we should consider the percent of correctness of data in selecting a problem. **Consideration about the problem itself:** My project should not be about `data-mining` or something like this. It just should be an application of ML methods in real-world.


title: Selectively tweaking pairs() axes?

text: ### Question on pairs() I'd like to use `pairs()` to choose a functional form for modeling a set of data. I know that several of my independent and dependent variables are probably lognormally distributed, so I'd like to produce `pairs()` plots where some of my variables are plotted on log axes and some are not. It seems like the par() command might hold the key, but I don't know how to implement it.   * Can I specify log/linear axes behavior when I call pairs()?    * Do I plot `pairs()` first and then modify the resulting graph using par()? Am I on the wrong track entirely? I'm sure I could accomplish the task more easily by transforming some of my variables to log space, but I'm trying to learn something along the way here.


title: What models and software are suited to modelling human decision making?

text: I'd like to model a set of processes. The processes in question are related to human-decision making, so the model will need a measure of input, processing, and then, finally, output. Ideally, the model would be implemented in something like R (which I know quite well) or Python (which I know less well). ### Questions:   * Where is the best place to start with something like this?    * What tools are available?   * What software or language is suited to writing the model in?   * What method would be suited to testing that model against data I've collected from real humans?


title: Random walks and Markov chains

text: Is it incorrect terminology to say that data follows a t-distribution random walk, for example? Does the fact the increments have this underlying distribution mean it should be referred to as a markov chain (and random walk is just a type of mc with a uniform distribution of increments)?


title: Cluster analysis on weighted survey data with continuous and categorical variables

text: I am trying to perform cluster analysis on survey data where each respondent has answered several questions, some of which have categorical answers ("blue" "pink" "green" etc) and some of which have scale answers (rating from 1 to 10 etc). My problem is that certain age groups were over-sampled and I need to weight the data collected in order to accurately reflect the current population. Will it make a difference if I do the cluster analysis on the weighted data, and if so, how do I do cluster analysis on the weighted data? Any advice would be much appreciated! Thanks Emma


title: What methodology should I choose? If hiearchical, what model design?

text: I am working on a problem that has can beyond my level of understanding. I am quite familiar with R, so that would be my preferred choice but I also have access to SAS. **Data** I have created a fake dataset with a couple of rows here My dataset is at the level of insurance "plans". Each plan belongs under a "contract". Each plan has what is known as a "formulary", which details the drugs that the insurance plan covers. Patients belonging to a plan can submit an appeal if they wish for their insurance to cover a drug that is not covered under the formulary. What I am trying to model is the number of appeals requests a plan receives during the year given information about the plan and their formulary. So, for 3097 different plans, I have the following information:   * # of appeals requests received   * contract ID   * formulary ID   * enrollment (~ number of patients enrolled in the plan for the year) There are 584 different contracts. Some contracts have multiple plans under them, some only have one plan. There are 297 unique formularies. Again, some formularies belong to only one plan, some formularies are shared by multiple plans. Each unique formulary has two main metrics of interest:   * # of drugs covered   * % of drugs covered that require prior authorization **Model** I immediately thought of using Poisson regression, since the response variable is the number of events occurring during a fixed time interval. There also seems to be over-dispersion so I have looked into negative binomial regression instead. Also, since the units are of differing size (enrollment) I am using the offset option. The problem arises with the fact that the request rates of plans that share formularies or contracts are probably related to each other and I don't exactly know how to capture this relationship. There also may be issues with grouping using contract and/or formulary due to the fact that quite a few contracts and formularies belong to only a single plan. I initially tried designing a hierarchical model (I have limited experience with Gibbs sampling) but I can't seem to design a model that works. I have almost no experience with mixed-effects models and no matter how many times I read about them I can't see how exactly I would set it up for my problem. I am mostly interested in the effect that **# of drugs covered** and **% of drugs covered that require prior authorization** have on rate of requests for a plan, but I understand that many plans will share covariates and the errors of plans within the same contract may be correlated.


title: Evaluating mean differences between multiple groups based on categorical factors

text: I am working on a very large dataset of patients who belong to one of five groups. I want to test if groups are not so much different than each other in terms of age, race, sex, smoking level, weight, and some other factors. Some factors are numerical and some are categorical. I know I can use ANOVA to compare mean differences between groups for numerical factors. However, I do not know how I can test mean differences for categorical factors. Some references say that Chi-Squared test can be used for this purpose, but it seems that `chisq.test` in `R` is only for testing either independency or goodness of fit. Could someone write a sample code in `R` showing how the mean differences between more than two groups can be evaluated for a categorical factor? Is any pre-evaluation required to check if a specific test is valid to be used? Thanks in advance.


title: understanding of libsvm output

text: I applied libsvm to build a text classifier. The output looks like as follows:               w[a_b] = 0.006     w[a_c] = 0.032     w[a_ctif] = 0.000     w[a_cs] = 1.009     w[aa_d] = 0.000     w[a_e] = 0.001      The terms, e.g., a_b, and a_c, etc are the features generated for characterizing these documents. My question is how to understand these values corresponding to each term? The model output shows that all of these values are bigger than zero. Can I say that a positive file tends to have a term with higher weight?


title: Stata ANOVA Interaction term differences

text: The following two statements produce different results.               /* DV: "price" e.g., from 0 to 100 */     /* IV: "treatmentA", "treatmentB" both binary variables */          /* statement 1 */     anova price treatmentA##treatmentB          /* statement 2 */         generate A_x_B = treatmentA * treatmentB     anova price treatmentA treatmentB A_x_B      Do the two statements not represent the same model? Any suggestions, why generating an interaction term by myself and including it into an ANOVA produces different results than using the built-in option to include interaction terms? edit: some additional output:               . anova price treatmentA##treatmentB                                     Number of obs =    2160     R-squared     =  0.0481                                Root MSE      = 3.62039     Adj R-squared =  0.0468                            Source |  Partial SS    df       MS           F     Prob > F         ---------------------+----------------------------------------------------                        Model |  1427.66343     3  475.887809      36.31     0.0000                              |                   treatmentA |  590.372454     1  590.372454      45.04     0.0000                   treatmentB |  786.261343     1  786.261343      59.99     0.0000         treatmentA#treatme~B |  258.622454     1  258.622454      19.73     0.0000                              |                     Residual |  28259.1958  2156  13.1072337            ---------------------+----------------------------------------------------                        Total |  29686.8593  2159  13.7502822              . generate A_x_B = treatmentA * treatmentB          . anova price treatmentA treatmentB A_x_B                                     Number of obs =    2160     R-squared     =  0.0481                                Root MSE      = 3.62039     Adj R-squared =  0.0468                            Source |  Partial SS    df       MS           F     Prob > F                  ------------+----------------------------------------------------                        Model |  1427.66343     3  475.887809      36.31     0.0000                              |                   treatmentA |       33.75     1       33.75       2.57     0.1087                   treatmentB |  730.034722     1  730.034722      55.70     0.0000                        A_x_B |  258.622454     1  258.622454      19.73     0.0000                              |                     Residual |  28259.1958  2156  13.1072337                     ------------+----------------------------------------------------                        Total |  29686.8593  2159  13.7502822


title: R Package for Treatment Effect Analysis?

text: I could not find a package in R about Treatment Effect Analysis. So is there a R package about Treatment Effect Analysis? This means, estimating the average treatment effect, average treatment effect on the treated using different methods, including algorithms such as k-nearest-neighbour matching, caliper-matching and things like that? I found single packages about k-nearest neighbour matching and matching itself (e.g. matching of FNN), but I could not find a comprehensive package about treatment effect analysis.


title: Regression interaction w/ dichotomous predictors: Two levels, two coefficients?

text: this question is about a way of presenting the results of an interaction I thought of, and I'm looking to see if it makes sense. If there is a better way to do what I'm trying to do (or if what I'm trying to do is not advisable), please let me know. To simplify, let's say I've got a multiple linear regression equation with two dichotomous predictors (dummies) and an interaction between the two--let's say the DV is test score, predictor 1 is gender (M/F) and predictor 2 is course type (experimental vs. traditional). I understand that in such equations, the coefficient for each non-product term assumes the value of the other non- product term is zero. For the special case of dichotomous predictors only, this means that each coefficient assumes one of the two values of the other: so for example, the coefficient for gender applies to _either_ the experimental or the traditional course, but not both as it would in an equation without an interaction. Which value of the other term is assumed depends on how the data are coded. So, my question is: if gender is the focal predictor and course type is the moderator, and assuming the interaction is significant, would it be valid to run the model once with course type coded as (0=experimental, 1=traditional) and again with these values switched to (0=traditional, 1=experimental)? Then I could report two coefficients for gender, each of which corresponds to one value of course type. This would be helpful because it would demonstrate the distinct relationships of gender to test score within each course, whereas a single coefficient would describe this relationship in only one course. From what I understand about interactions, this seems to make sense. But I have no idea whether this is something that's ever done, or if there's a better way to show the coefficients of a single dummy variable at the two levels of another. Thanks in advance!


title: Interaction Test with Odds Ratio

text: Reading Pocock et. al on subgroup analysis, the authors describe calculating an interaction test given 2 odds ratios (and the count data used to get them). Specifically,               mat <- matrix(c(22,12, 11, 11), nrow=2, ncol=2, byrow=T)      Where `22, 12` correspond to deaths/not deaths in one group and `11, 11` are same for another group. They report ORs of 2 and 1 respectively, and then report that the interaction test comparing odds ratios has _p_ = 0.21. What is the interaction test they performed?


title: Linear regression coefficient limitation

text: I have built a linear model using `lm` in R               lm(age ~ sibsp + parch + pclass, trainset)      Here `sibsp` is the number of siblings/spouses aboard the Titanic, `parch` is the number of parents/children aboard and `pclass` is the passenger class (1st, 2nd, 3rd). My problem is that when I run a prediction on a x (from a testset) where variable `sibsp > max(sibsb in trainset)` I get a negative prediction for the `age`. Any idea how I could work around that problem? I thought of limiting any `x$sibsb` to the `max(sibsb in trainset)` in the model but I'm not sure how to implement that in R. Any comments or other solutions welcome.


title: Bioinformatical problem - specific word enrichment in a given sequence

text: My bioinformatical problem looks like this:   two sets of gene sequences (size of set A > 1000, size of set B > 1000; sequence length varies from 1000 to 100000). > Set_A_Sequence_1: ACGTACGTACGT...   >  Set_A_Sequence_2: ACGGAAGT **AAA** T...   >  ....   >  Set_B_Sequence_1: **AAA** G **AAA** TG **AAA** ...   >  Set_B_Sequence_2: **AAA** TC **AAA** C **AAA** ...   >  ... I want to see if Set_B sequences are enriched for specific word (for example: `AAA`).   How can I do this?   I came up with three solutions:   1. Count how many Set_A sequences have/don't have word;   Count how many Set_B sequences have/don't have word.   Apply Fisher's test.   2. Count how many word occurrences there are per sequence in Set_A;   Count how many word occurrences there are per sequence in Set_B.   (For a given set example that would be: Set_A:0,1; Set_B:3,3).   What statistical test I can use for such enrichment analysis?   3. Calculate percentage of sequence in Set_A covered with word;   Calculate percentage of sequence in Set_B covered with word.   (For example, data would look like this: Set_A:0%,25%; Set_B:75%,75%).   What statistical test I can use for such enrichment analysis? **Questions:**   Is it right to use Fisher test in solution 1 (Contain/Don't contain word)?   What statistical tests I could use for solution 2 (Number of words)?   What statistical tests I could use for solution 3 (Coverage with word)? **Edit**   Simplified data looks like this:               Sequence name   Length   Contain word(0/1)   Number of words   Coverage with word(%)            Set_A_seq_1     1000               0                0                 0     Set_A_seq_2     2000               1                1                 15     Set_A_seq_3     3450               0                0                 0     Set_A_seq_4     10000              0                0                 0     Set_A_seq_5     25000              1                2                 5     ...          Set_B_seq_1     20000              1                3                 25       Set_B_seq_2     100000             1                3                 30     Set_B_seq_2     9000               1                5                 70     Set_B_seq_2     10000              1                10                85     Set_B_seq_2     12000              1                7                 60     ...                **EDIT** I wasn't able to find a lot of published methodology for genomic site enrichment, but this figure suggests perfect way of solving problem that I have. Figure A. - Enrichment in three different genomic sites compared using permutation and odd enrichment to permutated data. ![Relationship of repetitive elements to specific genomic sites \(22948768\)](http://i.stack.imgur.com/KfDJh.jpg)


title: One vs All and One vs one in svm?

text: What is the different between onee vs all and one vs one SVM classifier?? Is One vs All mean = 1 classifier to classify all types /categories of the new image and one vs one mean= each type /category of new image classify with different classifier (each category is handled by special classifier)?? For example if the new immage to be classified into circle, rectanglar, triangular..


title: The organization as a unit of analysis and sampling issue

text: I'm doing a quantitative study entitled "Technology Transfer and Competitive Advantages in Oil and Gas Companies". In my study, the unit of analysis is the organization and there are only 9 oil and gas companies in my country. The respondents will be randomly selected from the managerial level (stratified) based on the size of each company. The total number of employees in managerial level is 3,336 and the sample size is 959. I'm going to use Structure Equation Modeling for the analysis. My question is, taking into consideration the unit of analysis is the "organization", does the sample size "959" represent the population (9 companies), and is it enough? Is this the right methodology? My friends told me that if the unit of analysis is organization and I have only 9 companies, so one person from each organization will represent the organization. That means I will have only 9 respondents, which is not enough for the study.


title: Hierarchical regression using residuals

text: Edit: Perhaps my first asking (see below) was a little too ambitious. The crux of the question is this: in a linear regression setting, is there anything statistically invalid about regressing the residuals of a previous model on new predictors (assuming the observations are IID, plenty of degrees of freedom, etc)? There's a bit more about why one might want to do this below, along with some other possible alternatives to achieving the same aim; but since the above is something I actually did, I'd be happy to start with knowing how wrong it was. * * * Suppose a standard regression situation, if there is such a thing: IID data, a set of (not necessarily orthogonal) predictors, and the desire to figure out how some of the latter relate to the former. For sake of argument, let's say there are 9 predictors and 90+ independent observations, with the 9 predictors falling into 3 separate groups, based on some theoretical ordering of their relative importance, $y = a_1 + a_2 + a_3 + b_1 + b_2 + b_3 + c_1 + c_2 + c_3$. A standard approach to assessing the impact of the $b$ group after controlling for the $a$ group would be to use hierarchical regression, and compare the fits with an ANOVA; and likewise for the $c$ group. However, when you do this, the actual parameter estimates you end up looking at are those that come from whichever stage makes the most sense--so, if none of the $c$ group variables were significant (or more precisely, if there was a nonsignificant change in model fit when adding those variables), but some of the $b$ group were, you might consider your "official" estimates for the $b$ variables as those coming from the model $y = a_{1:3} + b_{1:3}$ (adopting a bit of shorthand). The issue is that those $b$ estimates were estimated simultaneously with the $a$ variables, so the $b$ estimates themselves and their individual significances aren't being strictly controlled for $a$, as far as I can tell. In any situation where you want to be fairly conservative in ascribing explanatory power to those $b$ variables, it seems like this is only giving you partial information--that is, you know the $b$ group as a whole is improving the fit, but not much beyond that. Off the top of my head, I can think of three possible solutions (there are surely more, and some or all of these might be completely idiotic). The first is to run multiple hierarchical models: for instance, maybe $y=a_{1:3}$ followed by each of $y=a_{1:3}+b_i, i\in \{1:3\}$ (or $y=a_{1:3}+b_{i}+b_{j}, i\ne{}j \in \{1:3\}$ followed by $y=a_{1:3}+b_{1:3}$) and drawing conclusions about $b_i$ from the model in which $b_i$ was first added. However, that doesn't solve the fact that the coefficients are all estimated simultaneously, so the $b$ estimates aren't "pure". Second is to orthogonalize groups of regressors w.r.t. one another: for instance, $y=a_{1:3}+b_{1\perp{}a_{1:3}}+b_{2\perp{}a_{1:3}}+b_{3\perp{}a_{1:3}}$. Now the estimates for the $b_i$s must reflect variance distinct from the $a_i$s, and presumably their significances are correct. For some reason this approach worried me, just because it feels a little voodoo to see your predictors changing around (especially when it comes to interpreting results), but the more I've thought about it, the more it's seemed possibly okay. Third, and the approach I settled on because of those concerns, is to fit $y=a_{1:3}$, then fit $y-\hat{y}_{a}=b_{1:3}$, and so forth--that is, to fit each subsequent model to the residuals of the previous model, of course always including only the next group of predictors. I can imagine there might be a problem with the degrees of freedom, insofar as each model after the first presumably has less and less residual variance (assuming earlier stages explain at least something), but the degrees of freedom only reflect the current batch of variables. That issue notwithstanding, I have three basic questions: 1) is there anything fundamentally wrong with this procedure? 2) does it address the question we're trying to address? 3) how does it compare to the previous two suggested approaches? Thanks!


title: Difference between a SVM and a perceptron

text: I am a bit confused with the difference between an SVM and a perceptron. Let me try to summarize my understanding here, and please feel free to correct where I am wrong and fill in what I have missed.   1. The Perceptron does not try to optimize the separation "distance". As long as it finds a hyperplane that separates the two sets, it is good. SVM on the other hand tries to maximize the "support vector", i.e., the distance between two closest opposite sample points.    2. The SVM typically tries to use a "kernel function" to project the sample points to high dimension space to make them linearly separable, while the perceptron assumes the sample points are linearly separable.


title: How do I sample multiple normally distributed items and maintain a fixed sum?

text: Suppose jars are being filled with Red, Green and Blue balls in an approximately fixed ratio. The total number of balls that fits in a jar is fixed (say N=1000). I know the mean and standard deviation of the ratios; let's say means for (R,G,B) are 0.5,0.4,0.1 and stdevs are 0.05,0.05,0.02. They appear to be approximately normally distributed. How can I construct samples from the distributions of R+G+B balls such that they sum to N? If I sample R and G independently and calculate B = N-R-G, then G may not follow the original distribution and may even be negative. If I sample independently from R, G and B, the total will not be N; if I renormalise the sample after, they will no longer follow the original distribution. Any thoughts welcome.


title: Generalized LM or LM in ecological dataset

text: I have previously posted a question for the same dataset but now I had somme issues with the models and I wanted to re-phrase my question. My dataset contains 50 morphometric characters (which we reduced with factor analysis or pca to few common components) measured on the roe deer skull. I want a model for predicting skull dimensions from the absolute areas of the forest and plowland in habitats they live in. My first idea was to enter absolute areas as predictors but I was advised to use them as proportions and to add population abundance as a predictor in the model. I know that GLM (with binomial error structure) is great for proportion data but I am unsure how to specify such a model in R. In fact I am also unsure how to enter proportions (I can calculate percentages, but am wondering if they must total to 1, 100%). Any ideas? **SAMPLE DATASET**                  Factor1 population manage foraging height biome abundance area  forest plough       -0.6033788 ADA_BEC   best    fields   low    agS    1500    73154  61154  12000       0.3250981 ADA_BEC   best    fields   low    agS    1500    73154  61154  12000       0.5577059 ADA_BEC   best    fields   low    agS    1500    73154  61154  12000      -0.1596194 PM        good    plains   med    kgS    23980   856251 89499  579870      -1.3089952 PM        good    plains   med    kgS    23980   856251 89499  579870      -2.1693392 SP        poor    mount    high   hgS    2500    65872  38000  47098      -0.9669080 SP        poor    mount    high   hgS    2500    65872  38000  47098      -1.8857842 SP        poor    mount    high   hgS    2500    65872  38000  47098       0.7242678 DKN       best    fields   plain  agS    65908   989981 181133 12400       1.6815373 DKN       best    fields   plain  agS    65908   989981 181133 12400      Area is the total area and forest is the forest and plough area (which don`t add up to the total area as there are meadowlands and urban areas, we have this data too). These areas are in expressed in Ha. Factor1 are factor axis scores and abundance is the total number although we have density (individual per area) also.


title: mob model tree algorithm

text: I am trying to figure out the inner workings of the mob function in the party package. I can't figure out how the splitting variable is selected when it is a categorical variable. In the publications by the authors, they say that it is done by a chisquare test of association between the residual deviances and the categorical variable. I can't see how this works. In the toy example below, the mob function works fine. The variable X is the splitting variable.               b <- (1:1000) / 100 - 5     a <- c(b,-b)     b <- c(b, b)     X <- as.factor(c(rep(FALSE, 1000), rep(TRUE, 1000)))     mob(a ~ b |  X)      However, if I fit a linear model without splitting, there is (by construction) no association between the residuals and the values of the splitting variable:               res <- sign(residuals(lm(a ~ b)))     table(res, X)      The resulting table is:               X     res  FALSE TRUE      -1   499  500     1    501  500      This table does not show that the X variable is important, still the algorithm can figure it out. How is that done?


title: Repeatedly measured, longitudinal and ordinal data

text: Description of the problem: The effect of a certain treatment is tested as follows: For n subjects, the endpoint is measured three times before and after the treatment is administered, respectively. The endpoint is modeled as an ordinal variable with outcomes A,B,C and order $A<B<C$. There are two questions to be answered: 1) How precise is the method measuring the endpoint? 2) Is the treatment effective? Statistical description of the problem: We consider a random variable $Y_{itk}$ with $i=1,\ldots,n$ the subject index, $t=1,2$ the time point index, and $k=1,2,3$ the index of the repeated measurements. The random variable $Y_{itk}$ has the three outcome possibilities A,B,C which can be ordered as mentioned above. Also, the random variables are independent in $i$ but, of course, neither in $t$ nor in $k$. Regarding the two questions: 1) What would be a good precision measure? I thought about the probability that a) all three results are the same, b) two are the same, c) all three measurements are different. For example the first probability could be estimated as follows: $\frac{1}{2n}\sum_{i,t}\mathbf{1}_{\{Y_{it1}=Y_{it2}=Y_{it3}\}}$. However, do I need to adjust of the dependence of the measurements at times $t=1$ and $t=2$? What other precision measure could be considered? 2) Concerning the comparison of the endpoint between the two time points, I have several ideas. 2.1) Defining a ranking of the outcomes of $(Y_{it1},Y_{it2},Y_{it3})$ and comparing the the distribution of the ranks for both time points. The things that bothers me is that the ranking system is up to a certain degree arbitrary and the results might depend on the ranking system. 2.2) Taking the median of the outcomes $\{Y_{it1},Y_{it2},Y_{it3}\}$. In the end, one could compare the distributions of the medians at times $t=1$ and $t=2$. 2.3) Multinomial regression with a random predictor. Eventually, the sample size $(n=10)$ might be to small for a fancy model. Any advice is highly appreciated. Thanks! By the way, a similar question has already been asked here. An answer hasn't been given though.


title: P values seem to be the wrong way around

text: Say the null hypothesis is that two samples emanate from two distributions with either equal medians or equal variances. And p values are the probability that a value or more extreme comes from the original distribution given the null hypothesis holds. Therefore, why do correlations with low p-values mean that they are strong? It seems that holds if the null hypothesis is that they do NOT come from similar distributions. Can you please put me straight? Regards,


title: Repeated measures ANOVA with replicated measurements.

text: I have performed an experiment testing 3 different modes of operation on 16 participants. Each mode of operation was tested twice per user.               User    Mode    Trial1Time  Trial2Time     1       1       20          30     1       2       5           7            1       3       40          25           2       1       10          20     2       2       15          17           2       3       30          35     3       1       13          26     3       2       11          11           3       3       35          38      I have run repeated measures ANOVA on my data using only the Trial1Times. I would like to use both Trial1 and Trial2 in the analysis however as I have such a small number of participants.   1. Is this possible with repeated measures?   2. If so, does this add statistical power to the results?   3. (Bonus) How would I go about doing this in SPSS? Appendix: The data formatted in a different way (The actual numbers are different too, but I just made it up as an example)               User    Mode1   Mode1   Mode2   Mode2   Mode3   Mode3             Test1   Test2   Test1   Test2   Test1   Test2     1       10      15      20      25      30      35     2       12      12      23      21      31      33     3       14      16      27      26      35      39


title: Point predictions and prediction intervals

text: ![enter image description here](http://i.stack.imgur.com/Sgd9z.png) I am having trouble getting started with this question. Is there a formula for finding point predictions and prediction intervals? Any advice or hints would be great! Let $\varepsilon$ be defined as white noise. Does $X_{81}$= 0.8*$X_{80}$ -0.2* $\varepsilon_{80}$?   Does $\varepsilon_{81}$ =$X_{81}$-0.8*$X_{80}$- 0.2*$\varepsilon_{80}$? Does $X_{82}$ = 0.8*$X_{81}$ -0.2*$\varepsilon_{81}$?


title: Interpretation with training and test set with standardized variables

text: I've standardized all the variables (even the response variable) and then I've split my data into a training and test part. And for example, I've got the following model based on my TRAINING set: y = 0.5x_1 - 0.2x_2 Now I shall get values for y like -0.4,0.7,... but I want to say something about the original response variable. What can I do to get this?


title: How to represent this problem in a useful way?

text: How can I represent this problem as a function. Bear with me because it may seem odd. Say, for example I make cakes. There are 30 variations of cake in one batch. For each batch there are 19 luxury set ingredients and a lot of secondary ones, of which there must be 100% of the luxury ingredients in at least 70% of the cakes to pass quality control. Who carry out random checks. But procurement tell me that the ratio of luxury ingredients to cake is too high. There needs to be less luxury ingredients in each cake, but still maintain evidence that 100% of the luxury ingredients go into 70% of all the cakes. So in order to solve the problem, I propose to make more variations of cakes in each batch, incorporating less luxury ingredients into each one. But in a random test of my batch, there must be 100% of the luxury ingredients found during the test, regardless of the cakes selected. SO could anyone help me with finding: The ratio of ingredients to cake when there are 30 variations and 100% of the luxury ingredients need to be evidenced in 70% of the cakes. In order to understand the nature of the complaint. If I increase the variety of cakes in a batch by some variable x, how many variations do i have to invent to make the ratio come down by 1 luxury ingredient, and express this as a function so I can experiment with bringing down the number of luxury ingredients in each variety to something more manageable. Now if I've left any glaring information out there, let me know.


title: Question on Machine Learning Overview?

text: I am a researcher from India and working in the field of Computational Linguistics for quite some time. I have lately started working in the field of Machine Learning based algorithms. To do this I tried to master over statistics, and tried to read the book by Ethern Alpaydin, attended lectures in Coursera by Professor Andrew Ng Micael Collins, discussed the problems with my teachers and colleagues. As they kindly suggested me I tried to do that, get a good overview on four or five algorithms and work out one or two in detail. So I tried to know Regression-Univariate/Multivariate, Decision Tree, SVM, CRF and worked out Naive Bayes and HMM with practical examples and trying to juggle around by changing variables and parameters. Of late I was thinking what may be said as the command over these two models and bit of ML? What questions should I address and how they may be tackled? I know NLTK, Scikit learn and more or less fluent in Python. As this room full with ML experts if any body may kindly guide me, an online resource may be of great help. Thanking in Advance, Subhabrata Banerjee.


title: How many different approaches (or schools, or methods) is in statistics?

text: I am beginning to learn statistics and made some searches in Internet about statistics. Most of the articles I found says that there are two different approaches (or schools, or methods) in statistics Frequentist and Bayesian. But one article says that there are also Information-Theoretic and Likelihood. So how many different approaches (or schools, or methods) are in statistics, 2, 3, 4 or may be more?


title: Include interaction in multiple imputation - r

text: I'm doing some imputation models of time until recurrence of tuberculosis (Cox model). This model should include an interaction between the time and the outcome of the previous episode of disease (0- cured, 1-default). I realized that one should include all information before imputation, in order to built a more robust imputed model. However, the best way to deal with it seems quite challenging. One of the option I'm trying is to split the dataset (in this case, by the categorical variable), impute each subset alone and append both subsets to further analyze. My question is **after imputing each subset separately, I found several differences between the model without interaction and the model imputed separately. Am I doing everything right?** I was not expecting a variable really significant to not be significant anymore (just an example). Can anyone tell me if this is the correct way to include an interaction? **And why are the results so different?** I've compared the imputed values for each variable in each of these two methods and they are not very different, this is one of the reason that puzzles me for having such a different model.. Here is my code to impute each subset separately:               dataCured <- data1[data1$Sit == 0,]         dataDefault <- data1[data1$Sit ==1,]          impCured <- mice(dataCured, m=70, maxit=10, pred=quickpred(dataCured, method="spearman"), method="pmm", seed=71152, printFlag=TRUE)          impDefault <- mice(dataDefault, m=70, maxit=10, pred=quickpred(dataDefault, method="spearman"), method="pmm", seed=71152, printFlag=TRUE)          imp <- rbind.mids(impCured, impDefault)


title: Standard Error of the Correlation Coefficient

text: As defined here, the estimation of the coefficient of correlation is $$r = \frac{\Sigma (X_i-E[X])(Y_i-E[Y])}{\sqrt { \Sigma (X_i-E[X])^2 \Sigma (Y_i - E[Y])^2}}$$ and the standard error of $r$ is $$SE_r = \sqrt{\frac{1-r^2}{n-2}}$$ **How was $SE_r$ calculated?**


title: Equivalence tests for non-normal data?

text: I have some data that I can't necessarily assume to be drawn from normal distributions, and I would like to conduct tests of equivalence between groups. For normal data, there are techniques like TOST (two one-sided t-tests). Is there anything analogous to TOST for non-normal data?


title: Correlation coefficient: If $\rho = 0$, then $r$ is normally distributed with mean 0. Why?

text: From this source, the estimation of the coefficient of correlation is $$r = \frac{\Sigma (X_i-E[X])(Y_i-E[Y])}{\sqrt { \Sigma (X_i-E[X])^2 \Sigma (Y_i - E[Y])^2}}$$ If the coefficient correlation is null ($\space\rho = 0\space$), then its estimator $r$ is normally distributed with mean equal 0. **Is this valid whatever is the distribution of $X$ and $Y$? Why?**


title: Deterministic components in covariates/exogenous variables in time series models

text: Actually, I have read a pair of books about time series analysis, but I am still not sure about how to treat deterministic components, like trend and seasonality, in the exogenous variables in a time series model. Do I have to detrend and deseasonalize the covariates before I use them as axplanatory variables in a time series model? I would also be thankful for a reference. Thank you in advance!


title: Correlation between matrices in R

text: I have problems in using the `cor()` and `cor.test()` functions. I just have two matrices (only numerical values, and the same number of row and columns) and I want to have the correlation number and the corresponding p-value. When I use `cor(matrix1, matrix2)` I get the correlation coefficients for all the cells. I just want a single number as result of cor. In additon when I do `cor.test(matrix1, matrix2)` I get the following error               Error in cor.test.default(matrix1, matrix2) : 'x' must be a numeric vector      How can I get p-values for matrices? You find the simple tables I want to correlate here: http://dl.dropbox.com/u/3288659/table_exp1_offline_MEANS.csv http://dl.dropbox.com/u/3288659/table_exp2_offline_MEANS.csv


title: Autocorrelated predictors in linear models

text: I need to predict the outcomes of a time-series variable $Y$ based on two time-series predictors $X1$ and $X2$. For simplicity I will only illustrate $X1$ in the rest of this question. The correlation between the predictors and $Y$ at lag 0 is relatively high as illustrated by the cross-correlation function. However the lag 1 correlation is a lot lower and does not seem to have much prediction power:               > ccf(as.numeric(X1), as.numeric(Y), lag.max = 3, plot = FALSE, na.action = na.pass)          Autocorrelations of series 'X', by lag              -3     -2     -1      0      1      2      3       0.015  0.008 -0.001 -0.415 -0.018 -0.002 -0.002      Indeed when I regress $Y_{t+1}$ against $X1_{t}$, $X2_{t}$, $R^2$ is very low. Now, looking at the $X1$ partial auto-correlation indicates that the predictor itself could be estimated using an AR model. The significance level for this series is around 0.01:               > pacf(X1, plot = FALSE, na.action = na.pass, lag.max = 15)          Partial autocorrelations of series 'X1', by lag             60   120   180   240   300   360   420   480   540   600   660   720   780   840   900      0.171 0.103 0.067 0.055 0.046 0.049 0.027 0.032 0.029 0.012 0.020 0.037 0.020 0.026 0.020      To resume, the correlation between $X1$ and $Y$ at lag 0 is high (but not for $X1_{t}, Y_{t+1}$) and $X1$ is significantly auto-correlated up to 10/15 periods. What would be the best approach to predict $Y$ and why?   1. Estimate $X1_{t+1}$ using an AR(15) model, then use the estimated value as a predictor for $Y$.   2. Fit all $X1_{t-15, ..., t}$ against $Y$.   3. Use an exponential moving average of $X1$ ($\alpha = ~0.125$) as the predictor.


title: Is there a rule of thumb on the relationship between the number of instances and the number of features?

text: If we build a classifier based on a very small number of instances (say, fewer than 300) and the number of features we are using is very large (say, larger than 100k features). If we decide to introduce a feature selection step before building the classifier, is there a rule of thumb on how many features we should choose?


title: Test of multicollinearity among independent variables in logistic regression

text: I am using 10 independent variables in building logistic regression model. I am sure that some of these variable are correlated. Can anybody tell me how to check for multicollinearity among independent variables in this case. Thanks!


title: How to implement a two-stage hierarchical model of time series data in R?

text: I'm currently working with a data set that consists of a monthly case count for several sites, along with a number of site-specific covariates. We're trying to estimate the effect of one of them on the case-load (either in terms of numbers or a percent). At the moment, I've been using Peng and Dominici's _Statistical Methods for Environmental Epidemiology with R_ as a reference for the approach I'm considering using. For a single site, the idea is to fit a Poisson model with roughly the following form: $$\log ({\mu _t}) = \alpha + \beta {x_t} + \eta {z_t} + f(t) + {\varepsilon _t}$$ Where `mu` is the monthly case count, `beta` is the log-relative increase in cases due to a one unit increase in the variable of interest, `eta` is the vector of other measured covariates and f(t) is some sort of smoothed spline over time. So far, I'm comfortable enough with this. My problem is implementing a hierarchical model that gives an estimate of the effect pulling from information from each site. Is there a particularly good way to approach this in R? A well-documented package for example? Is it sufficient to essentially treat this as a meta-analysis with each site as a "study" and use something like `metafor`? And is there a way to do it while modeling a site-specific covariate? The authors seem to be pushing for the use of a Bayesian hierarchical model, but the implementation of something like that is admittedly a bit of unexplored territory for me.


title: BIC and AIC(c) and group data

text: I would like to compare two models using the BIC and AICc. Doing so seems fairly straightforward if both models are fit to only one dataset. However, I have data from 10 participants, and there is no reason to assume that they are all well-described by the same set of parameters. Is it appropriate to sum the log-likelihoods of a model across all participants and multiply the number of parameters by the number of participants to obtain a "group BIC/AICc" for that model, or is there another way of obtaining the BIC/AICc in such a situation?


title: ARIMA, adjustments and intervention analysis

text: I have very little knowledge of time-series analysis (despite my stat master - didn't do anything else than an introductory course) but now I'm facing a statistical problem whose answer is this very kind of analysis - so would really need a helpful hand. In a nutshell, I have account monthly sales. In only _some_ of them (1 to 10 or so), a marketing project has been conducted (I call them the Affected accounts); the other accounts remain unaffected (called the Control accounts). The question that I have to answer is basically "did the project have an impact in the Affected sales (but here's the part that annoys me) compared to the Control sales?" So, as I have a time-series analysis (and a seasonality) I guessed that I had to perform a proc ARIMA (I'm using SAS - and its fantastic feature that I discover and displays an automatically fitted model, Time Series Forecasting System). Here are my questioning:   * How is an adjustment variable accounted for, in an ARIMA model ? If I make an educated guess from the linear regression model, I assume I just have to enter the variable in the INPUT, as independent / explanatory variables ? Also, if I am to enter my variable "Control_Sales" like that, is my adjustment correct ? Shouldn't I do a preliminary transformation of my dependent variable, for example redefining it as (Affected_Sales/Control_Sales) or (Affected_Sales - Control_Sales) ?   * Besides, I have another series of questions around the Intervention analysis - because the conduct of the projects is never a straightforward, pulse intervention and I cannot really pin them down to a precise type. Here's the best explanation I can provide: The projects, as I already said, are not carried out in all accounts, but only in a very small portion of them. Moreover, the projects can be regarded as marketing campaigns that last several months. They have a start date and an implementation(=end) date. But even from my supervisors, there is disagreement on when we should be starting to assess an effect on the sales - sometimes, the sales start to increase _during_ the campaign, most of the time, the effect is seen after the end of the project, but not straight after, perhaps from 2 to 9 months after... So I'm getting very confused when it comes to modelling the type of intervention (a ramp ? several pulses ?) and I am wondering if it wouldn't be better if I performed an Intervention Detection, to see where the impact actually took place, and in the discussion part of my analysis compare it to the reported, "theoretical" implementation date. I deeply believe this idea would be the most appealing, due to the fact that nobody is able to be clear about the Intervention start point (statistically speaking). What do you think ? How would you do ? But there comes the most challenging part: I have no clues whatsoever how to perform an Intervention Detection... I could really do with some piece of advice ! PS. I also meant to say that I am happy to provide the data if it can help with the discussion. * * * **Update** ![So here is the data](http://i.stack.imgur.com/avcFk.png) Here is my data, if that can help. Ideally I would stick to the ARIMA (or other, best-fitting time-series model) because I not comfortable at all with this field and this is something I am currently capable of doing (don't know anything about VAR or Granger causality). Moreover, I have only collected the sales from Jun-13 onwards, as I had fixed myself a 6-month limit for the "sale sbefore", but I can have figures back to 2011... I don't know however how relevant it is to include them in the analysis - and if that couldn't "dilute" the impact of the intervention ? * * * **Update 2** Data from 2011 onwards                   Date    Sales_nat   Sales_affected         Jan-11  13535.04614 10564.2         Feb-11  12255.18701 6338.52         Mar-11  15504.88513 16902.72         Apr-11  13259.76914 14085.6         May-11  15967.85091 13381.32         Jun-11  15351.15898 9859.92         Jul-11  16001.81365 16902.72         Aug-11  20151.51071 23692.09         Sep-11  21533.29437 30507.47         Oct-11  21122.32893 19242.99         Nov-11  22350.66487 25579.51         Dec-11  21707.95193 15019.31         Jan-12  23225.30391 28394.63         Feb-12  22782.53005 23466.67         Mar-12  24346.6397  30030.61         Apr-12  23093.62005 21361.83         May-12  26336.53924 22530.96         Jun-12  22695.90797 18770.13         Jul-12  26825.00843 21824.68         Aug-12  26202.68137 23225.24         Sep-12  24917.01741 23929.52         Oct-12  30170.13777 20649.55         Nov-12  28223.8397  27215.49         Dec-12  28165.34954 19713.84         Jan-13  30716.72604 20182.69         Feb-13  26684.0364  28867.48         Mar-13  28277.98102 15019.31         Apr-13  30858.22655 28159.2         May-13  31845.99066 24871.22         Jun-13  29444.00066 31425.17         Jul-13  34896.64914 41989.37         Aug-13  30925.31929 23945.52         Sep-13  32861.02081 25577.51         Oct-13  34976.22452 10795.63         Nov-13  32547.14685 31668.6         Dec-13  38381.84523 35435.43         Jan-14  40211.59741 4221.68         Feb-14  31925.53772 29569.76         Mar-14  37865.2967  49251.59         Apr-14  39391.11061 28865.48         May-14  35614.36797 7042.8         Jun-14  41398.94482 37316.84


title: Regression with interaction term that is a function of other regressors

text: Consider the simple interaction equation where moderator M determines the effect of X on Y. $$ Y=\beta_1X+\beta_2M+\beta_3XM \quad (1) $$ then the slope on X is $$ \frac{dY}{dX}=\beta_1+\beta_3M \quad (2) $$ Fine. Now suppose that $M=X/Z$. You might say that M is just a linear transformation of X, but if X and Z are different for every observation in a data set, M provides some additional information beyond X (I think). But now our interaction equation is $$ Y=\beta_1X+\beta_2\frac{X}{Z}+\beta_3\frac{(X^2)}{Z} \quad (3) $$ and the slope on X is $$ \frac{dY}{dX}=\beta_1+\beta_2\frac{1}{Z}+\frac{(2\beta_3X)}{Z} \quad (4) $$ which is very different than in eq (1). My (maybe ill formed) question is whether it is even acceptable to formulate a regression in this way, and what the implications on interpretation will be. I have a fairly well-formed theory that goes along with the equation, but I just can't seem to make sense of the math. Any insights are appreciated.


title: Struggling to make headway with Introduction to Probability by Bertsekas & Tsitsiklis

text: I'm an academic in biology and I'm trying to improve my working knowledge of probability and probabilistic methods. I was recommended the Introduction to Probability textbook as an excellent introduction to the subject. I am, however, finding it extremely hard to make progress with the exercises. I've been working through the exercises at the end of first chapter and it has taken me several weeks to get to the end. I usually spend a couple of hours on each question and invariably still end up with the incorrect answer. I studied A-level mathematics and statistics so I feel confident with the algebra. I also understand the examples given in each chapter. It is just the problems at the end of chapter. I find them much harder than the corresponding examples on the same subject. Can anyone offer any suggestions for my predicament?


title: Power of the t-test under unequal sample sizes

text: I heard that the power of a t-test used with unequal samples is limited by the size of the smaller sample. Can I take this to mean that the power of a t-test with unequal sample sizes is equal to the power of a t-test used on equal sample sizes where $n$ equals the size of the smaller sample?


title: Clustering 1D data

text: I have a dataset, I want to create clusters on that data based on only one variable (there are no missing values). I want to create 3 clusters based on that one variable. Which clustering algorithm to use, k-means, EM, DBSCAN etc.? My main question is, in what circumstances should I use k-means over EM or EM over k-means? Thanks Ali


title: How to interpret the coefficients for binormal part in zero-inflated regression?

text: use `zeroinfl {pscl}` I specified offset only for the poisson part. The regression estimate zeroes using regressor `zeros`               zeroinfl(formula = response ~ 1 | zeros, offset = predictor, dist = "poisson")          Pearson residuals:           Min        1Q    Median        3Q       Max      -1.412429 -0.704399  0.003589  0.711585  5.667481           Count model coefficients (poisson with log link):                  Estimate Std. Error z value Pr(>|z|)     (Intercept) -0.002522   0.007094  -0.356    0.722          Zero-inflation model coefficients (binomial with logit link):                 Estimate Std. Error z value Pr(>|z|)     (Intercept) -10.7082    17.1139  -0.626    0.532     zeros        -0.8981     6.4082  -0.140    0.889          Number of iterations in BFGS optimization: 22      Log-likelihood: -1.702e+04 on 3 Df      Now the p-value for binomial coefficients is to test against null parameter = 0, which is equivalent to 50% (`logit(0.5) = 0`) Since it is insignificant, it means the regressor is not good enough. What I am interested to know is, does it also imply there is excessive zeroes or not at all if the coefficients are not significant? Because the null model is 50%, when the alternative hypothesis fails, should I fall back to the conservative claim and auto assume that the zeros are at least 50%? (which makes no sense at all)


title: Variance of sample Variance

text: In page 292 of Introduction to Mathematical Statistics by Hogg and Craig it is stated that in order for the variance of the sample variance, i.e. $\text{var}(S^2)$ to exist we need to assume that $E[X^4]<\infty$, that is the fourth (uncentered) moment of our RV is finite. My question is, how do these two connect? I tried to derive a formula for the variance of the sample variance that does not reveal the connection as far as I can see. Since we know that $S^2=\frac{1}{n-1} \sum_{i=1}^n \left( X_i-\bar{X} \right)^2$ is an unbiased estimator of $\sigma^2$, i.e. the population variance, what we need to find is $$E\left[ \left( S^2-\sigma^2 \right) \right]^2$$ which upon simplifying becomes $$E\left(S^4\right)-\sigma^4$$. From now on the algebra will get horrific though and I am not even sure whether we can proceed without stating the distribution of $X$. Is there perhaps another way to make sense of this statement? Thank you in advance.


title: How to predict and optimize a queue?

text: I am trying to determine if I can model my system as a M/M/1 queue and if so do the numbers I get from it help me at all. I can model my system like this: ![block diagram](http://i.stack.imgur.com/dLpOm.jpg) **System Description** **A.** I want to spec out the max number of clients I can have connect to my system. Each client can generate at max 40 messages/second that are placed in the buffer. If link B goes down it will buffer all messages. For the sake of my math I am going to assume that it is actually generating a fixed 40 messages/second. **B.** I am specifying that this link should be able to go down for 3 days, but normal operation spec this as a network connection that all clients will have to communicate on. This network connection may also have traffic unrelated to my system on it. The link is fast enough that it isn't a bottle neck at all. I have found that it can easily pass 500,000 messages in less then 10 seconds. **C.** Allows for the clients to dump any data as fast as link B will allow. This queue can be as big as needed. **D.** This isn't really a link per say, block E pulls directly from the queue. **E, F, and G.** This is where the fun comes. Initially every packet I received I placed a single SQL query to the database. I found that each query took about 5-15ms with a few taking 1ms and a few taking 600+ms. I don't have a model for this delay and suspect I will need one. I realized that this delay was going to drastically hurt my systems performance so I started looking into how to improve this performance. I found that if I grouped multiple messaged from C together and did a larger SQL query I was able to improve performance. When I grouped 10 messages per query I was observing 10-20ms, again with a few very short and a few very long. So 10-20ms/10messages gives me an effective 1-2ms per message instead of the 5-15ms I had before. **My Objectives**   * Determine how to model E, F, and G. I figured the easiest and most accurate method would be by experimentation, but I am not sure what data I should be looking to collect.   * Determine how many messages I should group together in order to minimize latency while maximizing the number of clients in stage A that I can have. I figured I could have a variable number of grouped messages based off of load, but I don't know how to come up with what this should look like.   * Finally be able to list a max number of clients I can have in stage A and say with some amount of certainty that their messages can be stored to the database in x time once it gets to stage C. Any help would be great. Thanks!


title: Modeling a gamma-mixture waiting model in BUGS

text: I'm analyzing a noisy time series where where the inter-event interval is known to follow a two-gamma mixture distribution. If there was a simple model that would generate that kind of thing, it would be pretty simple to implement into BUGS. But otherwise, I can't think of anything that wouldn't be prohibitively kludgey. Can anyone think of a model that induces mixture-gamma waiting times?


title: Methodology for data transform that best fits normal

text: SAS EM has a node that generates the transformation for a dataset that makes it best fit a normal distribution. I can't imagine the people at SAS are the only ones to figure out how to do this. Is there any documentation out there on methodologies that do the same? Or an open source implementation?


title: Informative censoring caused by cesarean section

text: I study labor and delivery as an epidemiologist. It is well established that a large fetus has a higher risk of causing maternal birth trauma. But a large baby is also likely to be delivered by cesarean section due to difficult labor. Those delivered by cesarean section will not have any risk for maternal birth trauma and, therefore, excluded from an analysis using trauma as the outcome. Thus, informative censoring is a real issue to estimate an unbiased relative risk, particularly when the cesarean rate is high. Are there any methods to adjust for such a bias?


title: Determining variance of meta-analysis log-response ratio generated from fitted curve

text: I'm working on a meta-analysis and have generated a quirky question for which I'm at a bit of a loss. The MA is for a large set of factorial experiments. Calculating the Log Response Ratio (LRR) and variance in said ratio for the experimental data is a cinch, and we're comparing the effects of one type of treatment to the other (and any interactions). However, my group is curious at examining the effect of HALF the level of one of the treatments (it's a continuous treatment, and we've been looking at the highest versus the lowest level of the treatment). We have fit nonlinear curves for each experiment that describe how the treatment effects the response over a wide range of treatment levels. We've got the coefficient error and residual error for each of these curve fits. And how we want to calculate log(full treatment) - log(half treatment) from the fitted curve for the this half-treatment log response ratio. Easy. But...how would we then calculate the variance in the half-treatment log response ratio? Would we use the SE estimates for the curve coefficients? The residual error? What would the sample size be for the variance calculation? Or is this unimportant? Thoughts?


title: Tips for training dropout neural network

text: I use NN for my mini project research, and I found out the newest trick for feed forward NN is using dropout for regularization instead of L1/L2 norm and rectified linear unit as an activation function. But when I tried it, I always got worse results compared to a standard NN with sigmoid / hyperbolic tangent activation function. Is there some rule of thumb or trick that we can use for training dropout ReLU NN? Thanks :D


title: Skin Detection with Gaussian

text: I'm reading this paper but i don't understand why the result of covariance and why the formula is different. I am using elliptical Gaussian joint probability density function (pdf) of this paper and the result covariance that i got is between <0-1>. Thank you


title: Simultaneous multiple perturbations in Markov chain Monte Carlo

text: I'm coding a **McMC** algorithm for geophysical applications. Using the Metropolis–Hastings scheme to accept/reject the proposed models is something that I thought I understood completely, but I don't. To be clear: normally you have a current model $m$, you perturb one parameter (in my case velocity $v_i$), and in this way you get a trial model $m'$. The acceptance probability for such a perturbation is: $$\alpha(m'|m)=min\left[ 1, \frac{p(m')}{p(m)} \cdot \frac{p(d|m')}{p(d|m)} \cdot \frac{q(m|m')}{q(m|m)} \cdot |J|\right]$$ Where $p(m)$ are the priors, $p(d|m)$ are the likelihoods, $q(m|m')$ are the proposals, and $J$ the Jacobian. So assuming that the priors are symmetric, and the Jacobian is 1, I remain with: $$\alpha(m'|m)=min\left[ 1,\frac{p(d|m')}{p(d|m)} \cdot \frac{q(m|m')}{q(m|m)} \right]$$ The likelihood ratio is not a problem, so now I want to deal only with the proposal ratio: Supposing that the way I make a perturbation is choosing one parameter, and perturb it with Gaussian probability, so that the new velocity value will be $$v'_i=v_i + u \cdot \sigma$$ ($u$ is normal distributed with mean = 0 and var. = 1)   This perturbation can be described by a proposal in the form: $$q(v'|v)= \frac{1}{\sigma \sqrt{2 \pi}} \cdot \exp \left[- \frac{(v'-v)^2}{2 \sigma ^2} \right]$$ It's easy to see that for this kind of perturbation the proposals are symmetric, so the acceptance probability is just the ratio of the likelihoods. **Questions:**   1. If instead of perturbing one parameter, I want to change 3 of them at the same time, what happens to my proposals, and so to my acceptance probability $\alpha$? How should I compute them? ($v_1, v_2, v_3$ are not correlated.)   2. If I want to include some deterministic constraints because I know that the 3 parameters I want to change are correlated, i.e. if I increase $v_1$ of $v_p$, then I want to decrease $v_2$ and $v_3$ of $v_p / 3$ (a sort of compensation)...how am I supposed to come up with a representation of the proposal?


title: R p-value half to SPSS rANOVA

text: I have done rANOVA in SPSS and then in R and got two different p-values from the same model. While SPSS gives 0.032, R gives 0.0162, which when rounded is a half. Other data are exactly the same. **1. Why?** Intuitively, one of them does a one-tailed, the other a two-tailed test. Other values are: Sum Sq `19.071` (SPSS) v. `24.14` (R) or F value `4.863` (SPSS) v. `6.156` (R). Other parts of the table are the same though, both `IV1 * IV2` and `residuals`. **2. Why? Which is _“more correct”_?** **3. How to get the SPSS-reported values from R?** SPSS approach:               Analyze –> General Linear Model –> Repeated Measures…      R approach:               summary(aov(DV ~ IV1 * IV2 + Error(subject/(IV1 * IV2)), data))      **Update to answer @PeterFlom and give additional details:** I am evaluating _the effect of an intervention (treatment: a lecture) on the ability to spot behavioral clues of emotions using visual stimuli_ using **pretest-posttest plus control group** experimental design. The experimental group is given the first set of stimuli, then there is intervention and the second set of stimuli follows. Naturally, the same without the intervention occurs in the case of the control group. **Two groups, two measurements each.** The **hypothesis** is that the intervention significantly raises the ability to spot behavioral clues of emotions. SPSS calculates the same interaction too. The data:                   subject   group   phase     value     1        A1     exp     pre        13     2        A2     exp     pre         7     .     35       B1    cont     pre         9     36       B2    cont     pre        14     .     57       A1     exp    post        11     58       A2     exp    post        12     .     91       B1    cont    post        13     92       B2    cont    post        12      Honestly, I am using it as is because a tutorial said so. I have got a very shallow understanding of the correct process in R. **Update in response to @ttnphns:** SPSS data:               group    pre    post       exp     13      11       exp      7      12      cont      9      13      cont     14      12      SPSS commands:               GET DATA        /TYPE=TXT        /FILE="F:\file.csv"        /DELCASE=LINE        /DELIMITERS=";"        /ARRANGEMENT=DELIMITED        /FIRSTCASE=2        /IMPORTCASE=ALL        /VARIABLES=        group A4        pre F2.0        post F2.0.      CACHE.      EXECUTE.      DATASET NAME DataSet1 WINDOW=FRONT.      GLM pre post BY group        /WSFACTOR=phase 2 Polynomial        /METHOD=SSTYPE(3)        /PLOT=PROFILE(phase*group)        /PRINT=DESCRIPTIVE ETASQ        /CRITERIA=ALPHA(.05)        /WSDESIGN=phase        /DESIGN=group.


title: What's good about I-projections?

text: There seems to be a large body of applied research where distribution q is picked to minimize KL(q,p) where p is empirical distribution. Are there theoretical reasons to prefer this estimator? For instance, a theoretical reason to prefer MLE for estimating mean of a one-dimensional Gaussian with mean between -1 and 1 because it's admissible, and also as sample size goes to infinity, it minimizes maximum risk and Bayes risk.


title: Overfitting in K-NN and Decision Trees?

text: To avoid over fitting for K-NN could you increase the value of K to reduce anomalous results etc. However, if the value of K is very large with respect to a sample, would this also incur in over fitting? To avoid over fitting in decision trees should you increase the depth of the tree to allow more rogue values to be classified correctly, whereas a lower depth could result in over fitting if there is not enough attributes to correctly classify the new instance.


title: PyMC beginner: how to actually sample from the fitted model

text: I'm trying a very simple model: fitting a Normal where I assume I know the precision, and I just want to find the mean. The code below seems to fit the Normal correctly. But after fitting, I want to sample from the model, ie generate new data which is similar to my `data` variable. I know I can use `trace("mean")` to get samples for the mean variable. But how can I get new samples from the model itself? I've looked at docs eg http://pymc- devs.github.io/pymc/database.html#accessing-sampled-data. I've also looked at quite a few examples, eg the mining disasters, and several from the Probabilistic Programming notebooks, and none mention this. I (more or less an MCMC beginner) expected that sampling from the fitted model was the whole point! What am I missing?               from pymc import *     data = np.array([-1, 0, 4, 0, 2, -2, 1, 0, 0, 2, 1, -3, -1, 0, 0, 1, 0, 1])     mean = Uniform("mean", -4, 4)     precision = 2.0**-2     obs = Normal("obs", mean, precision, value=data, observed=True)     model = Model( {"mean": mean, "obs": obs})     mcmc = MCMC(model)     mcmc.sample(10000, 1000, 1)     # I can get samples for the "mean" variable     mean_samples = mcmc.trace("mean")[:]     hist(mean_samples)     # but how can I do the equivalent of mcmc.trace("obs")?


title: Does a cointegration model exist for irregularly spaced time series?

text: It isn't clear to me how to calculate cointegration with irregular time series (ideally using the Johansen test with VECM). My initial thought would be to regularize the series and interpolate missing values, although that may bias the estimation. Is there any literature on this subject?


title: How to make a time series stationary?

text: Besides taking differences, what are other techniques for making a non- stationary time series, stationary? Ordinarily one refers to a series as "integrated of order p" if it can be made stationary through a lag operator $(1-L)^P X_t$.


title: Regression on Negatively Skewed Data

text: I am trying to explore the relationship between a customer satisfaction rating and a response time (in days). I have been using linear, logistic and quantile regression, to no avail. I believe the issues lies with the distribution of the data. 85% of respondents give a 9 or 10 rating (regardless of response time) and 70% of the data lies in the 9/10 rating and the 0-4 day response time. Since the data is so negatively skewed, I'm looking for some advice on whether or not to normalize the data and if regression is the right method to define the relationship between the variables (and if not, what would be a better method). I've attached a summary of the data distributions percents. The range for the rating is 1-10 inclusive and the response rate is continuous, but we limit to 10 days because of low volume for 10+ days. Any assistance is appreciated. ![Distribution by Percent](http://i.stack.imgur.com/cy1it.png)


title: When is it incorrect to compute factor scores by summing (or averaging) raw variable scores?

text: I understand that a problem will be inevitable if the variables have different scales of measurement. My inclination is to think that even if the variables have the same scale of measurement it would be _incorrect_ to simply sum (or to simply average) the raw scores of the variables that are "enough loaded" by the factor, since doing so would weight each variable equally, while in fact the factor loads some of those variables more and some of them less. Is my reasoning correct? Are there any other circumstances in which it is incorrect to compute factor scores just by summing (or averaging) the raw scores of the factor's variables? The factor scores are being computed with a view to them being used as predictors in a subsequent regression analysis.


title: Similarity between two sets of points as a level of a path overlay

text: Say I have two sets, A and B, with a different amount of X,Y points Is there a way to calculate how similar the data is? So something like this would be ~65%?, because they partial overlay.               |                                  B8     |                                       |                                 B7     |                               [A7]     |B1                        B6[A6]     |   B2               B5[A5]     |        B3   B4[A4]     |        [A3]     |     |     [A2]     |        |     |   [A1]     |____________________________________      And data like the following would have a low percentage even though they are similar. they are not overlaying.               |                B6     |             B5     |          B4     |       B3     |    B2     | B1     |                         [A7]     |                     [A6]     |                 [A5]     |             [A4]     |         [A3]     |     [A2]     | [A1]     |____________________________________      Also this data would have a low percentage, even though they "intersect"               |  B1     |                     [A6]     |       B2     |                 [A5]     |           B3     |             [A4]     |               B4     |         [A3]     |                  B5     |     [A2]     |                      B6     | [A1]     |____________________________________      And this **data would be almost 100%, because their path overlays**               |     |                             B5     |                            [A5]     |     |     |                            B4     |                         [A4]     |                  [A3]B3     |         [A2]B2     |     B1     |  [A1]     |____________________________________      Another example of the issue is. Say I record a set [A] of points (latitude, and longitude) while I'm driving. I drive from my house on to the highway and to work. A friend also records a set [B] of points from his house to the **SAME** strip of highway I drove, to his office. The data would look like the **first graph**. Is there a way to calculate how similar the data (geolocation) is? Is it possible to calculate a percentage of similiarity between the two sets of points? I've looked into the following methods, but can't get it to work:   * Spearman/Pearson's Correlation: I think this is the closest thing I'm looking for, but it only compares two variables, instead of sets of points   * Linear approx.: The problem with this is that the data is definitely NOT linear   * Canonical Correlations Analysis: Not sure how to use it...


title: Implementing ANCOVA via residuals

text: While it is often mis-used, I'd like to provide users of `ezANOVA()` (from the ez package for R) with the ability to specify covariates for ANCOVA. Since `ezANOVA()` uses a combination of `stats::aov()` and `car::ANOVA()`, I think the simplest way for me to achieve ANCOVA within `ezANOVA()` is to simply compute two anovas, one that includes just the covariates as predictors from which I would obtain residuals that would be used as the predicted data for a second anova that uses the regular non-covariate predictors. Does anything look inherently wrong about this approach? (P.S. I thought I'd also include some evaluations of ANCOVA's assumptions, like checking that the covariates don't correlate with any of the non- covariate predictors, and checking that the covariates' effects on the DV don't interact with those of the non-covariate predictors.)


title: How to do piecewise linear regression with multiple unknown knots?

text: Are there any packages to do piecewise linear regression, which can detect the multiple knots automatically? Thanks. When I use the strucchange package. I could not detect the change points. I have no idea how it detects the change points. From the plots, I could see there are several points I want it could help me to pick them out. Could anyone give an example here?


title: Connection between bootstrap (both parametric and non-parametric) and maximum likelihood

text: Can anyone show me what the connection between bootstrap methods and maximum likelihood is? I was told that as the number of bootstraps goes to infinity, the statistics like confidence intervals, mean, variance and stability are same values that maximum likelihood would estimate.


title: Is the sample variance a useful measure for non-normal data?

text: Does it ever make sense to compare the variance for two sets of data, neither of which are even approximately normal (e.g. bimodal)?