title: Nested vs random effect

text: I'm having a little trouble analysing a data set from the following experiment: we applied 4 treatments (`control`, `concentration1`, `conc2`, `conc3`) to 5 plants each for a grand total of 20 plants. We repeated the experiment 3 times independently (`bloc`). When I visualised the data, I saw that there was important variability between the bloc on the response variable, so I decided to control the bloc effect as following:               mymodel <- aov(y ~ treatment + Error(bloc), data = mydata)      But I am not sure it is what I am trying to do. What is the difference between the above and:               mymodel <- aov(y ~ treatment + Error(treatment/bloc), data = mydata)     mymodel <- aov(y ~ treatment + treatment/bloc, data = mydata)      And which is the most appropriate for my analysis? I read the "nested effect model" part in the R Book, but it didn't help me. Any help is welcome. Thanks in advance EDIT: We had a total of 60 plants, 20 for each bloc


title: Constructing confidence intervals for predictive model

text: I'm constructing a model to predict the weight of a species of insect given a set of other variables. The plot below shows the performance of my model using a set of test data where the true weight of the insects is known. The x-axis is the true weight of the insect and the y-value is the associated error of my model—the absolute value of the predicted weight - true weight: ![enter image description here](http://i.stack.imgur.com/S2DWZ.jpg) From this visual, you can see that I have many insects with relatively low weights. In these cases, my model has a relatively lower predictive error. In contrast, I have relatively few insects which are heavy and the error associated with these predictions displays much more variance along with higher error. Given this model and the test data, I'd like to find a way to construct confidence intervals for new predictions. For example, if my model predicts a given insect is relatively heavy, the confidence intervals around this prediction would be large. My questions is, how can I do this? A linear model seems inappropriate for this data since most of the points are clustered near the origin. I'm at a loss for how I can construct confidence intervals for my predictions.


title: Choosing the number of components of a linear model via cross-validation

text: **Background** I have a series of Principal Components which I'm regressing on some data (it's a Blind Source Separation problem; I assume that my data is a linear combination of signals which I obtained via PCA), and I have to decide on the number of components for the model; I chose to do this via $K$-fold Cross- Validation. **The problem** I'm having problems because the number of components chosen via CV is highly sensitive to the choice of $K$, the number of folds of the CV procedure. I have $n=40$ datapoints and a maximum of 8 components to regress the data on. Here's an example of a $10$-fold Cross-Validation procedure: ![enter image description here](http://i.stack.imgur.com/DTjgX.png) In this case, the minimum CV error occurs using 8 components, but the value 1-SE away from this value is 3, so I might be tempted to use 3 components in this case. However, in this case the validation folds are too small, and that explains the high variance on the CV-error. If I instead use, e.g., $5$-fold CV, then the value 1-SE away from the minimum is 7 components; this number oscillates for each value of $K$ and so does the number of components where the global minimum is. **My approach** After looking at several CV plots for different $K$'s, I noted that there is always a local minimum at 4 components, followed by an increase in the CV- error on components $>4$, followed by a global minimum on components $7-8$. Considering that I'm using principal components, it seems reasonable to me to use this number as the 'optimal' number of components because this is the first minimum I encounter counting from the lowest to the highest number of components (considering that the number of components are ordered in terms of their eigenvalues). I find that this reasoning, however, is by no means 'clean' or an 'elegant' solution/explaination to the problem at hand; it just seems intuitive to me...do you have any advice on expanding this solution or proposing another one?


title: Data normalization and classification

text: I hope this clearly states the problem I have in hand. Here goes: I've trained a neural network with one initial data set that was normalized in order to guarantee an equal participation of each variable in the learning process. Once the the neural network is trained I'll have new sets of data coming to be classified. Q1: How should I proceed in terms of normalizing a new coming data set? Is there any chance that a different mean and standard deviation will make this new data set significantly different from the one I originally used to train my classifier? Can this generate "wrong" results? (of course I'm considering the new income samples to come from the same source as the original training set). Q2: Am I supposed to follow the same normalization procedure if instead of just classifying I'm first clustering the training data set and them classifying new income samples (into the clusters/groups found by my clustering algorithm)?


title: Pairwise Correlations--Multiple Comparison Correction

text: I have a dataset with two sets of variables defining each sample. I am performing pairwise correlations between each variable in the first set with every variable in the second set. I essentially have $n$ samples. For each sample, I have two sets of variables $[a_1, ..., a_{m1}]$ and $[b_1, ..., b_{m2}]$ that define that particular sample. Correlating all the variables in $a$ with those in $b$ gives me a final correlation matrix of size $m_1 \times m_2$. In addition, some of the variables within set $a$ may be correlated with each other and some of the variables within set $b$ may also be correlated with each other, but I am not assessing those correlations. I'm trying to figure out:   1. if I need to correct for multiple comparisons, and    2. if I do, what method I should use.  I have tried Bonferroni, but with the large number of comparisons I get an extremely small adjusted p-value.


title: Pros/Cons of recoding ordinal/nominal variables to target mean for logistic regression?

text: Say I have an independent variable with the following relationship to the binary dependent variable, DV:                ___________________________________________________________________________     |verx_s                          |  # Recs |    % Recs   |  # DV   | DV Rt  |     |________________________________|_________|_____________|_________|________|     |0                               |   75,700|      6.4467%|      941| 1.243% |     |1                               |  277,129|     23.6009%|    1,471| .5308% |     |2                               |   51,662|      4.3996%|      219| .4239% |     |3                               |  769,737|     65.5526%|    2,269| .2948% |     |All                             |1,174,228|    100.0000%|    4,900| .4173% |     |________________________________|_________|_____________|_________|________|      It's common practice at my company to recode the values of verx_s to the value of DV Rt and treat it as a continuous variable when modeling logistic regression. Confidence intervals are not important. All we care about is whether the model validations on an out-of-time sample. Is there anything inherently wrong with taking this shortcut? It should also be mentioned that in most cases the independent variable is crafted in such a way that it makes intuitive sense to our customers. Therefore the _ordering_ of the target mean is important. Hence, why we can't use simple dummy vars.


title: How to derive the $\alpha$ for the Pareto rule

text: Suppose we have the CDF for the Pareto Distribution given by: $$ P(X \leq x) = 1-\left(\frac{x_m}{x}\right)^\alpha \;\;\;\;\;\;\;\;\;\; x \geq x_m$$ What is the intuitive way to find the alpha for which the 80/20 rule holds?


title: Matlab or R: relationship between categorical vs continuous variable

text: I have a six groups=[1:6] ACTid and a resulting continuous variable observed for each of the members is the group: Y               ACTid(1:10)  Y(1:10,:)          ans =              3.0000   36.8791         2.0000   71.5823         6.0000  104.3075         3.0000   33.4862         1.0000    1.5298         5.0000    0.1063         6.0000  157.3075         4.0000   11.4117         2.0000    3.8024         1.0000   52.3496      What's the best way to check the relationship of the groups with the continuous variable? The Continuous variable is derived from this model: **UPDATE** My model is of following multiplicative form: > Let us simplify notation by allowing the number of germs (cfu) units picked > up by nurses' hands to be called $Y$. **ACTid** determines **n** in, $i=1:n > $ > > \begin{equation} Y_{i}=\displaystyle \lambda_i V_i A_i-\beta_iY_{i-1} > \end{equation} > > Where values are sampled via Monte-Carlo sampling for each run for: > $h=lognormal~(1.5,0.1)$, $\lambda=\Gamma(15,3)$, $\beta=$empirical non- > negative. $A=lognormal~(7,1.9)$ and $V=$empirical non-negative. The choice > of these parameters is further dependent on the type of care $ACTid=1:6$;


title: Simulation of Fisher's exact test is underestimating power

text: Summary: My simulations don't match my power calculation. (This question has become rather long, so you might not read to the bottom. I think it's a sample size problem.) I'm running some simulations to help design a study, and I'm running into something that's confusing me. There are two groups, treated and untreated, and the outcome is crashing. We have an estimated crash rate of 0.1 in the untreated group, and think we can bring this down to 0.08. A small (and non-feasible) sample size is 100 per group. So I run power.prop.test() to estimate the power.               power.prop.test(p1=0.1, p2=0.08, n=100)          > power.prop.test(p1=0.1, p2=0.08, n=100)                                                                       Two-sample comparison of proportions power calculation                         n = 100                  p1 = 0.1                  p2 = 0.08           sig.level = 0.05               power = 0.07122973         alternative = two.sided           NOTE: n is number in *each* group       Then I ran a simulation, because I wanted to explore how often we'd make the wrong decision as to the best treatment.               library(Exact)          #Create a data frame called d, populate it with the numbers above.     set.seed(12345)     nTreated <- 100     nUntreated <- 100     probTreated <- 0.1     probUntreated <- 0.08          d <- data.frame(id = 1:10000)     d$nTreated <- nTreated     d$nUntreated <- nUntreated     d$probTreated <- probTreated      d$probUntreated <- probUntreated               #Generate some random results using rbinom()     d$treatedCrashes <- apply(cbind(d$nTreated, d$probTreated), 1,                            function(x)  sum(rbinom(x[1], 1, x[2])))          d$untreatedCrashes <- apply(cbind(d$nUntreated, d$probUntreated), 1,                            function(x)  sum(rbinom(x[1], 1,                                                    x[2])))               #Do fisher's exact test on each replication:     d$fisher <- apply(cbind(d$nTreated - d$treatedCrashes,                              d$treatedCrashes,                             d$nUntreated - d$untreatedCrashes,                              d$untreatedCrashes), 1,                                function(x)  fisher.test(matrix(x,                                                              nrow=2))$p.value)     #test power     mean(d$fisher < 0.05)      And I get 4.8% power, which is lower than the power.prop.test function said, AND is less than 0.05 - which seems kind of wrong.               >     mean(d$fisher < 0.05)     [1] 0.0478      Is this about the small sample approximation? Is it a stupid coding error? (I don't think it is, but I've often been wrong about that before). Is it something I haven't thought of? In response to the suggestion that it's because the Fisher's exact test is conditioned on the margins, I reran the model with Barnard's test (in the Exact library). (But I reduced to 1000 replications, as this took 40 minutes).               d$exact <- apply(cbind(d$nTreated - d$treatedCrashes,                               d$treatedCrashes,                              d$nUntreated - d$untreatedCrashes,                               d$untreatedCrashes), 1,                             function(x)  exact.test(matrix(x, nrow=2), to.plot=F, cond.row=T)$p.value)     d$exact <- lapply(d$exact, function(x) x[1][[1]])          mean(d$exact < 0.05)                > mean(d$exact < 0.05)                           [1] 0.049                      I have pretty much the same power. However, I also ran an exact test for power using the power.exact.test() function, also in the Exact library, which gives a very similar level of power:               > power.exact.test(p1=0.1, p2=0.08, n1=100, n2=100,simulation=T, nsim=1000, method="Boschloo")     $power     [1] 0.045          $alternative     [1] "two.sided"          $method     [1] "Boschloo"      Which makes me think that it's a sample size problem.


title: Is a diffuse prior always weak?

text: There are countless examples of diffuse priors being used to 'allow the data to speak.' However, what if one's past experience leads you to be skeptical of new data, without necessarily having a prior as to where the bulk of the density is located. Is there some way to have a diffuse prior that nevertheless strongly influences the posterior? **Edit for clarification:** I'm thinking of a situation in which you want to capture some of the "wisdom" of an experienced researcher who knows nothing about a particular domain but has priors about the amount of evidence required before he/she accepts a finding. Something like the Ioannidis story. So it's diffuse because you don't have any reason to believe one value over another, but it's informative in that you want to down-weight the data.


title: Terminology used in normal probability plots

text: I'm vaguely familiar with normal probability plots (non-statistician who knows enough to get by). The graph in this article: What are the two dotted curves called, what do they signify, and how are they calculated? ![](http://www.public.iastate.edu/~wrstephe/stat403/plotnorm.gif) **edit:** My question is software-independent. I just want to know what side bands used in a typical normal probability plot are called and how they are defined. * * * I don't have a good answer but after stabbing around in the dark on Google I finally found this lecture notes on quantile-quantile plots: https://www.stat.auckland.ac.nz/~ihaka/787/lectures-quantiles2-handouts.pdf


title: Probability calculation in survival analysis

text: I have 13 similar independent units and I need to calculate the probability at least 6 of them to survive for time less than 1.10. From a table with the Kaplan-Meir estimates I get the following;               Time          1.05614 , 1.31581         number at risk  19   ,   18         number failed 1,1         survival Prob. 0.474771 , 0.448395      I have omitted the previous and the next values because `Time=1.10` is between those two values and `P(T<1.10)=F(1.10)=F(1.05614)=1-survival prob.=1-0.474771= 0.525229` Which is the probability for 1 unit to survive for less than 1.10 hours? I don't know to continue from that point.


title: Testing for differences, would I use ANOVA?

text: I have four age groups all measuring different heights: How would I test for difference between the groups at 0.01 level of significance? Would I use ANOVA? Also - how could I establish which group is different from the others at the 0.01 level of significance?


title: What is the name of this equation to combine probabilities?

text: I am doing some work that involves combining the probabilities of multiple independent, non-disjoint events into an overall probability of any event occuring. We are using an equation like this: $P(A_1 \cup A_2 \cup \dots \cup A_N) = 1 - \prod_{i=1}^N{(1 - P(A_i))}$ Is there a name for this equation? I can't seem to find where it comes from in probability theory.


title: How to form a predictive model in R?

text: I have two data sets whose structure is like this: **DATA SET 1:**                       month_year  sales      [1,]  "Jan 2000"  "30000"      [2,]  "Feb 2000"  "12364"      [3,]  "Mar 2000"  "37485"      [4,]  "Apr 2000"  "2000"      [5,]  "Jun 2000"  "7573"               .     .      .               .     .      .      **DATA SET 2:**                         month_year    profit          [1,] "Jan 2000" "84737476"          [2,] "Jan 2000" "39450334"          [3,] "Jan 2000" "48384943"          [4,] "Feb 2000" "12345678"          [5,] "Feb 2000" "49595340"          [6,] "Jan 2001" "36769493"                   .     .      .                   .     .      .      As it can be seen the **first data set has one sales value for each month of each year while in second data set I have`n`(say 100 but not constant all the time) number of values for each month of each year.** Now I want to develop **a predictive model** which can provide **predictions for profit for future coming months**. What is the best approach to develop such a predictive model? I searched online and found that we can form **a regression** or correlation **fit** model and based on that can do predictions. Which way is better- Regression or correlation? And **how to form a regression fit for such data sets** where one data set has multiple data points for each month of each year while the other has a single data point for each month of each year? **NOTE:** I don't know whether it is the right approach or not but if we can get the regression/correlation value for like the month of Jan,Feb,Mar...and so on and then based on sales value we can do prediction of profit in general. I am new to this and hence I can lack in having certain concepts clear and thus request for some guidance for this issue.


title: Autoregressive distributed lag models ADL(p,q): seeking "how to" preferably in Matlab (Stata/R/Python/C# etc)

text: Could anyone provide me the details of how to determine the lag order of the distributed lags for an $\text{ADL}(p,q)$ model in Matlab or another statistical package (and very much preferably in combination with the autoregression lags)? Full working examples with model selection criteria ($\text{AIC}$ and $\text{BIC}$) seem to be available on the Matlab website for $\text{VAR}$ models, $\text{ARMA}$ models etc. but I can't find one for the $\text{ADL}(p,q)$ model. I would not have a clue how to rewrite those models to $\text{ADL}(p,q)$ myself but I have a vague feeling that such a thing would be possible. In the end I want to automate this analysis by first checking the lag orders $p$, and $q$ and then using these numbers automatically to create the regressions out of this. So basically I'm looking for a fully working example. (I want to skip part of adding and deleting regressors by hand as much as possible to get a quick idea of the distributed lags of several assets).


title: Good reference on analysis of dependent samples

text: Is there a good reference (book if exists) about mathematical statistics of dependent random variables? Something along the lines of Casella and Berger's textbook, but for dependent samples. I'm not interested in a specific application, just theory, of course assuming independence sometimes seems to be good enough, but a lot of of algorithms in machine learning assume iid samples, however i would like to pose what if the data is not independent what is the body of theory for this?, or whether there is such a book or resource ever written about this.


title: How does one go about clustering data?

text: (I have updated the question following conversation with @whuber in the comments). My case is as follows: I have around one thousand row vectors of dimension $1 \times 8$. These row vectors are filled with known and non-negative numbers (they are all rational numbers), that can also be zero. In some positions I have $0/1$ variables, but not in all. In principle, they may exist identical vectors. Through utility maximization theory, each row vector is linked to a corresponding profile of preferences (of $8$ dimensions). Each position in the vector represents an observable aspect of a product characteristic synthesized with the product selling price, and the corresponding position in the preferences vector represents the utility effect/hedonic price of this characteristic. What I need to do is to somehow group these vectors to represent preferences profiles much fewer than one thousand. I could then treat each sub-group as a collection of realizations from the same multivariate random vector, and proceed to characterize the distributions. Strictly speaking, if two vectors differ even in one of the eight positions, and even in the slightest, they represent a different preferences profile.   But I am willing to accept that, for each position, a preference profile can produce a range of values. That is, I accept that "indifference curves are thick" as we say in Economics, or that preferences have some randomness in them. From what I understand this is what we call "data clustering", and there are various methods to do this, but this is new territory for me, so I do not trust me to search in past CV posts. This is related to research, so please treat it like homework: names of methods and some literature to get me started will really be enough. If any additional info is needed, please ask.


title: Removing the intercept term in a dynamic regression justified?

text: I'm trying to build a cross-sectional prediction model (dynamic panel) with the following form: (including a LDV) $Y_{i,t+1}=a+bY_{i,t}+cX_{i,t}+e_{i,t+1}$ As the sample contains for example countries that differ greatly in terms of size and as all variables are in absolute values (e.g. absolute GDP in USD) the intercept term severely biases my forecasts. As my intercepts are often extremely negative, the forecast for a small country with a low gdp will be negative. I' am wondering now whether you can improve the model (in terms of higher forecast accuracy) by simply removing the intercept. Doing so implies that we go through the origin but I don't see how that should not fit my model. What are other consequences? I've read that R^2 can be significantly higher but R^2 is not the best way to evaluate a model anyway and as I include a lagged DV I should not pay too much attention on R^2 anyway.


title: Learning association between multiple effects and causal outcomes

text: We have some data represented as follows: on each row there are 10 observable effects, and 3 non-observable causes out of 8 that were guessed by scholars. Effects are real numbers, and causes are discrete categories. We want to train a system that would accept observed effects, and guess three most likely causes. So how to represent "3 out of 8 causes", and what would be the best approach and how to use such a learning algorithm? Edit: Observables are characteristics of cracks in a structure. Causes are things like load, wind, temperature. Example sample is like 1.2 2.3 45 6.7 5.5 12 3.4 1.1 5.6 2.3 load weather temperature It would be better if we could determine the rank of causes. E.g. most important is load, then weather, etc. Because the sample is gathered this way.


title: How to calculate the log likelihood in HMM from the output of the forward algorithm in R?

text: How can I get the log likelihood or the probability value from the forward algorithm for an observed sequence? For example, when I executed the forward algorithm in `R` on a sequence of length 3 for a trained Hidden Markov Model of 3 states, it gave me the probability values for each observed symbol on each state: $$ \begin{matrix} ~ &O1 & O2 & O3 \\\ S1 & -143.2 & -500.0 & -231.5 \\\ S2 & -212.2 & -231.4 & -200.0 \\\ S3 & -112.4 & -115.6 & -118.5 \end{matrix} $$ In the above matrix, rows are states and columns are observed symbols. Now how do I calculate the final log likelihood value from this matrix?


title: Problems with Poisson Regression

text: I have been studying a linguistic construction (let’s call it `C`) in a language `L1` trying to individuate several factors (`F1`, `F2`, ..., `Fn`) which can influence / trigger the presence of `C` in `L1`. All factors are categorical variables with several levels. I ran basic significance tests (Goodness of fit) which revealed a high level of significance for each factors involved. That said, I would now like to measure the interaction of all factors / all levels. As for the method, I think that Poisson Regression is what I am looking for since “C” has no levels and can be represented in terms of frequencies. I thus tried Poisson regression in R but I have soon run into several problems. I have two questions in particular:   1. The data frame I used looked like that (here you can find only a tiny part of it as an example):              Subj (Factor1)   NP (F2)    New (F3)    EXP_CON (F4)      To get the frequencies of all variable combinations, I used the function `ftable` and then I went on manually computing the frequencies. Needless to say, it was a painful procedure and it took a lot of time! Does R provide an easier and quicker way to get the frequencies for all the interacting factors / levels (possibly eliminating automatically the interactions whose frequency is zero)? In other words, I would like to obtain something like that:              Sub (F1)     NP (F2)      EXP_CON (F3)         New (F4)        (Freq) 6        2. After that, I tried to run a Poisson regression model but the result was puzzling. Here is an example of what I got:              Syn_FunOther:Focus_CCIMP_CON:IS_CCGiven        NA         NA      NA       NA     Syn_FunSub:Focus_CCIMP_CON:IS_CCGiven          NA         NA      NA       NA     Syn_FunOther:Focus_CCNO_CON:IS_CCGiven         NA         NA      NA       NA     Syn_FunSub:Focus_CCNO_CON:IS_CCGiven           NA         NA      NA       NA     Syn_FunOther:Focus_CCNOV:IS_CCGiven            NA         NA      NA       NA      I do not really understand. Why are there so many `NA` values?


title: How to write a loop in R to select multiple regression model and validate it?

text: I would like to run a loop in R. I have never done this before, so I would be very grateful for your help !   1. I have a sample set: 25 objects. I would like to draw 1 object from it and use it as a test set for my future external validation. The remaining 24 objects I would like to use as a training set (to select a model). I would like to repeat this process until all 25 objects are used as a test set.    2. For each of the training sets I would like to run the following code: $\,$               library(leaps)     forward <- regsubsets(Y ~.,data = training, method = "forward", nbest=1)     backward <- regsubsets(Y ~.,data = training, method = "backward", nbest=1)     stepwise <- regsubsets(Y ~., data = training, method = "seqrep", nbest=1)     exhaustive <- regsubsets(Y ~.,data = training, method = "forward", nbest=1)     summary(forward)     summary(backward)     summary(stepwise)      summary(exhaustive)      I would like R programme to select the best model (with the highest adjusted $R^2$) using each of the selection methods, so there are 4 final best models (e.g. the best model selected with forward selection, the best model selected with backward selection and so on...). Afterwards I would like to perform internal cross validation of all four selected models and choose one of the four which has the lowest average mean squared error (MSE). I used to do it using the code below:               library(DAAG)      val.daag<-CVlm(df=training, m=1, form.lm=formula(Y ~ X1+X2+X3))      val.daag<-CVlm(df=training, m=1, form.lm=formula(Y ~ X1+X2+X4))      val.daag<-CVlm(df=training, m=1, form.lm=formula(Y ~ X3+X4+X5))      val.daag<-CVlm(df=training, m=1, form.lm=formula(Y ~ X4+X5+X7))       For the best selected model (the lowest MSE) I would like to perform an external validation on 1 object left on the site at the beginning of the study (please refer to point 1.). 3\. And loop again using different training and test set .... I hope that you could help me with this. If you have any suggestions how to select the best model and perform validation more efficiently, I would be happy to hear about that. Thank you !


title: Name for the Bayesian posterior probability that a regression coefficient is larger than zero

text: I have the following logistic regression: $$ \text{logit} (y) = \beta_0 + \beta_1\, x $$ from which I can estimate the following posterior probability (using a Bayesian approach): $$ P(\beta_1>0\,|\,\text{Data}). $$ Is there a particular name for that probability (something like Bayesian one-sided p-value)?


title: What would be a suitable way to present a prevalence rate?

text: Suppose we say that the prevalence rate of some disease is $1$ in $50$. Does this make sense? Or should be convert this into a $ \%$? That is, should we say that the prevalence rate is $2 \%$?


title: Computing prediction intervals for logistic regression

text: I would like to understand how to generate _prediction intervals_ for logistic regression estimates. I was advised to follow the procedures in Collett's _Modelling Binary Data_ , 2nd Ed p.98-99. After implementing this procedure and comparing it to R's `predict.glm`, I actually think this book is showing the procedure for computing _confidence intervals_ , not prediction intervals. Implementation of the procedure from Collett, with a comparison to `predict.glm`, is shown below. I would like to know: how do I go from here to producing a prediction interval instead of a confidence interval?               #Derived from Collett 'Modelling Binary Data' 2nd Edition p.98-99     #Need reproducible "random" numbers.     seed <- 67          num.students <- 1000     which.student <- 1          #Generate data frame with made-up data from students:     set.seed(seed) #reset seed     v1 <- rbinom(num.students,1,0.7)     v2 <- rnorm(length(v1),0.7,0.3)     v3 <- rpois(length(v1),1)          #Create df representing students     students <- data.frame(         intercept = rep(1,length(v1)),         outcome = v1,         score1 = v2,         score2 = v3     )     print(head(students))          predict.and.append <- function(input){         #Create a vanilla logistic model as a function of score1 and score2         data.model <- glm(outcome ~ score1 + score2, data=input, family=binomial)              #Calculate predictions and SE.fit with the R package's internal method         # These are in logits.         predictions <- as.data.frame(predict(data.model, se.fit=TRUE, type='link'))              predictions$actual <- input$outcome         predictions$lower <- plogis(predictions$fit - 1.96 * predictions$se.fit)         predictions$prediction <- plogis(predictions$fit)         predictions$upper <- plogis(predictions$fit + 1.96 * predictions$se.fit)                   return (list(data.model, predictions))     }          output <- predict.and.append(students)          data.model <- output[[1]]          #summary(data.model)          #Export vcov matrix      model.vcov <- vcov(data.model)          # Now our goal is to reproduce 'predictions' and the se.fit manually using the vcov matrix     this.student.predictors <- as.matrix(students[which.student,c(1,3,4)])          #Prediction:     this.student.prediction <- sum(this.student.predictors * coef(data.model))     square.student <- t(this.student.predictors) %*% this.student.predictors     se.student <- sqrt(sum(model.vcov * square.student))          manual.prediction <- data.frame(lower = plogis(this.student.prediction - 1.96*se.student),          prediction = plogis(this.student.prediction),          upper = plogis(this.student.prediction + 1.96*se.student))          print("Data preview:")     print(head(students))     print(paste("Point estimate of the outcome probability for student", which.student,"(2.5%, point prediction, 97.5%) by Collett's procedure:"))     manual.prediction     print(paste("Point estimate of the outcome probability for student", which.student,"(2.5%, point prediction, 97.5%) by R's predict.glm:"))         print(output[[2]][which.student,c('lower','prediction','upper')])


title: Is cross validation for validating a model or for selecting best model in different kinds of models?

text: I am confused about the concept of cross validation and its usage. As I read about cross validation before, it is a way of validating a model. I did cross validation in my project (developing different regression models on a dataset, model validation and finally choosing best model). Cross validation just give me a model but no statistical criteria that shows ability of the model. By the way each time that I run cross validation (in R program) the result is different because train dataset is changed. For selecting the best model, I calculate AIC for the model obtained by cross validation but now I think it is wrong. Because first, train data used for each model is different, and second, even for a certain model, AIC will change by repeating cross validation. Lately I read that cross validation is used for selecting the best model! All these confused me about utility of cross validation and the way of interpreting the result. Could you please help me on figuring out all these?


title: Does pooled variance correct for/protect from unequal variance when calculating effect size?

text: This may be a lame question, but I got stuck and can't get my head around it. I am running a gene expression analysis, comparing $\sim 10,000$ genes between two groups, $n=6$ samples per group. My pipeline goes like this:   1. I check whether the two groups have equal variances using Levene's test for each gene.   2. For those genes that do have equal variances, I run a moderated t-test (as in limma), as it assumes equal variances.   3. For the rest of the genes I run Mann-Whitney U. (Before you tell me about multiple test correction, I do calculate FDR later on by permuting the population labels on samples and re-running the above tests). Now, I want to calculate effect size for all the genes. I am using Cohen's d with pooled variance (see Wikipedia), but I wonder whether I can do it for both groups of genes (ie ones with and without equal variance between groups). I think I can't, but it's a guess really (I think I can only run it for the equal variance genes). If I can't, is there a more appropriate way to do this?


title: Standardising data within-subjects when analysing distributions

text: Analysing behavioural data, Spivey, Grosjean and Knoblich, (2005) wanted to show that their results came from a unimodal distribution, rather than from averaging over subpopulations in a bimodal distribution. Their data was structured $42(participants) \times 2(conditions-\ control\ \&\ cohort) \times 16 (trials\ per\ condition)$, and to examine the distributions they standardised scores within-subjets, and then pooled them: > With too few trials within a participant to provide an adequate measure of > the unimodality or bimodality of the distribution of these trajectory > deflections, the values for both cohort and control conditions were together > converted into z scores within a participant and then pooled across > participants. > > [...] > > We also z-scored, within each participant, the area under the trajectory > separately for cohort and control trials, and we then pooled across > subjects. These pooled distributions were then analysed for bimodality, using the bimodality coefficient (not cited in the paper), among other measures, from which they concluded that their data was from a unimodal distribution. This standardise-pool-analyse approach has become the standard in the mouse- tracking literature inspired by this paper. However, given that no previous work is cited to justify this method, and we know that analyses reported in PNAS aren't always above reproach, I'd like to ask the opinions of greater statistical minds than mine. **Is this approach justified? Is it an objective method of assessing bimodality in the data, or does standardising make data appear more unimodal than it actually is?** I should also note that Freeman and Dale (2012), in a technical report on statistical tools for assessing bimodality, standardise their behavioural data in the same way, without discussing the matter: > Using the z distributions of the trajectories’ AUC values [...]


title: How to model two variables with a cumulative exponential relationship?

text: Suppose I have a response variable(`y`), which is normally distributed, and `y` is generally changed with `x` in a way like a cumulative exponential distribution, which means `y` may reach a plateau at some turning point of `x`, such like that:               set.seed(123)     df<-data.frame(y=pexp(x)+runif(length(x),0,0.2),x=1:901)     df<-df[sample(1:300,900,T),]     hist(df$y)     plot(df$x,df$y,xlab="x",ylab="y")      ![y](http://i.stack.imgur.com/TMd7j.png) ![xy](http://i.stack.imgur.com/Rw3Df.png) My question is what's the best way to model this sort data? Spline, polynomial or something others?


title: What is a mathematical way to define a point on a scatter plot as an outlier?

text: I have a graph and there are two points that could be two potential outliers. I'm trying to create a polynomial line of best fit with undefined order. I believe I could use a >2 standard deviations exclusion rule, but I'm sure how exactly this is applied. Do I determine the qualifiying standard deviation by using all the points initially? That's what I did, but I have two points that qualify for exclusion. Do I remove both at once or one at a time? If I remove only the one with the greatest deviation from the line of best fit, I can recalculate the line and it will give me a new line of best fit of new polynomial order. If this new fit is used, the other point which was previously an outlier is no longer an outlier. What is the correct procedure?


title: SPSS Repeated measures with different pre-post measures

text: I have the following data: 20 participants, each one going through 2 types of treatments (A & B) at 3 different sessions (times). For each treatment session, i am collecting 6 "normal" measurements, having just one value per session. In addition, for each treatment session, i have 2 different tests with pre and post values, for which i am interested in the difference between pre and post values. I would like to check which of the measures or the test differences is different when getting treatment A compared to treatment B. Table for example (the full data has Measure 1..6):               | Participant Number | Treatment | Session | Measure1 | Measure2 | Test1_pre | Test1_post | Test2_pre | Test2_post |     |--------------------|-----------|---------|----------|----------|-----------|------------|-----------|------------|     | 1                  | A         | 1       |          |          |           |            |           |            |     | 1                  | A         | 2       |          |          |           |            |           |            |     | 1                  | A         | 3       |          |          |           |            |           |            |     | 1                  | B         | 1       |          |          |           |            |           |            |     | 1                  | B         | 2       |          |          |           |            |           |            |     | 1                  | B         | 3       |          |          |           |            |           |            |      Etc.. for each participant. This is the first time i am using SPSS, and i am trying to understand how to do such analysis. So far, i have managed to use a repeated measures analysis, defining the treatment and session factors and all measurements, but i still could not understand how to set the pre-post tests values so that the difference between pre and post values for each test is considered. It seems to me like the pre values should be defined as covariates, but i cannot see a way to connect each pre value to its correct post value. I could put a new column in my data with the calculated difference for each test, but i understood that this is less recommended since it neglects the actual values of the tests. I would appreciate any help or guidance. Thank you very much!


title: Proper analyses for 2x2 contingency tables

text: So I have multiple questions about the different types of analyses that can be used on a 2x2 contingency table and when different analyses should be used. To open, I'll note that this stems from my efforts to find the appropriate analyses to run for research using categorical data in a 2x2 contingency table. In my case I have two conditions: control and experimental and within the experimental group there are two items of interest (say, Item A and Item B). In all cases I'm investigating whether an outcome occurs or not: How many experimental participants show the outcome on Item A? How many experimental participants show the outcome on Item B? How many control participants show the outcome? So in addition to comparing Item A with the Control group (A vs. C) and Item B with the Control group (B vs. C), I also want to compare Item A with Item B (A vs. B). My understanding is that for between-subjects comparisons (e.g., A vs. C & B vs. C) a Pearson's $\chi^{2}$ test is traditionally used and that for within- subjects comparisons (e.g., A vs. B) that McNemar's test is a common choice. **Question 1:** is this an accurate understanding of the common choices for statistical tests using 2x2 contingency tables? * * * My follow-up questions, and where most of my confusion lies, are related to what to do when multiple _observed_ cell counts are low (less than 10). I've explored numerous threads here, the Wikipedia pages for different topics (Yates's correction, Fisher's exact, etc.), a few journal articles, and descriptions in R-packages, but I have found conflicting views. Some sources indicate that you should use Yates's correction, others say it's out-dated because, in addition to being overly conservative, it was only preferred because the other options were computationally intense. Others say that you should use Fisher's exact test, but I've seen additional information suggesting that Fisher's exact test is only appropriate if you know both the row and column totals and that other options (Boschloo, Barnard, etc.) are generally more powerful. However, which of those to use is also not clear. Finally, even others suggest that the Pearson's $\chi^{2}$ test is appropriate as long as the _expected_ cell counts are not low. I've seen some indication that these arguments occur for McNemar's test as well (e.g., continuity correction & McNemar's exact test). **Question 2:** Is the assumption regarding cell counts about _expected_ cell counts or _observed_ cell counts? If the _observed_ cell counts are low, does that matter to the type of test that should be used? **Question 3:** What _is_ an exact test and when should exact tests be used? * * * In addition to answers to these questions, I would also greatly appreciate any suggestions as to books or articles I can read to learn more about this. Part of my struggle with all of this has been just general difficult finding clear information about these topics and how they relate to one another.


title: Clustering variables with outliers

text: I am performing a clustering analysis in SAS and some of the variables that I am trying to cluster contain outliers. I've tried to transform the data (log and/or standardize them) but didn't quite work out. So, for example, let's say I came down to 9 clusters, then one or two clusters will have just one value in them. Deleting outliers is not optional due to the nature of my work. I want to assign a set of values into groups so that values in the same clusters are more similar (in some sense or another) to each other than to those in other clusters however, if I decrease the number of clusters and basically merge the outliers group with other groups I'm afraid that the merged group won't be homogeneous. Thus, to sum up my questions, how can I deal with outliers for cluster analysis? Thanks in advance!


title: Hidden Markov models and expectation maximization algorithm

text: Can somebody clarify how hidden Markov models are related to expectation maximization? I have gone through many links but couldn't come up with a clear view. Thanks!


title: Interpreting exp(b)

text: With respect to (binary) logistic regression and a categorical IV, State, if a given exp(b) is .08 and I'm using indicator coding such that my reference category for the variable, State, is, say, "Wyoming," how do I interpret .08? Presumably I contrast the .08 with the "mean" of Wyoming (with "mean" here being the proportionality, or number of positive instances, of Wyoming)...So if Wyoming = .06, do I multiply .08 * .06 = .0048? That doesn't quite feel right...


title: Distributions similar to Normal distribution

text: Can anyone give me examples of distributions symmetric around the mean but that are more densely concentrated around the mean than the Normal distribution for the same mean and variance? Thanks


title: Three questions about support vector regression: feature pre-processing, time-series issues, and marginal accuracy contribution of each feature

text: I am familiar with using support vector machines as the base classifier in a Python re-implementation of Poselets. But I am new to the use of support vector regression (not machines for classification) for time-series data. In my data set, I have a vector of target variables $y_{i,t}$ and a matrix of predictor variables (features) $X_{i,t}$ where $i$ denotes the individual in the population (rows of observations), and I have these cross-sections of targets and predictors for many time periods $t$. I use the notation $x_{i,t}$ to specify the entire set of features (entire row) for individual $i$ at time $t$. The notation $x^{j}_{t}$ would denote the column with the $j$-th feature across all individuals at time $t$. So $X_{i,t}^{j}$ is the $(i,j)$ entry of the feature matrix at time $t$. In my case, all of the predictors are continuous variables, though I can imagine some categorical predictors that could be useful. We can also assume the time frequency is fixed for all target and predictor variables. My question is threefold: (1) Is there a set of generally accepted best practices for pre-processing the predictors (scoring, normalizing, rescaling, etc.) to provide more sensible inputs to the support vector regression step? (E.g. can you point me to a reference for this. I've already looked in the Smola and Scholkopf tutorial, the Bishop book, and Duda, Hart, and Stork, but found no useful guides or discussions of practical pre-processing concerns). (2) Are there similar references or best-practices for handling time-series inputs to support vector regression? For example, in a linear regression setting, sometimes people are worried about cross-sectional correlation of error terms or predictor variables. This is treated by taking averages across the time dimension, generally under assumptions of no autocorrelation, and (e.g. Fama-MacBeth regression) these can be shown to be unbiased estimators of the coefficients even in the presence of cross-sectional correlation. I've never encountered places where support vector methods were described as sensitive to cross-sectional correlation between predictors or error terms. Is this a concern? How is it handled in support vector modeling? At one extreme end, I could imagine doing the following: For a given time period $t^{*}$, not only use the predictor variables that are formulated at the same time period, $X_{i,t^{*}}$ but also augment this set of feature vectors with lagged versions of the variables, for whatever trailing time periods I might believe are relevant. And in the scoring step of (1) I could perhaps apply a decaying weight function so that newer observations are given more weight. But to me this sounds like a bad idea. For one, the predictors $X_{i,t}$ are generally chosen because of domain-specific knowledge, intuition, or some other method of getting prior information about what might be a successful predictor. Just dumping all lagged versions of the same thing feels a bit too much like data mining and I worry it would lead to overfitting big time. Secondly, even with support vector methods, there is still a goal of having a parsimonious model, so surely I would not want degrees of freedom for a huge set of lagged versions of predictor variables. Given all of this, what are some examples of accepted or practically useful support vector regression treatments of time series data? How do they address these problems? (Again, a reference is desirable. I'm not looking for someone to write a graduate thesis as an answer here.) (3) For the model being fitted, there are definitely going to be some comparisons drawn with more classical models, like plain OLS or GLS. In those models, there are two very nice features. First, the degrees of freedom is very transparent because there are generally no hyperparameters to tune and no pre-processing steps beyond linear scoring steps on the input variables. Second, after the regression functions is fitted, it is straightforward to attribute the accuracy of the classifier to the different linear components of the model. Are there analogues for these with support vector regression? When using kernel functions, like a Gaussian RBF, to transform the inputs, is it fair to say that the only extra degrees of freedom introduced are whatever hyperparameters govern the kernel's functional form? That feels a bit wrong to me, because you're effectively allowing yourself to explore a whole space of unvetted transformations of the data, which isn't really captured in just the functional form of the kernel function. How can you fairly penalize a model like this for having much more freedom to overfit the data with non-linear transformations? And lastly, are there analogous ways to decompose the accuracy of a fitted support vector model? My first thought here was that you would somehow need to measure something like the mutual information of the classifier trained inclusive of predictor $x^{j}$ (the $j$-th column of the feature vector inputs) with the classifier trained without that one predictor. And this measure would in some sense be the marginal contribution of accuracy when that predictor is used. But then this would require fitting at least $K+1$ models (the full model and $K$ different leave-one-predictor-out models, if there are $K$ predictors). And if you include lots of lagged variables, $K$ might be very large, and you'd be performing some type of cross-validation at each step. It seems like getting that kind of marginal contribution statistic is vastly more expensive for a support vector model. Are there references showing different summary statistics that capture the same kind of marginal contribution statistics?


title: The correct probability distribution / way to identify large deviations in a set of daily changes to portfolio value

text: I am working on a report which is being sent through to end users that should flag to them any "large changes" in the day-to-day values for the past 30 days. These values are day-to-day differences which we assume to have a mean of zero. To shed a bit more light on the data, the values are proportional returns to a portfolio. They are used to measure the change in the riskiness of different aspects a portfolio which is constantly managed, so the riskiness should always stay at "round about the same level" based upon the investing strategy associated with that given portfolio. That being said, the risk will fluctuate a bit from day-to-day, but it gets re-balanced and so the change in daily risk should fluctuate around zero - positively for roughly half the time & negatively for roughly half the time. This analysis is to simply point out any large daily deviations so they can be examined. So, assume we have the following data:               Day:    Change from previous day:     1       -40     2       30     3       15     4       12     5       -34     6       -2     ...     30      12      And, as I stated earlier, they don't care about the direction of the change, just to flag any day where the change is of a "statistically significant" magnitude. I'm curious as to how to calculate the statistics for this properly / what distribution would best be used to assign probability levels for this set of data. What has been proposed is to, firstly, look at the square of the changes (rather than the original values since all we care about is magnitude):   1. To calculate the std dev as $\sqrt{ \dfrac{\sum \left( x_i - 0 \right)^2 }{29}}$ - In other words **assume** a zero expected value   2. To Calculate the std dev as $\sqrt{ \dfrac{\sum \left( x_i - \bar{x} \right)^2 }{29}}$ - In other words, use the mean we would calculate for this sample And then to take each day's value and divide it by this std_dev to calculate it's magnitude. I have quite a few problems with a few things here, but especially idea #2 since I believe the mean SHOULD BE ZERO. Generally speaking, what is the correct way to analyze deltas with an expected value of zero and what distribution would you use for a 30-sample data-set to find statistically significant values? Note that historical data is available / can be easily gathered. If we assumed the daily changes to be ROUGHLY normally distributed, would the square of them (only measuring deviation from zero, regardless of whether it's positive or negative) be distributed via Chi-Square? ... Or could a folded- Normal work if we took their absolute value? I'm just looking for something to show how change magnitude might be distributed... Nobody would take any issue with assuming the changes to be normally distributed if necessary... I hope this makes sense, but, if not, please let me know where I could clarify more. THANKS!!! (I originally asked this question here and was recommended to ask it here - Hopefully someone here can give me a good response)


title: Classifying handwritten digits using PCA

text: > Classify handwritten digits using PCA. Use 200 digits for the train phase > and 20 for the test. I have no idea how PCA works as a classification method. I've learned to use it as a dimension reduction method where we subtract the original data from its mean, then we calculate the covariance matrix, eigenvalues and eigenvectors. From there, we're able to choose principal components and ignore the rest. How should I classify a bunch of handwritten digits? How to distinguish data from different classes? Or does it mean something totally different, that I must use PCA for feature extraction purposes and use a classification method afterwards?


title: Implications of doing a confirmatory factor analysis with a correlation matrix as input instead of a variance-covariance matrix?

text: Is this possible, and if so what are the implications of doing things one way rather than another. Is one approach generally preferable? So far I have only been taught to use a variance-covariance matrix.


title: Standard error of a sum of two Poisson variables

text: I have 2 Poisson variables, I know mean and standard error for each. How do you calculate the Standard error of the sum?


title: Nonparametric Time Series Forecasting

text: I am trying to understand how Kernel Density Estimation (KDE) or (nonparametric) Quantile Regression can be used to forecast values given historical observations. For example, consider the following nonparametric regression model $Y_t = g(X_t) + u_t$ where $X_t = (Y_{t - 1}, \ldots, Y_{t - q})$. In the following equation (one- step-ahead forecasting) $\hat{Y}_{n + t,1} \equiv \hat{\text{E}}(Y_{n + t}|X_{n + t}) = \dfrac{\sum_{j = q + t}^{n + t - 1}Y_jK_h(X_j - X_{n + t})}{\sum_{j = q + t}^{n + t - 1}K_h(X_j - X_{n + t})}$ where $K_h(X_j - X_{n + t}) = \prod_{s = 1}^q h_s^{-1}k((Y_{j - s} - Y_{n + t - s})/h_s)$, it doesn't seem to me that the order of the observations is taken into account when evaluating the kernel function at previous observations. Shouldn't it matter? Otherwise, it seems to me that observations in previous time periods are just treated as i.i.d. data (no time structure). Similarly to the prediction of future quantiles if I were to use the function npqreg in the R np package for time series data. Thank you in advance for any comments.


title: Evaluation indexes hypothesis for clustering

text: I read on the cluster analysis page of wikipedia: > For example, k-means clustering can only find convex clusters, and many > evaluation indexes assume convex clusters. On a data set with non-convex > clusters neither the use of k-means, nor of an evaluation criterion that > assumes convexity, is sound. I can't see how an evaluation index assumes convex clusters. Can anyone illustrate this idea?


title: What best practices should I follow when preparing plots?

text: I usually make my own idiosyncratic choices when preparing plots. However, I wonder if there are any best practices for generating plots. Note: Rob's comment to an answer to this question is very relevant here.


title: How to tell if something happened in a data set which monitors a value over time

text: I have a data set where a series of measurements are being taken each week. In general the data set shows a +/- 1mm change each week with a mean measurement staying at about 0mm. In plotting the data this week it appears that some noticeable movement has occurred at two points and looking back at the data set, it is also possible that movement occurred last week as well. What is the best way of looking at this data set to see how likely it is that the movements that have been seen are real movements rather than just some effect caused by the natural tolerance in the readings. **Edit** Some more information on the data set. Measurements have been taken at 39 locations which should behave in a similar way although only some of the points may show signs of movement. At each point the readings have now been taken 10 times on a bi-weekly basis and up until the most recent set of readings the measurements were between -1mm and 1mm. The measurements can only be taken with mm accuracy so we only receive results to the nearest mm. The results for one of the points showing a movement is 0mm, 1mm, 0mm, -1mm, -1mm, 0mm, -1mm, -1mm, 1mm, 3mm. We are not looking for statistically significant information, just an indicator of what might have occurred. The reason is that if a measurement reaches 5mm in a subsequent week we have a problem and we'd like to be forewarned that this might occur.


title: Clustering time series of measurements in R

text: I have a dataframe consisting time series of measurements taken every hour for 366 days or a year. Below is shown a sample of hourly measurements for the first two weeks. I want to cluster days with similar profile of hourly measurements. By using R, how can I group days with similar profiles? Thanks.                      Day1 Day2 Day3 Day4 Day5 Day6 Day7 Day8 Day9 Day10 Day11 Day12 Day13 Day14     Hour1    11   14   20   23   51   13   18   15    8    25    13    24     8    48     Hour2     2   13   21   19   21   17    8    6   11    21    16    25    11    18     Hour3     5    6    7   13   14    6    8    9    7    15    12    11     2    11     Hour4     6    4    6   18   11    2   10    4    5    20    10    10     3     8     Hour5     4   12    5   13    7   13    7    9   17    15    22     9     1     4     Hour6    33   33   25   20   13   46   29   42   27    22    32    29    30    10     Hour7    65   72   70   19   12   72   73   66   75    21    80    74    62     9     Hour8    79  105   78   45   30   87  107   75   86    47    91    82    76    27     Hour9    96   64   75   48   34   89   96   91   96    50   101    79    93    31     Hour10   73   91   70   80   42  103   63   82   79    82    84    74    70    39     Hour11   61   65   71   82   48   64   65   70   66    84    71    75    58    45     Hour12   54   59   99   71   62   68   85   65   58    73    63   103    51    59     Hour13   85   69  106   95   83   57   67   74   73    97    78   110    82    80     Hour14   78   75  108   72   82   58   80   78   70    74    75   112    75    79     Hour15   65   93   95   89   66   72   91   75   78    91    83    99    62    63     Hour16   95   98   92   88   94  103  108   88   88    90    93    96    92    91     Hour17   97  118   95   95   85  125  113  106  106    97   111    99    94    82     Hour18   86   86  114  106   83  104   98   84  109   108   114   118    83    80     Hour19  101   95  105   90   71   84   89   85  103    92   108   109    98    68     Hour20   60   55  102   56   55   46   72   61   71    58    76   106    57    52     Hour21   55   51   54   57   65   42   58   42   36    59    41    58    52    62     Hour22   45   34   23   30   33   51   45   56   37    32    42    27    42    30     Hour23   44   31   36   43   34   29   30   44    0    45     5    40    41    31     Hour24   25   17   30   40   14   24   17   15   28    42    33    34    22    11


title: Optimization in R vs Python, constrained, unconstrained and automatic differentiation?

text: I am an economics/stat guy who uses quite a bit of optimization (maximum likelihood, simulated maximum likelihood), constrained optimization (mathematical programming w/ equilibrium conditions), dynamic programming, etc. I am wondering how Python compares to R for optimization. In my opinion (feel free to argue this and show me if I am wrong), R is not great for optimization in that it does not have automatic differentiation nor great solvers for higher dimensional problems. Programs like AMPL (can hook up to any world class solver) or even Matlab with KNITRO solver won't make you explicitly enter a horribly complex jacobian, hessian, nor sparsity pattern, but R (even IPOPTR) makes you do this. As such, I find R poor for optimization. Is Python any better? Does it have the automatic differentiation feature? Can it link up to world-class solvers? Any information on this would be great. Otherwise, I fear I will have to do my work in AMPL or CPLEX, which I do not know, and learning yet another language isn't my highest priority unless I need to.


title: What does correlated data mean and how can I visualize it with scatterplot?

text: I'm playing around with the Abalone dataset in R and following along with this article. The dataset has 8 variables that are taken into account to predict the number of `rings`. To find the pairwise correlation, the blog post does this:               as.matrix(cor(na.omit(abalone[,-1])))      and comes to conclusion that the data is heavily correlated. My question is how do they come to this conclusion? What information should I be looking for to come to this conclusion? Here is the code               > aburl = 'http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data'     > abnames = c('sex','length','diameter','height','weight.w','weight.s','weight.v','weight.sh','rings')     > abalone = read.table(aburl, header = F , sep = ',', col.names = abnames)     > as.matrix(cor(na.omit(abalone[,-1])))                       length  diameter    height  weight.w  weight.s  weight.v weight.sh     rings     length    1.0000000 0.9868116 0.8275536 0.9252612 0.8979137 0.9030177 0.8977056 0.5567196     diameter  0.9868116 1.0000000 0.8336837 0.9254521 0.8931625 0.8997244 0.9053298 0.5746599     height    0.8275536 0.8336837 1.0000000 0.8192208 0.7749723 0.7983193 0.8173380 0.5574673     weight.w  0.9252612 0.9254521 0.8192208 1.0000000 0.9694055 0.9663751 0.9553554 0.5403897     weight.s  0.8979137 0.8931625 0.7749723 0.9694055 1.0000000 0.9319613 0.8826171 0.4208837     weight.v  0.9030177 0.8997244 0.7983193 0.9663751 0.9319613 1.0000000 0.9076563 0.5038192     weight.sh 0.8977056 0.9053298 0.8173380 0.9553554 0.8826171 0.9076563 1.0000000 0.6275740     rings     0.5567196 0.5746599 0.5574673 0.5403897 0.4208837 0.5038192 0.6275740 1.0000000     > pairs(abalone[,-1]      ![enter image description here](http://i.stack.imgur.com/xD98j.png) **Questions**   * What is the problem when the data is heavily correlated?   * How can one come to conclusion that the data is heavily correlated by looking at the matrix or by looking at the scatterplot? Is it because the scatterplot shows linear lines in almost every row-column?


title: Bayesian Spatial Prediction: What is f?

text: I'm working my way through Gelfand, A.E. & Li Zhu, B.P.C. (2001). On the change of support problem in spatio-temporal data. Biostatistics, 2:1, 31-45. I'm stuck at: ![enter image description here](http://i.stack.imgur.com/09HUl.png) This is absolutely the first mention of f anywhere in the paper. Can someone help me to identify it?   * Y = spatial process    * s = measured locations    * Beta, Theta = presumed coefficients for mean function and covariance function    * H = pairwise (?) covariance function <\-- not too clear about that


title: Random Forest proof notation

text: I am having a bit of difficulty understanding the notation in equation (1) on page 4 of the following paper: https://escholarship.org/uc/item/35x3v9t4#page-4 Specifically, what do $E_{X,Y}$ and $E_\theta$ mean? It seems like it is a saying the $PE$ of the forest tends toward the $PE$ of a single tree as $k \to \infty$, but that doesn't make sense.


title: Does least squares regression imply normality of errors?

text: For the linear model $$y_i=\beta_0 +\sum_{k=1}^{n}\beta_k x_{ik} + \epsilon_i$$ the parameter estimates are the same for the maximum likelihood method and the least square method (minimizing $\sum_{i=1}^n\epsilon_i^2$) if the errors are assumed iid normal. So I am thinking that if we use least squares, we are perhaps implicitly assigning normality to the errors. Is this correct - can it be shown that the least square method does this (without appealing to the maximum likelihood solution)?


title: How do I choose between a simple and a mixed effect logistic regression?

text: I have a list of predictor variables to put in to a logistic regression model. How I know that should I do a simple logistic regression (using `glm` function in R) or a mixed effect logistic regression (using `glmer` function in R) analysis?


title: Power law fitting of data with uncertainties

text: I need to fit data points on a power law and each one of these carries an uncertainty. I've been using Python, more precisely scipy.optimize.curve_fit to get the job done but I don't know how to handle the uncertainties with it. I thought about using a linear fit in log-log scale but it seems less precise compared to fitting directly on a power-law (not taking care of the uncertainties, just testing both on generated data and added noise). I found that but it's not really helping me. I wouldn't mind using R to do the job if needed.


title: Sum of Squares Constraint Regularisation

text: Suppose I have a function $g_0\in L_2(\mathbb{R})$ such that we observe $(X_i,Y_i)$, $i=1,...,n$ such that $$ \begin{align*} Y_i & = g_0(X_i) + V_i \end{align*} $$ We wish to estimate $g_0$, and since $g_0\in L_2(\mathbb{R})$, we can approximate $g_0$ by using an orthonormal series $\\{h_k\\}_{k\in\mathbb{N}}$. So for each sample size $n$, we approximate $g_0$ by considering the truncated series $g_{p_n}(x) = \sum_{k=1}^{p_n}\alpha_k h_k(x)$, such that $\alpha_k = \langle g,h_k\rangle$, and my estimator $\hat g = \sum_{k=1}^{p_n}\hat\alpha_kh_k$ will be defined as $$ \begin{align*} \hat{\boldsymbol \alpha} = \arg\min_{\boldsymbol{\alpha}\in\mathbb{R^{p_n}}} \sum_{i=1}^{n}\left[Y_i - \sum_{k=1}^{p_n}\alpha_kh_k\right]^2\\\ \end{align*} $$ The point is I wish to constrain the fitted coefficients $\hat{\boldsymbol\alpha}$ such that for each entry of $\hat{\boldsymbol \alpha}$, $\alpha_k$, I would like $|\alpha_k| \leq C$ for some $C$ strictly larger than $\|g_0\|_{\infty}$. I understand that if we were to Heuristically argue this change, all we would do is threshold everything to be less than or equal to in magnitude to $C$, but how does this constraint change the standard "least squares" solution of $\hat{\boldsymbol \alpha}$ in closed form?


title: Cox proportional hazard model and interpretation of coefficients when higher case interaction is involved

text: Here is the summary-output of the Coxph-model I used (I used R and the output is based on the best final model i.e. all significant explanatory variables and their interactions are included):                coxph(formula = Y ~ LT + Food + Temp2 + LT:Food + LT:Temp2 +      Food:Temp2 + LT:Food:Temp2) # Y<-Surv(Time,Status==1)      n= 555                              coef         exp(coef)          se(coef)      z           Pr(>|z|)          LT             9.302e+02      Inf             2.822e+02    3.297        0.000979 ***      Food           3.397e+03      Inf             1.023e+03    3.321        0.000896 ***      Temp2          5.016e+03      Inf             1.522e+03    3.296        0.000979 ***      LT:Food        -2.250e+02    1.950e-98        6.807e+01    -3.305       0.000949 ***      LT:Temp2       -3.327e+02    3.352e-145       1.013e+02    -3.284       0.001022 **      Food:Temp2     -1.212e+03    0.000e+00        3.666e+02    -3.307       0.000942 ***      LT:Food:Temp2   8.046e+01    8.815e+34        2.442e+01     3.295       0.000986 ***      ---      Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Rsquare= 0.123   (max possible= 0.858 )      Likelihood ratio test= 72.91  on 7 df,   p=3.811e-13      Wald test            = 55.79  on 7 df,   p=1.042e-09      Score (logrank) test = 78.57  on 7 df,   p=2.687e-14       Question is: How to interpret coefficient and exp(coef) values in this case, as they are very large values? Also 3-case interaction is involved, which confuses the interpretation more. All the examples concerning Coxph-model I have found so far online have been really simple regarding the intercation terms (which have always turned out to be unsignificant) and also coefficient-values (=hazard rates) and exponentials of these (=hazard ratios) have been pretty small and "easy to handle" numbers, e.g. coefficient = 1.73 -> exp(coef) = 5.64. BUT mine are way bigger numbers as you can see from the summary output (above). And because they are so large vaues, they almost seem to not make any sense. It seems to be a bit ridiculous to think that survival is e.g. 8.815e+34 (hazard ratio taken from the interaction LT:Food:Temp2) times lower when the interaction increases by one unit (?). Actually I don't know how to interpret this 3-case interaction either. Does it mean that when all of the variables in interaction increase by one unit, the survival decreases by the certain amount (told by the exp(coef)-value)? Would be great if somebody can help me out here. :) Below is the part of my data sheet I used for the cox-analysis. Here you can see, that I have used many times same explanatory variabe value (i.e. LT, Food, and Temp2) for several "Time, Status response variable". These explanatory variable values are already the mean values of these variables (due to the field-work setup out in the nature, it was not possible to get individual explanatory variable value for each observed response individual, hence the mean values used already in this phase), and this would answer to suggestion 1 (?) (see the first answer). Suggestion 2 (see the 1st answer): I am using R, and I am not yet super god in it. :) Thus, if I use function predict(cox.model,type="expected"), I get a huge amount of different values and have no clue to which explanatory variable they are referring to and in which order. Or is it possible to highlight certain interaction term in predict function? I am not sure if I am making myself very clear here. Suggestion 3 (see 1st answer): in the part of the data sheet below, one can see the units of different explanatory variables. They are all different and include decimals. Can this have something to do with the cox outcome? Part of the data sheet:               Time (days)     Status      LT(h) Food (portions per day) Temp2 (ºC)     28                0         14.42        4.46             3.049     22                0         14.42        4.46             3.049     9                 1         14.42        4.46             3.049     24                0         15.33        4.45             2.595     24                0         15.33        4.45             2.595     19                1         15.33        4.45             2.595      Cheers, Unna


title: Multivariate analysis and PCA

text: I observe y as a function of x. y is a scalar, x is a vector of dimension n. I have k observations of y and x. Let $$X = \begin{bmatrix}x_{1,1} & \cdots & x_{1,n}\\\\\vdots & \vdots & \vdots \\\ x_{k,1} & \cdots & x_{k,n}\end{bmatrix}$$ $$Y = \begin{bmatrix}y_1 & x_{1,1} & \cdots & x_{1,n}\\\\\vdots &\vdots & \vdots & \vdots \\\ y_k & x_{k,1} & \cdots & x_{k,n}\end{bmatrix}$$ The variables in x are supposed to be roughly independent (I can do a PCA to verify that, what would be the implications of an non-identity diagonalized matrix?), and I want to quantify which variables in x explain y the best. The idea is to look at the principal axes in the PCA of Y, do you think this is the right way? How do I interpret the results then?


title: solving an equation with more than two variable algebrically

text: I want to solve an equation with Lagrange multiplier. I know elimination method but I need to know this method. I wonder if anyone let me know how to estimate $\mu$ and 3 unknown variables ($\hat A$) algebraically . $n.\hat\mu$ + $n_1\hat A_1$ + $n_2\hat A_2$ + $n_3 \hat A_3$ + $0 \lambda$ = $y_{..}$   $n_1\hat\mu$ + $n_1\hat A_1$ + $1 \lambda$ = $y_{1.}$   $n_2\hat\mu$ + $n_2\hat A_2$ + $1 \lambda$ = $y_{2.}$   $n_3\hat\mu$ + $n_3 \hat A_3$ + $1 \lambda$ = $y_{3.}$   $0\hat\mu$ + $1\hat A_1$ + $1\hat A_2$ + $1 \hat A_3$ + $0 \lambda$ = $0$


title: Correlation plots with missing values

text: I have 9 vectors of numeric data. I want to make correlation plots for these 9 vectors. However some vectors are of a different length to others - some have up to 240 observations while others have only 159 observations. R can only make correlation plots for vectors of the same length. To make the vectors the same length, I simply need to use:               length(x) = length(y)      The problem is that this chops of up to 81 observations from some vectors. Is there any way round this, other than simply chopping off observations? If I do need to chop off observations, how do I determine which observations to chop off?


title: How to normalize two time series for comparison?

text: I have two time series a and b, which I want to compare. Due to their range difference I normalize them first.               a(i), b(i) are natural numbers for i=1,...,N      two different normalizations: ( mean and std_dev both refer to the whole time series ) 1) `a'(i) := a(i) / mean( a )` _goal_ : mean( a' ) = 1 2) `a'(i) := [ a(i) - mean( a ) ] / std_dev( a )` _goal_ : usual normalization what confuses me is how do the meanings after those transformations differ? does the first transformation make any sense at all?


title: Bayesian bandits with delayed rewards

text: I looked into the topic of bayesian bandits in order to create a simple testing tool for headline optimizations. UCB1 seemed easy enough until I discovered that there is probably a problem with the delayed rewards in the real life setup. I am in no way a mathematician. Maybe anyone could point me in the right direction on how to handle this? Thanks Valentin


title: Correcting sample bias

text: I am studying the correlation of two observed variables (call them $A$ and $B$). the underlying distribution for $A$ is symmetric around $0$ (for sure), however in my sample I have $411$ observations where it is negative and $470$ where it is positive due to sample bias. Is it a correct way to remove the bias to simply filter out $470-411 = 59$ random positive observations from variable $A$ and the corresponding observations of variable $B$ before carrying on my analysis? This would rebalance the distribution of observations for $A$. Ok, let's take an absurd example to make my point. Suppose you have a process following a perfect normal distribution centered at 0. A student draws 10000 samples one at the time and writes down the result. However, for reasons known only to her, she likes to write all the positive results on sheets of red paper and all the negatives on sheets of violet paper. If she can write 100 numbers on each sheet, at the end she will have _roughly_ 50 red sheets and 50 violet sheets. She takes back all the sheet to her office desk to carry on the analysis but something terrible happens: a mutant paper-eating creature attacks her and eats only 23 of her violet sheets (the creature despises red paper). Now she has still around 50 red sheets, ie. around 5000 positive numbers, but only around 27 violet sheets, ie. 2700 negative numbers. She cannot repeat the experiment. What can she do? Can she choose randomly 23 red sheets and throw them away to rebalance the samples?


title: K-fold cross validation for neural networks: separate validation set also needed?

text: When using k-fold cross validation in a neural network, do you also need a separate validation set? Or is the use of the k-fold on its own good enough to minimise the possibility of over-training?


title: Cluster analysis on time series samples

text: In the follow-up to this Ways to understand 2-dimensional time-series data I'm working on 2D time series data where two attributes are depth and temperature. When I plotted depth-vs-temp curve and saw its variation over time, the fluctuation occurs at few places only. Lets say temperature is dependent on depth and it varies at few depths over time. Would doing cluster analysis on each time sample and comparing it with the next one be good? The point is to statistically identify changes across timestamps. And then decide which are significant timestamp samples. What models/papers should I study to get insights out of the data?


title: factor analysis for given data with help of matlab

text: suppose that we have following data ![enter image description here](http://i.stack.imgur.com/Nd9Ux.png) i have done covariance matrix and eigenvalue decomposition               >> means=mean(X);     >> centered=X-repmat(means,10,1);     >> covariance=(centered'*centered)/(9);          >> [V,D]=eig(covariance);     >> [e,i]=sort(diag(D),'descend');     >> sorted=V(:,i);          e          e =             1.0e+03 *              3.8861         0.4605         0.1162         0.0465          >> sorted          sorted =              0.5080    0.5042    0.3747   -0.5893         0.0090   -0.6091    0.7930   -0.0092         0.8560   -0.3490   -0.2747    0.2647         0.0955    0.5030    0.3941    0.7633      i have calculated also percentage distribution               (e./sum(e))*100          ans =             86.1796        10.2125         2.5776         1.0303      for instance i want to choose first two,because in sum they have more then 96% of total variance,now please how can i continue for factor analysis?matlab commands also visual help will be very good,according to this answer Steps done in factor analysis compared to steps done in PCA i was confused,so please for given concrete numerical example,how can i continue,not just definitions,i need concrete numerical steps.thansk in advance EDITED : because i have choose first two eigenvector,i have done PCA also               factors=sorted(:,1:2)          factors =              0.5080    0.5042         0.0090   -0.6091         0.8560   -0.3490         0.0955    0.5030          >> PCA=X*factors          PCA =            299.5739   24.8602       150.1380   65.5846       183.9319   36.7606       156.7749   11.6325        58.1864   17.0585       236.6723   35.3584       176.0291   68.7008       208.4237   45.3761       190.6039   33.7278       206.2131    3.4430      how can i calculate loadings?how can i do orthogonal rotations?


title: Lag order in ur.df or adf test [R]

text: 1. I want to test stationarity using adf test or ur.df function on R ur.df(y, type = c("none", "drift", "trend"), lags = 1,selectlags = c("Fixed", "AIC", "BIC")) My question is when using adf.test the function selects the lag order on its own using some criterion but how we use ur.df if i am unsure what lag order to use?   1. When using adf.test i am facing a peculiar error timeseries -> stationary , p<0.01 LN(timeseries) -> non-stationary, p=0.13 Is there anyway to explain this? Thanks


title: Looking for appropriate time series econometrics materials, with Stata examples

text: I want to know whether my data are trending, co-integrated, have unit roots, have autocorrelation, or have heteroskedasticity. I want to do all this using Stata, so please suggest material or direct me a book on these topics


title: Probability function of three Random Variables multiplied, solidifying intuition

text: I ran across this exercise: > Let $T$ be a random variable distributed as a $\text{Bernoulli}(p)$, $U$ be > a random variable distributed as a $\text{Bernoulli}(q)$ and $W$ be a random > variable distributed as a $\text{Poisson}(\lambda)$. Find the probability > function of $X =T \cdot U\cdot W\,$ . I first calculated: $ \ \ P(T\cdot U = 1) = P(T=1)\cdot P(U=1) = pq \\\ P(T \cdot U = 0) = 1 - P(T\cdot U=1) = 1 - pq $ Then I made a mistake thinking that: $ P(X=0) =P(T\cdot U = 0 , W> 0) +P(T \cdot U = 1, W= 0) $ was the correct formula for the probability of $ X = 0$, this seemed reasonable and symmetrical to me, my intuition lead me astray. I then thought that the $ P(X = 0) = P(T\cdot U= 0) + P(W = 0)$ For the second time, my intuition lead me astray. The solutions highlight as the correct solution : $P(X= 0) = P(T\cdot U = 0) + P(T\cdot U = 1, W = 0) $ Intuitively I explained this result to myself as: if $ P(T\cdot U = 0) $ then I do not need to check the other random variable because anything times 0 is 0. If $ P(T\cdot U = 1) $ then I need the Poisson to be 0 to have $ P(X = 0) $. But because my intuition failed me I would like to understand which probability law the solution is using to solve this exercise, It seems to me not the Total Probability law or any form of Bayes rule. A Formal proof of the equality being used would be helpful to me. Also If you want to point out in witch way my first two attempts where wrong I would be grateful.


title: Usage of HMM in quantitative finance. Examples of HMM that works to detect trend / turning points?

text: I am discovering the marvellous world of such called "Hidden Markov Models", also called "regime switching models". I would like to adapt a HMM in R to detect trends and turning points. I would like to build the model as generic as possible so that I can test it on many prices. **Can anyone recommend a paper?** I have seen (and read) (more than) a few but I am looking for a simple model that is easy to implement. Also, **what R packages are recommended?** I can see there is a lot of them doing HMM. I have bought the book "Hidden Markov models for time series: an introduction using R", let see what's in it ;) Fred


title: Finding the average GPS point

text: I need to write a program to find the average GPS point from a population of points. **In practice the following happens:**   * Each month a person records a GPS point of the same static asset.   * Because of the nature of GPS, these points differ slightly each month.   * Sometimes the person makes a mistake a records the wrong assest at a completely different location.   * Each GPS point has a certainty weight (HDOP) that indicates how accurate the current GPS data is. GPS points with better HDOP values are preferred over lower ones.. **How do I determine the following:**   * Deal with data with 2 values vs. a single value like age. (Find the average age in a population of people)   * Determine the outliers. In the example below these would be [-28.252, 25.018] and [-28.632, 25.219]   * After excluding the outliers, find the average GPS point in this it might be [-28.389, 25.245].   * It would be a bonus if can work the "weight" provided by HDOP value for each point. ![alt text](http://i.stack.imgur.com/ecYHJ.png)


title: Is there a problem with the random effects structure in this lme model

text: I have a situation were I've conducted 18 sampling sessions. Each session is done within one of two streams. In each stream I take ten water samples at evenly spaced distances downstream of a DNA source. These samples are taken at different flows and there are different biomasses in the DNA source for each session. So, there are ten repeated measurements within a session. I've set it up so that the ten distance measurements are allowed to have random slopes and intercepts within each session. The continuous Flow and Biomass variables are treated as fixed as well as the categorical stream variable. The model looks like this:               Mod1 <- lme(DNACount ~ Distance + Flow + Stream + Biomass +                   Distance*Flow + Distance*Stream + Flow*Stream +                    Distance*Biomass + Flow*Biomass,                 random = ~Distance | SampleSession,                 data = DNASig)      If I plot predicted values at different flows within the two streams (using new data), one stream makes sense. I see different slopes and intercepts for the plotted lines. However, the other stream does not seem to have a random slope. It is just plotting parallel lines. Any thoughts on why this is?


title: Conditional probability distribution estimation from numerical data

text: I am trying to come up with a conditional probability distribution $P(B|A)$, over the continuous variables which look something like this:               A    B         ----------     2    5         5    7         8    10        10   15     12   25      ..........      etc. How do I go about finding what distribution will fit my case? From what I understand, it is a 2-step process. First is that I try to find what would be a good distribution to use to model this and then estimate the parameters assuming that distribution which can then be used in my calculations.


title: How do I state an ARIMA(0,1,4)x(0,1,1)12 in terms of Yt etc and not in terms of the backshift operator

text: I would like to state an ARIMA(0,1,4)x(0,1,1) model (written in format (p,d,q)x(P,D,Q) model without using the backshift operator. Thanks in advance


title: Interpretation of simple predictions to odds ratios in logistic regression

text: I'm somewhat new to using logistic regression, and a bit confused by a discrepancy between my interpretations of the following values which I thought would be the same:   * exponentiated beta values    * predicted probability of the outcome using beta values.  Here is a simplified version of the model I am using, where undernutrition and insurance are both binary, and wealth is continuous:               Under.Nutrition ~ insurance + wealth      My (actual) model returns an exponentiated beta value of .8 for insurance, which I would interpret as: > "The probability of being undernourished for an insured individual is .8 > times the probability of being undernourished for an uninsured individual." However, when I calculate the difference in probabilities for individuals by putting in values of 0 and 1 into the insurance variable and the mean value for wealth, the difference in undernutrition is only .04. That is calculated as follows:               Probability Undernourished = exp(β0 + β1*Insurance + β2*Wealth) /                                  (1+exp(β0 + β1*Insurance + β2*wealth))      I would really appreciate it if someone could explain why these values are different, and what a better interpretation (particularly for the second value) might be. * * * **Further Clarification Edits**   As I understand it, the probability of being under-nourished for an uninsured person (where B1 corresponds to insurance) is:               Prob(Unins) = exp(β0 + β1*0 + β2*Wealth) /                   (1+exp(β0 + β1*0+ β2*wealth))      While the Probability of being under-nourished for an insured person is:               Prob(Ins)= exp(β0 + β1*1 + β2*Wealth) /                (1+exp(β0 + β1*1+ β2*wealth))      The odds of being undernourished for an uninsured person compared to an insured person is:               exp(B1)      Is there a way to translate between these values (mathematically)? I'm still a bit confused by this equation (where I should probably be a different value on the RHS):               Prob(Ins) - Prob(Unins) != exp(B)      In layman's terms, the question is why doesn't insuring an individual change their probability of being under-nourished as much as the odds ratio indicates it does? In my data, Prob(Ins) - Prob(Unins) = .04, where the exponentiated beta value is .8 (so why is the difference not .2?)


title: For Logistic Regression, how do you interpret predicted probabilities?

text: Hi I am unsure how to interpret the predicted probabilities for my logistic regression model. I am trying to work out if these variables effect the chances of becoming obese. How would I interpret these results? REGION   Reference- East -North 0.62 -South 0.54 -West 0.47 The link suggested that makes my question a duplicate is insightful however it does not explain what they mean compared to reference category. For my odds ratio is know that the West were (odds=0.5) 'less' likely than east to be obese. what can I infer from the predicted probability? or is it not like that? do you just say the probability of having obesity and living in the west is 0.47? Many thanks :-)


title: Running a Probit on Survival-Time Data?

text: Can I run a probit on survival time data? It's discrete-round, and I want to look at whether lagged variables affect the failure event. I am, however, getting negative coefficients for a probit regression, and I'm not sure that should be happening. Help?


title: How to find P-value of Mutual Information Correlation (MIC) in R?

text: With reference to the article. I have a query regarding p-value calculation. For example, I have MIC strength value 0.1643 with sample size 30 and the rest parameters are default. As authors given p-values here How can we compute the p-value for the above MIC strength value 0.1643?


title: What exactly is the equation for SVM classification for new example?

text: I understand that in the case of Logistic Regression, we simply multiply our weights with Input example for classification. But what exactly is the equation that we calculate in the case of SVM to predict? I have already gone through Matlab's svmclassify, and it doesn't help. I will be absolutely grateful for a response.


title: multiple imputation with binary variables

text: I have 54 missing values in my dataset of 459 cases. Variables are all binary (0-1). I want to try a multiple imputation to avoid a listwise deletion, using the `mi` package for R. However, I ran a first analysis and the message was `mi did not converge`. I simply used the command:               mi(mydata)      leaving the default 30 iterations and 3 imputations. **Question: what are the best practices to do a multiple imputation for binary variables using the`mi` package?**


title: Transforming Rank-Ordered Variable

text: I am currently working on an archival dataset (General Social Survey 1985). In a measure, participants are asked to rank 5 job characteristics by degree of preference. In other words, participants decide what characteristic is the most important, what is the second most important, and so on, in a hypothetical job. I am particularly interested in how they ranked 2 of the 5 characteristics, because they match very well with a construct I am studying. These 2 characteristics actually tap into both extremes of my construct. As an illustration, if my construct of interest was job autonomy, one of the characteristics would be "Making important decisions in my job" while the other one would be "being closely supervised." In my model, the construct "autonomy" is the predictor of another construct in this dataset. The other construct is continuous. If the data for characteristic preference was simple likert scales, things would be MUCH easier. I would probably try to combine the two items into a single score. The problem is I do not know how to re-code / analyze this type of ranking for my two job characteristics. What I have done so far is calculate a difference score between "Making important decisions in my job" and "being closely supervised", which enables me to capture a relative preference for one or the other. What would you suggest me to do?


title: A binomial random number generating algorithm that works when n*p is very small

text: I need to generate binomial random numbers: > For example, consider binomial random numbers. A binomial random number is > the number of heads in N tosses of a coin with probability p of a heads on > any single toss. If you generate N uniform random numbers on the interval > (0,1) and count the number less than p, then the count is a binomial random > number with parameters N and p. In my case, my N could range from 1*10^3 to 1*10^10. My p is around 1*10^(-7). Often my n*p is around 1*10^(-3). There is a trivial implementation to generate such binomial random number through loops:               public static int getBinomial(int n, double p) {       int x = 0;       for(int i = 0; i < n; i++) {         if(Math.random() < p)           x++;       }       return x;     }      This native implementation is very slow. I tried the Acceptance Rejection/Inversion method [1] implemented in the Colt (http://acs.lbl.gov/software/colt/) lib. It is very fast, but the distribution of its generated number only agrees with the native implementation when n*p is not very small. In my case when n*p = 1*10^(-3), the native implementation can still generate the number 1 after many runs, but the Acceptance Rejection/Inversion method can never generate the number 1 (always returns 0). Does anyone know what is the problem here? Or can you suggest a better binomial random number generating algorithm that can solve my case. [1] V. Kachitvichyanukul, B.W. Schmeiser (1988): Binomial random variate generation, Communications of the ACM 31, 216-222.


title: Random intercept for different groups of subjects

text: Is there a way to specify two different random intercept terms for two sub- groups of subjects within a single mixed model? The reason why I'm interested is that it seems that the inter-personal variability of one group of subjects on my task is greater than the inter-personal variability of the other. I'm using lmer and currently specifying the subject factor as a random intercept, which pools all my participants into a single estimate, but I would like to have two separate estimates for the two sub-groups. Thank you!


title: Multinomial logistic regression vs binary logistic regression

text: Lets say we have a dependent variable $Y$ with few categories and set of independent variables. What are the advantages of Multinomial logistic regression over set of binary logistic Regressions? By set of binary logistic regression I mean that for each category $y_{i} \in Y$ we build separate binary logistic regression model with target=1 when $Y=y_{i}$ and 0 otherwise.


title: ROC curves and AUC in simulations to compare models

text: I am using ROC curves to compare different methods but not sure if I need to re-simulate datasets using different seeds in R in order to reduce the "by- chance" issue for a particular output. Here is a brief outline of my simulation:   1. The function `generate.data` is used to simulate data of some distribution, and by simulation, I know which data are true positives. **The random number generator is controlled by fixing the`seed` in R**.   2. The function `check.models` is used to test a total of 5 methods, and return the quantities used to draw a ROC curve for each method. Also for each curve (method), the AUC is reported.   3. The function `plot.roc` is used for plotting. In step #1, there are some other factors to change so that the data are under different "alternatives". When I run steps #1 and #2 above using `seed=123` and pick up the method with the highest AUC, I got one set of results. However, when I re-run using a different seed (say `seed=456`), I got another set of results **not identical to the first run**. Therefore, I think rigorously I should run my simulation across different `seed`'s in R to generate data in step #1, so that the "by-chance" issue of using a particular dataset is reduced. Am I correct? If so, then I should report the average of the AUC's for each method across (say, 1000) simulations, and pick up the highest among the methods compared? Thanks!


title: Statistical test on percentage data with replicates

text: Imagine I am doing an experiment where I want to look at the cellular composition of tumour samples I have 3 samples with phenotype A (e.g. metastasizing) and 3 samples with phenotype B (e.g. not metastasizing) The samples with phenotype A consist of (respectively) 30%, 32% and 33% "green cells", 28%, 27% and 28% "blue cells" and 42%, 41% and 39% "red cells". The samples with phenotype B consist of (respectively) 10, 9 and 11 "green cells", 50%, 49% and 47% "blue cells" and 40%, 42% and 42% "red cells". How can I tell if there is a significantly different % of the different cell types between the different types tumours (e.g. between phenotypes)? Normally it would require 3 t-tests. However because they are %s, I suppose that the assumptions of a t-test are not met. Is there an analogous test for %s?