title: Specifying a correct full model for an unbalanced repeated measures design in R lme4

text: Recently, I have done a fairly complex experiment, and I am having trouble coming up with a model that is suitable for the data. I have spent a few days reading about, e.g., when random effects should be nested or crossed, and which variables should be included in a full model. Yet, the literature that is readable for a non-statistician like myself is usually limited to two nesting levels, whereas I (may) have more. The literature that does seem to apply, is too complex for me to understand. I hope you can help me. My aim is to specify a full model, and start model simplification from there. The experiment I have done, was as follows:   * All participants completed three tasks, in a fixed order (Task1, Task2, Task3).   * Task 1 consisted of 30 trials, Tasks 2 and 3 both consisted of 10 trials.   * Each task was completed three times by each participant (Round1, Round2, Round3).   * The manipulation consisted of one factor with three levels (Condition1, Condition2, and Condition3).   * The conditions were tied to the rounds, so each participant completed Task 1 in each condition (e.g., 30 trials in Condition1, followed by 30 trials in Condition2, followed by 30 trials in Condition3), before moving on to Task 2.   * I used six different stimuli to manipulate these conditions (two for each condition; Stimulus1.1, Stimulus1.2, Stimulus2.1 .. Stimulus3.2)   * The dependent variable is binary To complicate matters even more, we used 6 different orders of presenting the stimuli. I think it's easiest to demonstrate the data structure using an example:               order <- rbind(     c('Stimulus1.1', 'Stimulus2.1', 'Stimulus3.1', 'Stimulus2.2', 'Stimulus3.2', 'Stimulus1.2', 'Stimulus3.1', 'Stimulus1.2', 'Stimulus2.1'),     c('Stimulus1.1', 'Stimulus3.1', 'Stimulus2.1', 'Stimulus3.2', 'Stimulus2.2', 'Stimulus1.2', 'Stimulus2.1', 'Stimulus1.2', 'Stimulus3.1'),     c('Stimulus2.1', 'Stimulus3.1', 'Stimulus1.1', 'Stimulus3.2', 'Stimulus1.2', 'Stimulus2.2', 'Stimulus1.2', 'Stimulus2.1', 'Stimulus3.1'),     c('Stimulus2.1', 'Stimulus1.1', 'Stimulus3.1', 'Stimulus1.2', 'Stimulus3.2', 'Stimulus2.2', 'Stimulus3.2', 'Stimulus2.1', 'Stimulus1.1'),     c('Stimulus3.1', 'Stimulus2.1', 'Stimulus1.1', 'Stimulus2.2', 'Stimulus1.2', 'Stimulus3.2', 'Stimulus1.1', 'Stimulus3.2', 'Stimulus2.1'),     c('Stimulus3.1', 'Stimulus1.1', 'Stimulus2.1', 'Stimulus1.2', 'Stimulus2.2', 'Stimulus3.2', 'Stimulus2.1', 'Stimulus3.2', 'Stimulus1.1'))          test <- expand.grid(trial=1:30, task=c('Task1', 'Task2', 'Task3'), round=c('Round1', 'Round2', 'Round3'), pp=1:12)     test <- test[! (with(test, trial > 10 & task == 'Task2')),]     test <- test[! (with(test, trial > 10 & task == 'Task3')),]     test$taskround <- factor(paste(test$task, test$round, sep=':'))     test$task <- factor(test$task)     test$round <- factor(test$round)     test$stimulus <- factor(unlist(lapply(1:nrow(test), function(x) {order[1 + (test[x, 'pp'] %% 6), as.numeric(test[x,'taskround'])]})))     test$condition <- factor(paste('Condition', substr(as.character(test$stimulus), 9,9), sep=''))     test$response <- factor(rbinom(nrow(test),1, prob=.95))     test$pp <- factor(test$pp)      The resulting data structure is:                  trial  task  round pp    taskround    stimulus  condition response     1      1 Task1 Round1  1 Task1:Round1 Stimulus1.1 Condition1        1     2      2 Task1 Round1  1 Task1:Round1 Stimulus1.1 Condition1        1     3      3 Task1 Round1  1 Task1:Round1 Stimulus1.1 Condition1        1     4      4 Task1 Round1  1 Task1:Round1 Stimulus1.1 Condition1        1     5      5 Task1 Round1  1 Task1:Round1 Stimulus1.1 Condition1        1     6      6 Task1 Round1  1 Task1:Round1 Stimulus1.1 Condition1        1     7      7 Task1 Round1  1 Task1:Round1 Stimulus1.1 Condition1        1     8      8 Task1 Round1  1 Task1:Round1 Stimulus1.1 Condition1        1     9      9 Task1 Round1  1 Task1:Round1 Stimulus1.1 Condition1        1     10    10 Task1 Round1  1 Task1:Round1 Stimulus1.1 Condition1        1      I'm mainly interested in the main effect of Condition. I am anticipating that the response tendency may differ between tasks and participants. This would translate into a different intercept for each task and participant, correct? I am also anticipating that the manipulation may affect some participants more than others; and I suspect that different participants may approach different tasks in a different way (e.g., Task 1 may elicit more successes that Task 2, but participant 1 may be more sensitive to this task aspect than participant 2). My current model is: `glmer(response ~ condition + (1 | pp/task/round) + (0 + task + condition | pp), data=test, family=binomial)` However, I am not sure if this model is correct. My questions are:   1. I am uncertain about the nesting of the random effects in the first clause. I believe `task/round` is correct, but I am not sure whether I should think of these as being nested under participants, or crossed with participants as follows: `glmer(response ~ condition + (1 | task/round) + (1 | pp) + (0 + task + condition | pp), data=test, family=binomial)`.   2. I am not sure at all about the `(0 + task + condition | pp)` random effects. I get correlations between the task and condition random effects, and I am not sure whether I would want/need these.   3. I removed the intercept here, because it was already included in the first random effects clause, but I'm not sure about that either.   4. Whether this model takes into account any order effects, or whether I should explicitly model those?   5. Whether I should take into account differences between stimuli (and/or the possibility that different participants react differently to different stimuli), and how I should model these (crossed, or nested)?   6. I'm worried about fatigue effect (lineair or quadratic) over the course of the experiment AND within tasks. I would like to include these effects into the model as well, but again, I'm not sure how to do that. So, in summary, I have lots of questions about analyzing these data. I think I am mainly interested in recommendations on how to approach this dataset, and suggestions for a suitable full model. Any other tips/suggestions are very welcome as well.


title: Rejection regions nested or not?

text: When varying the significance level, the rejection regions can be chosen to be nested or not nested. I was wondering what some theoretical and practical considerations are in using either nested or non-nested rejection regions? Thanks and regards!


title: Data precision independent of the data range

text: I am working with the execution times of an algorithm and I need to evaluate their variability. Standard deviation is a metric that tells you how much a given distribution variates from its mean. According to "Definition of "precision" of an estimator" question, in stpasha's answer he says that > Precision can be defined as 1/variance, or 1/standard deviation But standard deviation depends on the magnitude of data. So, what I am looking for is the process used to extract how much the data variates, independently from its magnitude. For instance, a distribution like {100,101,102,99,101} has bigger precision than {0.1,0.2,0.3,0.5} but its standard deviation is way bigger. Should I just divide by the distribution's average? This might be just a very trivial case of statistics (my fault for not remembering those classes), but I would like to know what is the standard procedure to evaluate this and how is it called.


**Background:** I would like to offer readers a statistic that conveys the similarity of two sets of numbers. I thought that I had what I needed with correlation coefficient (indeed, I have a coefficient of determination/r2 of 98.7). But I've just realized that r and r2 conveys only that the two sets of numbers are correlated in their upward and downward movement; and what I want (also) to convey is that they are also very nearly identical. But yet I feel like regression analysis is more than I need (or, more precisely, more than I can do with Excel). To be a little more concrete, imagine a data set with monthly earnings for two stores in a strip mall. Correlation coefficients would capture whether the stores' earnings were covariate. But I want to demonstrate also that the monthly earnings also happen to be very nearly the same actual number. On a graph, they are visually one is on top of the other. But is there a statistical metric that captures that? **Question:** What is a measure of association that conveys actual similarity of values between two variables and not just correlation?


title: Probability of throwing same sequence of 13 heads or tails

text: Simple (I hope) probability question I'd love it if someone could answer for me. You flip a coin 13 times (two outcomes: heads or tails). Before the coin tossing, I've written down the 'winning sequence', eg. a random sequence of 13 outcomes, heads or tails for each toss. For example: Heads, Heads, Tails, Heads, Tails, Tails, Tails, Heads, Tails, Heads, Heads, Tails, Heads. What are the chances of someone tossing the coin and getting the same sequence (heads or tails in the exact same sequence) as me?


title: Can I insert an observation (evidence) to a Winbugs model?

text: Greetings, Is it possible to use evidence in a Winbug model? For example, a random variable in a model has been observed, and I'd like to update the other variables in the model, pretty much the same update perfomed in tools like Smile, or other inference software. Gibbs sampling is supposed to use observed values in full conditional distribution when there observations/evidence are entered into the model, but I am not sure if Winbugs allows this. Regards edit to clarify: The winbug documentation says a stochastic node is a data node if it has been observed, but a stochastic node is described via a distribution as in some_var ~ dbin(theta,n) if it has been observed, then I'd like to tell this to winbugs without losing the semantics of the stochastic node, expressing something like "some_var has this distribution, and it has been observed to have this particular value". so how do I do that? by declaring some_var as I've done above and than setting a value to it as in some_var = 5 ? Would that express what I want to express? In this case, in each observation of a node in a bayesian network, I'd need to redefine the winbugs model, (quite likely) replacing initial values of the unobserved nodes with the outcomes of the previous simulation. In short, I'm trying to understand how to perform updates on a Bayesian network similar to message passing in exact inference, but using Gibbs sampling instead, via Winbugs.


title: How to analyse a study where measures are taken repeatedly from right and left legs in five walking conditions?

text: There are 13 subjects in this study, each subject was "repeatedly measured" right and left legs on five walking conditions. The variables of this data set is: ID, Y, Leg, Conditions. The research questions are to find the difference between legs and conditions. At beginning, I plan to use RM-ANOVA: $Y = \text{Leg} + \text{Conditions} + \text{Leg}\times\text{Conditions}$. But the RM-ANOVA might not be appropriate because right and left legs are also correlated (from same subject). Does anyone have an idea how to analyze this type of data?


title: Difference between full house and two pair in consideration of numbers chosen?

text: So I am a bit confused about why the number of ways to get a full house is: $ 13 \choose 1 $$ 12 \choose 1$$ 4 \choose 3$$ 4 \choose 2$ as opposed to: $ 13 \choose 2 $$ 4 \choose 3$$ 4 \choose 2$. Alternatively, with 2 pairs, why it is: $ 13 \choose 2 $$ 4 \choose 2$$ 4 \choose 2$$4 \choose 1$$11 \choose 1$ over: $ 13 \choose 1 $$ 12 \choose 1 $$ 4 \choose 2$$ 4 \choose 2$$4 \choose 1$$11 \choose 1$


title: Jensen-Shannon divergence for bivariate normal distributions

text: Given two bivariate normal distributions $P \equiv \mathcal{N}(\mu_p, \Sigma_p)$ and $Q \equiv \mathcal{N}(\mu_q, \Sigma_q)$, I am trying to calculate the Jensen-Shannon divergence between them, defined (for the discrete case) as: $JSD(P\|Q) = \frac{1}{2} (KLD(P\|M)+ KLD(Q\|M))$ where $KLD$ is the Kullback-Leibler divergence, and $M=\frac{1}{2}(P+Q)$   I've found the way to calculate $KLD$ in terms of the distributions' parameters, and thus $JSD$. My doubts are:   1. To calculate $M$, I just did $M \equiv \mathcal{N}(\frac{1}{2}(\mu_p + \mu_q), \frac{1}{2}(\Sigma_p + \Sigma_q))$. Is this right?    2. I've read in [1] that the $JSD$ is bounded, but that doesn't appear to be true when I calculate it as described above for normal distributions. Does it mean I am calculating it wrong, violating an assumption, or something else I don't understand?


title: ANCOVA with multiple instances of the between-subject covariate

text: If I find that my covariate (reaction time) alters over the length of my experiment (e.g. due to fatigue), can I somehow build that into my model?   So what I am saying is that the effect of my covariate is not constant (between subjects and within subjects).


title: What's the difference between "measurement method" and "measurement procedure"?

text: The ISO VIM defines them as: > **measurement method:** generic description of a logical organization of > operations used in a measurement. > > **measurement procedure:** detailed description of a measurement according > to one or more measurement principles and to a given measurement method, > based on a measurement model and including any calculation to obtain a > measurement result. But I'm still not sure what their difference is. Can anyone explain this to me, maybe with some examples?


title: How would life expectancy impact the calculation of disease prevalence?

text: I want to understand how to make calculations on the prevalence of a disease in a country population and the impact that the element of average life expectancy (of those suffering with the disease at time of diagnosis) has on this calculation. Are there any 'best practice' papers on making epidemiology calculations?


title: Principal component analysis -- backwards?

text: I have carried out a principal components analysis of six variables $A$, $B$, $C$, $D$, $E$ and $F$. If I understand correctly -- unrotated PC1 tells me what linear combination of these variables describes/explains the most variance in the data and PC2 tells me what linear combination of these variables describes the next most variance in the data and so on, I'm just curious -- is there any way of doing this backwards? Lets say I choose some linear combination of these variables -- e.g. $A+2B+5C$, could I work out how much variance in the data this describes?


title: How to interpret this residuals vs fitted plot for logistic regression using R

text: I am working on a logistic regression on some fundraising data where "gave" is a rare event (approx 3.5%). My current model has 64% accuracy on test data and an AUC of .604. When I run the standard plot(model) I am seeing some outliers on each plot that are segmented completely away from the rest of the data set. What is your first thoughts when you see segregation in the plot like this? Is this indicative of a problem or just the nature of the nature of such a rare event? ![Residuals vs fitted plot](http://i.stack.imgur.com/209Ri.png)


title: Should there be an "i" in a regression equation?

text: I have developed a general linear model to predict my dependent variable Y. I am unsure on how to present my equation. Should it be: (not sure if the i should be there or not) $Yi= .432 + .320 Age_i + .520 WE_i + .300 JP1_i + .210 JP2_i$ or $Y= .432 + .320 Age + .520 WE + .300 JP1 + .210 JP2$ Where: WE= work experience; JP= Job Position; JP1= floor level staff; JP2= Managers;


title: Evolution of slopes over time

text: Participants in a study are tested on their ability to recall items depending on the number of times these items are presented by day. We follow their performance over a time period of 4 days. So I have two independent variables: number of presentations and days. Basically, I want to know if there is a linear effect of the number of presentations and how this effect evolves over time. At first I was thinking about a linear model with one beta per dependent variable. But the issue is that doing that won't give account about an interaction, since if I understood correctly, the interaction term is entered as another variable formed by multiplying the values of the two other regressors. Finally, I'm thinking of using a model where I would fit the number of presentations regressor over each day independently and run an ANOVA on the betas to test the interaction effect. Would I be correct by doing this, or is there better approach to analysis?


title: Grouping data values so each group has the minimum range

text: Apologies in advance if this is very obvious or has been asked elsewhere already - I'm not very 'maths-y' so not sure how to ask this intelligently. What is the best way to group a set of data values into _n_ groups, such that each groups has the smallest possible range? Lets say my sample values are as follows:               {1, 2, 7, 8, 9, 15, 16, 16, 17, 18, 18, 20}      and I am wanting to create 3 groups with the least amount of spread in each group. We could do three groups like so:               {1, 2, 7}                         spread = 6     {8, 9}                            spread = 1     {15, 16, 16, 17, 18, 18, 20}      spread = 5     ---------------------------------------------                                        total = 12      Or, if we rework the first two groups we can reduce it:               {1, 2}                            spread = 1     {7, 8, 9}                         spread = 2     {15, 16, 16, 17, 18, 18, 20}      spread = 5     ---------------------------------------------                                        total = 8      So my question is: is there a mathematical way of determining how the groups should be formed? As a reverse version of the question, is there a way to specify not the **number** of groups, but the **maximum spread** a group can have, and figure out the groups such that they don't exceed the spread? This may have some crossover to programming algorithms so I will post it on Stack Overflow as well.


title: Nearest Neighbor Algorithm for Prediction

text: I am building a prediction model based on k-th-nearest-neighbor (KNN) method, as I have similar events. As the algorithm doesn't classify points to Sets, I was wondering if there is a way or criteria to set the number of neighbors. How to estimate the probability of predicted values based on the nearest neighbors.


title: What is the most likely meaning of a confidence interval?

text: I'm using a service (Optimizely, as it happens), which show a confidence interval but doesn't explain what it is. I'm assuming it is a Binomial Confidence Interval, but I have no idea what the % certainty of the range is, or the detailed meaning. Searching the Internet doesn't help - there is nobody with a clear explanation of what the default for this kind of thing is. For example, for 97 success out of 939 trials, Optimizely says "10.33% (±1.95)". What does that ±1.95 likely mean, _exactly_? And yes, I've already complained to them that they should document it. But clearly, there is some default that statisticians habitually use. What is it?


title: What is sparse regression model

text: I am learning the concepts of Sparse regression and facing initial hurdles in terminology. sparse regression model explains the definition of what is meant by sparse. When the number of samples $n$ is less than the signal dimension $p$ then we say it is sparse regression model.   1. For a model, $x_t = a1x_{t-1} + a2x_{t-2} + white gaussian noise$, the parameters $(a1,a2)$ do not vary with time and for $n=t$ samples we get only two parameters. Then, how come the paper says that $A \in R^{n \times p}$? What does this mean? $A$ is a vector of coefficients and not a matrix. Then what does $R^{n \times p}$ mean?   2. What I have understood is that there are more number of zero coefficients, but if that so then is there an upper bound which will indicate that the signal is sparse?   3. I am looking for a reference where I can find such a sparse AR and MA model. Can somebody please point out a link or help in creating such a model? Thank you


title: What is signal dimension

text: Estimating Unknown Sparsity in Compressed Sensing is a paper about sparse signal. I am just learning the concepts. In the first paragraph, it says that when the number of observation data samples $n$ is less than the signal dimension $p$, then the desired signal $x$ is sparse. I have never come across he term signal dimension. Can somebody please explain what is signal dimension ?


title: Inference test when data points are means

text: ### General background I've administered two types of questions to students multiple times. Type 1 I've administered six times, and type 2 I've administered five times. The two questions get at the same basic idea but in different ways. Each "question" is just a single multiple-choice question; students provide just one answer (i.e., a single bubble on a Scrantron sheet). I know the percentage correct on each administration of each question. That is, my data looks something like (35%, 32%, 24%, ...) for one of the questions and (20%, 24%, ...) for the other. I have a total of about 11 data points split into two categories. ### Question I want to know if students perform better on one type of question as compared to the other. Is there an inference test or class of tests that are appropriate here? ### Possibly important details   * A **different number of students** is involved with each administration, but they're all between 60 and 180.   * From other information I know that repeated administrations of the same question result in a gaussian distribution when looking at percentage correct.   * The standard deviations associated with the distributions just mentioned depend greatly on the particular question being asked (presumably due to differences in student backgrounds, etc.) Thus, I don't want to assume or prescribe some type of $\sqrt{Np(1-p)}$ or other standard deviation to each administration. ### What I've considered My initial inclination is to do some type of "weighted" t-test, but I couldn't find anything. Then I went down the route of looking at meta-analysis tests, but I'm quite unfamiliar with these. Are there inference tests that are meant to address this particular situation?


title: Confusion related to Kalman filters density view

text: I was reading this book related to Kalman filters and I didn't understand a couple of things. I have also attached the screenshot of the pages from the book where I had confusion. The book is Shumway, R.H. and Stoffer, D.S. Time Series Analysis and Its Applications: With R Examples From the first page, I didn't understand how this was derived ![enter image description here](http://i.stack.imgur.com/wEHbF.png) I mean how come we have $p(x_t|x_{t-1}) = f_w(x_t-\phi_{x_{t-1}})$ and has mean 0 and covariance Q. I mean it should be for the $w_t$ which has mean 0 and covariance Q not $p(x_t|x_{t-1})$ From the second page, I didn't understand the following ![enter image description here](http://i.stack.imgur.com/B0v7Q.png) I mean obviously $p(\mu_t|Y_{t-1} => N(\mu_t^{t-1},P_t^{t-1}))$. But I didn't get how we get that from the integral From the third page, I didn't get how this was derived ![enter image description here](http://i.stack.imgur.com/pKfGh.png) I know I might be asking a lot, but I have spent several days to exactly grasp the derivation of the equations in Kalman filters. I would really appreciate if someone could help me. Thanks in advance ![enter image description here](http://i.stack.imgur.com/Ru9wU.png) ![enter image description here](http://i.stack.imgur.com/whVb4.png) ![enter image description here](http://i.stack.imgur.com/A4diL.png)


title: What is a stationary Gaussian field?

text: I know what a Gaussian field is. However, I am not quite sure what is meant by stationary. I have seen this stationary thing at lots of places like stationary autoregressive processes etc but don't actually know what is meant by stationary.


title: How can I sample from a correlated multivariate Bernoulli distribution with known covariances?

text: I have $N$ Bernoulli random variables $X_1, ..., X_{N}$ with known parameters $p_1, ..., p_{N}$. They are dependent with known covariances. How can I sample from the joint distribution of the $X_i$? I would ideally like a more general approach than this.


title: Which test to use to see if tumor weight increases at different stages?

text: I have 60 cancer patients. I have checked their staging and the weight of their tumour in gram. One variable is cancer staging - Stage 1, Stage 2, Stage 3 and Stage 4 (in increasing order). Another variable is Tumour weight - a continuous variable. I want to check as stage increases from Stage 1 to 4, whether tumour weight also increases. Could you please suggest which test should I use?


title: What is the relationship between correlation coefficients and regression coefficients in multiple regression?

text: Consider the model $ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \epsilon$.   1. What is the relationship between the correlation coefficients $r_{y,x_1}$, $r_{y,x_2}$ and the regression coefficients $\beta_1$ and $\beta_2$?    2. In particular, how to interpret a situation where a particular correlation coefficient is statistically significant but the corresponding regression coefficient is not statistically significant?   3. I recall reading a concept called "partial correlation coefficients". Is that in some way relevant in the above context?


title: Pairwise Mahalanobis distance

text: I need to calculate the sample Mahalanobis distance in R between every pair of observations in a n x p matrix of covariates. I need a solution that is efficient, i.e. only n(n-1)/2 distances are calculated, and preferably implemented in C/RCpp/fortran etc. **I assume $\Sigma$, the covariance matrix, is unknown and use the sample covariance matrix in its place.** I am particularly interested in this question since there seems to be no "consensus" method for calculating _pairwise_ Mahalanobis distance in the `R` package, i.e. it is not implemented in the `dist` function nor in the `cluster::daisy` function. The `mahalanobis` function does not calculate pairwise distance without additional work from the programmer. This was already asked in the posting Pairwise Mahalanobis distance in R, but the topvoted solution seems incorrect (see printouts below):               library("MASS")     library("ICSNP")     set.seed(0)          x0<-mvrnorm(33,1:10,diag(c(seq(1,1/2,l=10)),10))     x1<-pair.diff(x0) #C-implementation.     dM<-mahalanobis(x1,colMeans(x1),var(x1))      Compare this to a correct but terribly inefficient (since n x n distances calculated) method:               dM2 = as.dist(apply(x0, 1, function(i) mahalanobis(x0, i, cov = cov(x0))))     plot(dM, dM2)      ![enter image description here](http://i.stack.imgur.com/jbnC7.png)               dM[1:10]     # [1] 12.046 5.835 9.545 8.152 11.294 7.658 13.847 9.062 9.404 6.632          dM2[1:10]     # [1] 23.49 11.79 17.55 14.57 24.10 15.02 25.62 19.84 18.01 14.65          (dM/dM2)[1:10]     # [1] 0.5128370 0.4949167 0.5438379 0.5596519 0.4686635 0.5098689 0.5405308     # [8] 0.4566752 0.5222861 0.4527394      The method from the original posting is off by a factor of 2, plus some strange unknown variation. (1) Am I missing something in the original post? The discrepancies between the two methods are concerning (I'd post there but I ~~don't~~ didn't have enough reputation) (2) This is easy enough to code myself in C, but I feel like something this basic should have a preexisting solution. Is there one? (3) Other solutions that fall short: HDMD::pairwise.mahalanobis() calculates n x n distances, when only n(n-1)/2 unique distances are required. compositions::MahalanobisDist() seems promising, but I don't want my function to come from a package that depends on rgl, which severely limits others' ability to run my code. Unless this implementation is perfect, I'd rather write my own. Anybody have experience with this function?


title: How I can I characterize the treasury yield curve as one number?

text: We know when the treasury yield curve inverts it is a sign that the economy is headed for trouble. I would like to analyze the relationship between the yield curve and the economy. My problem is coming up with a method to describe the shape of the yield curve using one number. I'd like to do something a little more complex than just measuring the slope between the near and far term instruments. Any suggestions?


title: What are the differences between "Mixed Effects Modelling" and "Latent Growth Modelling"?

text: I'm decently familiar with mixed effects models (MEM), but a colleague recently asked me how it compares to latent growth models (LGM). I did a bit of googling, and it seems that LGM is a variant of structural equation modelling that is applied to circumstances where repeated measures are obtained within each level of at least one random effect, thus making Time a fixed effect in the model. Otherwise, MEM and LGM seem pretty similar (eg. they both permit exploration of different covariance structures, etc). Am I correct that LGM is conceptually a special case of MEM, or are there differences between the two approaches with respect to their assumptions or capacity to evaluate different types of theories?


title: Information on how value of k in k-fold cross-validation affects resulting accuracies

text: I've been doing some Machine Learning, and have been using k-fold cross- validation to assess the generalisation performance of the algorithm. I've tried k-fold cross-validation with k = 5 and k = 200 and get very different results for Support Vector Machine classification.               k    SVM accuracy     -----------------     5    75%     200  94%      This seems like a huge difference in accuracy caused by changing the number of splits we're doing for the k-fold cross-validation. Is there any reason for this? I can't seem to find any references on studies that have been done investigating the effects of using different k values. Obviously, which k value I decide to use in my report gives completely different impressions of the quality of my classifier!


title: The minimum number of time points required on repeated measurement

text: @JeromyAnglim, thanks for your answer to "How to determine the sample size needed for repeated measurement ANOVA". My question: Is there a formula to estimate the minimum number of time points required for repeated measurements when we do a clinical trial? I know that there's a formula to estimate the sample size when we plan to do a clinical trial; I'm hoping for something like that. I am now planning to do a clinical trial to evaluate the effect of a new device compared to the old one. The patients will be randomized 2:1 to a study arm or controlled arm. My outcome variable is continuous. If the study is conducted over 1 year, I want to know what is the minimum number of time points I should get? Could you give me some reference?


title: Testing for collinearity and parallelism with the axes

text: Cross-posted from math; let me know if you need me to wait there for a bit before posting here I have a set of datapoints in $\mathbb{R}^3$ (integer or real values of $X$, $Y$ and $Z$, for, say. 30 points). I need to do a least squares regression for linearity. If the data set is linear, I need to see if it is close to vertical or horizontal. How could I do this?


title: How can I use logistic regression betas + raw data to get probabilities

text: I have a model fitted (from the literature). I also have the raw data for the predictive variables. What's the equation I should be using to get probabilities? Basically, how do I combine raw data and coefficients to get probabilities?


title: Fisher Test in R

text: Suppose we have the following data set:                               Men    Women         Dieting         10      30     Non-dieting     5       60      If I run the Fisher exact test in R then what does `alternative = greater` (or less) imply? For example:               mat = matrix(c(10,5,30,60), 2,2)     fisher.test(mat,alternative ="greater")      I get the p-value = 0.01588 and odds ratio = 3.943534. Also, when I flip the rows of the contingency table like this:               mat = matrix(c(5,10,60,30), 2,2)     fisher.test(mat,alternative ="greater")      then I get the p-value = 0.9967 and odds ratio = 0.2535796. But, when I run the two contingency table without the alternative argument (i.e., `fisher.test(mat)`) then I get the p-value = 0.02063.   1. Could you please explain the reason to me?    2. Also, what is the null hypothesis and alternative hypothesis in the above cases?   3. Can I run the fisher test on a contingency table like this:  `mat = matrix(c(5000,10000,69999,39999), 2,2)` PS: I am not a statistician. I am trying to learn statistics so your help (answers in simple English) would be highly appreciated.


title: Distribution of ratio $X/Y$ where $X$ is normal, $Y$ is half normal

text: I know the ratio of 2 normally distributed random variables is a cauchy distribution...which, of course, has no variance. But, if I have $X$ and $Y=|X|$ and take this ratio: $Z=X/Y$, is it possible to find the variance of the ratio of the normal distribution to the half normal? * * * **Edit** : (More details) I am missing something a little in the answers so far, but I am extremly grateful. Let me provide a little more background. Measuring data from sensors includes the process value and some noise. Xk = Vk + Nk, where Vk is the process measurement and Nk is the noise. There are two issues, (1) Estimating Vk and (2) detecting a shift or change point. I can estimate Vk using an exponential filter. so F0k = a*(Xk) + (1-a)*F0k-1. The expected value of F0 -> E(F0) = E(X). This means exponential filtering will give me my expected value. The variance of F0 -> var(F0) = (a/(2-a))*var(X) At this point I am interested in residuals: RES = Xk - F0k. This quantity should be 0 in the steady state. That is to say, the expected value is 0. I also assume that the value will be normally distributed (don't yell at me about this assumption). I can estimate the expected value online using a simple exponential filter (basically the same as an average): F1k = b*(RES) + (1-b)*F1k-1. The expected value of F1 -> E(F1) = E(RES) = 0. The variance of F1 -> var(F1) = (2/(2-a))(b/(2-b))*var(X) Then I can take abs(RES) and use exponential filtering again: F2k = b*(abs(RES)) + (1-b)*F2k-1 The expected value of F2 -> E(F2) = E(abs(RES)) = [sqrt(2/(2-a))*sqrt(2)/(sqrt(pi))] _std(X). The variance of F2 -> var(F2) = (2/(2-a))_ (b/(2-b))*(1-(2/pi))*var(X) All of these expected values and variances are just from either the exponential filtering properies or the expected value and variance of the normal or half-normal distribution. My question was, what is the variance (and expected value) of Z = F1/F2 F1 will be some normally distributed value with expected value = 0. F2 will be some half-normally distributed value with expected value based on the standard deviation of x. Can you find the var of Z? I understand that X/abs(X) would just give you a bernoulli, but it seems it doesn't give me exactly what i'm looking for? * * * **Edit** : (More details) Here are some interesting simulation results: F1 - exponential filter of Residual (noise added to original signal was Gaussian - result is normally distributed) ![enter image description here](http://i.stack.imgur.com/3cKs2.png) F2 - exponential filter of abs(Residual) (half-normal distribution) ![enter image description here](http://i.stack.imgur.com/pEgGZ.png) Z - F1/F2 - simulated distribution of ratio of F1/F2 (notice the large variance around expected value of 0) ![enter image description here](http://i.stack.imgur.com/AfEP7.png)


title: Graph interpretation hetroskedasticity

text: ![enter image description here](http://i.stack.imgur.com/SeF6U.jpg) I ran a command of rvf plot . Does this graph show no heteroskedasticity?


title: How to use cross-section country samples of different years in multiple regression analysis?

text: I want to run a multiple regression analysis measuring the effect of various independent variables on a continuous dependent variable measuring the strength of a political institution. The problem is that the data is measured in different years. Thus I might have data for countries A and B only for the year 2009 and for countries C and D only for the year 2010. In order to increase my sample, I would like to use all countries. How can I do that, ie which statistical method is most adequate? Thank you sooo much!


title: Obtaining and interpreting bootstrapped confidence intervals from hierarchical data

text: I'm interested in obtaining a bootstrapped confidence interval on quantity X, when this quantity is measured 10 times in each of 10 individuals. One approach is to obtain the mean per individual, then bootstrap the means (eg. resample the means with replacement). Another approach is to do the following on each iteration of the bootstrapping procedure: within each individual, resample that individual's 10 observations with replacement, then compute a new mean for that individual, and finally compute a new group mean. In this approach, each individual observed in the original data set always contribute to the group mean on each iteration of the bootstrap procedure. Finally, a third approach is to combine the above two approaches: resample individuals then resample within those individuals. This approach differs from the preceding approach in that it permits the same individual to contribute multiply to the group mean on each iteration, though because each contribution is generated via an independent resampling procedure, these contributions may be expected to vary slightly from eachother. In practice, I find that these approaches yield different estimates for the confidence interval (ex. with one data set, I find that the third approach yields much larger confidence intervals than the first two approaches), so I'm curious what each might be interpreted to represent.


title: Libraries for forest and funnel plots

text: Can anyone recommend me an open source graphic library to create forest and funnel plots? I was aiming at using it on a Java desktop application.


title: Good Literature on effect sizes

text: reporting effect sizes is pretty common nowadays, but I couldn't yet find any good book that compares different effect sizes values (d, g, r, etc.). I know the article by Rosnow and Rosenthal (2009), but they have no examples of computation. For me sometimes it's not clear, which standard deviation or standard error is used and what range the effect size can take. To be more specific. Right now I want to report results of independent as wel as dependent sample t-tests. I know r for effect sizes, but I would like to compute Cohen's d for example. How do I do it? what is the range of it? Thanks in advance!


title: LIBSVM cross-validation - number of folds vs. sample size

text: I am using LIBSVM for classification and the default RBF Kernel. To get optimal parameters for C and gamma, I am using the grid.py tool included in the LIBSVM package. It is based on a k-fold cross-validation, as far as I know. I have a training and a testing set of 15 451 instances with 5 features. Should I run the cross-validation on the whole training set or split it into e.g. a half and to the validation on it. And is there a rule for how many folds k I should chose in respect to the size of the cross-validation sample? For example, I wanted to do the cross-validation on the whole training set (15451 instances) with k = 50. Is this a good approach for getting the best parameter values?


title: birthday problem via simulation

text: I would like to plot the cumulative density function for this problem by running the following experiment $p$ times: generate random numbers between 1 and $n$ and stop when I get a repeat. The x axis will be the number of trials it took. How do I compute the y axis?


title: Significance level on real life

text: Recently I was reading a problem for one-tailed test, a fragment in the resolution demanded my attention: > "As noted in the preceding section, the level of significance, denoted by > **_α_** , is the probability of making a Type I error by rejecting **_Ho_** > when the null hypothesis is true as an equality. The decision maker must > specify the level of significance. If **the cost of making** a Type I error > is high, a small value should be chosen for the level of significance. If > the cost is not high, a larger value is more appropriate". _("Statistics for > Business and Economics", Anderson, 11ed, page 357)_ We know the concept for level of significance, but the question is... based on what we assign a value for significance level? when a high value? when a small value? how do we determine the cost of making a certain type of error into business or research? In few words, how does the level of significance is worked on real life?


title: Fisher's Exact Test and Hypergeometric Distribution

text: I wanted to understand fisher exact test better, so I devised up the following toy example, where f and m corresponds to male and female, and n and y corresponds to "soda consumption" like this:               > soda_gender              f m       n 0 5       y 5 0      Obviously, this is a drastic simplification, but I didn't want the context to get in the way. Here I just assumed that males don't drink soda and females drink soda, and wanted to see if the statistical procedures come to the same conclusion. When I run the fisher exact test in R, I get the following results:               > fisher.test(soda_gender)     Fisher's Exact Test for Count Data          data:  soda_gender     p-value = 0.007937     alternative hypothesis: true odds ratio is not equal to 1     95 percent confidence interval:      0.0000000 0.4353226     sample estimates:     odds ratio               0       Here, since p-value is 0.007937, we would conclude that gender and soda consumption are associated. I know that fisher-exact test is related to hypergeomteric distribution. So I wanted to get the similar results using that. In other words, you can view this problems as following : there are 10 balls, where 5 are labeled as "male", and 5 are labeled as "female", and you draw 5 balls randomly without replacement, and you see 0 male balls. What is the chance of this observation? To answer this question, I used the following command:               > phyper(q=0,m=5,n=5,k=5,lower.tail=TRUE)     [1] 0.003968254      My questions are: 1) How come the two results are different? 2) Is there anything incorrect or not rigorous in my reasoning above?


title: Using mvmeta to perform a network meta-analysis

text: A few months ago, I posted these two questions about practical considerations in collecting data for a network meta-analysis and the available R packages that can analyze such data. For those unfamiliar with network meta-analysis, it is a method that compares a set of treatments against each other. Part of what's attractive about this method is that multivariate extensions of network meta-analysis can accommodate studies that compare more than two treatments against each other; traditional meta-analytic methods can only accommodate two treatments. I am currently attempting to use the mvmeta package to conduct the meta- analysis. However, I'm a little stuck in figuring out how to implement the meta-analysis. As I understand it, network meta-analysis can be implemented at either the arm level or the contrast level. In my own meta-analysis, however, authors rarely report the outcome at the arm level, and much of my data had to be extracted using test statistics from the comparison between arms (e.g., a t-test comparing treatment A to treatment B). Therefore, I'm using the contrast-based approach, which accommodates studies with more than three arms by modeling the different comparisons between the arms as separate outcomes. For a study with k arms, k - 1 comparisons fully represent the differences between the arms. So, one must choose a reference treatment and calculate the k - 1 comparisons to that reference treatment. For example, for a study comparing treatments A, B, and C to each other, one would calculate separate effect sizes for the AB and AC comparisons. What I don't understand is how to extend this framework to represent studies that don't include the reference treatment. So, in the above example, I don't understand how to accommodate a study that compares treatments B and C. One possibility is that I would simply run a new analysis with B as the reference treatment. Do any of you know if this is actually the case? For reference, I have included some R code that I wrote that I believe performs a network meta-analysis using mvmeta below. My question, restated in the context of the sample R code, is that I don't see how to represent studies that observe BC comparisons.               library(mvmeta)          # Construct data     AB <- c(.5, .3, .2, .2, NA, NA, .1)     AC <- c(NA, .3, .4, NA, .2, .1, .5) # Where are the BC comparisons?     d <- data.frame(AB, AC)          n1 <- 50     n2 <- 50     nT <- 150 # Total N for 3-arm studies in this example          # Create list of N = 7 within-study variance-covariance matrices (S matrices).     # These are a function of the effect sizes and study Ns specified above.     S <- list(matrix(c(1/n1 + 1/n2 + .5/(2*(n1 + n2)), NA,                         NA, NA), ncol = 2, nrow = 2),                matrix(c(1/n1 + 1/n2 + .3/(2*nT), 1/n1 + .3 * .3 / (2*nT),                         1/n1 + .3 * .3 / (2*nT), 1/n1 + 1/n2 + .3/(2*nT)), ncol = 2, nrow = 2),               matrix(c(1/n1 + 1/n2 + .2/(2*nT), 1/n1 + .2 * .4 / (2*nT),                        1/n1 + .2 * .4 / (2*nT), 1/n1 + 1/n2 + .4/(2*nT)), ncol = 2, nrow = 2),               matrix(c(1/n1 + 1/n2 + .2/(2*(n1 + n2)), NA,                        NA, NA), ncol = 2, nrow = 2),               matrix(c(NA, NA,                        NA, 1/n1 + 1/n2 + .2/(2*(n1 + n2))), ncol = 2, nrow = 2),               matrix(c(NA, NA,                        NA, 1/n1 + 1/n2 + .1/(2*(n1 + n2))), ncol = 2, nrow = 2),               matrix(c(1/n1 + 1/n2 + .1/(2*nT), 1/n1 + .1 * .5 / (2*nT),                        1/n1 + .1 * .5 / (2*nT), 1/n1 + 1/n2 + .5/(2*nT)), ncol = 2, nrow = 2))          mod <- mvmeta(cbind(AB, AC) ~ 1, S = S, data = d, method = "reml")     summary(mod)


title: Modelling Issues in panel regression

text: I came across a paper that uses a panel data of US states (from 2000-2010) with the following model: $y_{it}=b_{1}x_{it} +b_{2}x_{it}*D_{it}+other vars + \alpha_{t} + \gamma_{i}$ where, $x$ and $y$ are continuous vars and $D$ is a dummy variable that takes a value of 1 if there is a decrease in $x$ from one year to another year and zero otherwise (i.e., for an increase or no change in $x$). The paper then interprets coefficient on $x_1$ as the effect of the increase in $x$ on $y$ and $b_{1}+b_{2}$ as the effect of a decrease in $x$.I was wondering whether this is a valid interpretation.


title: Hazard Function - Survival Analysis

text: I just started taking survival analysis class and I'm stumped on this question. Let $T_1,...,T_n$ random independent continuous variables, with hazard function of $h_1(t),...,h_n(t)$. $T=min(T_1,...,T_n)$. And we need to show that the hazard function of T is $\sum_j{h_j(t)} $ Any help or direction are welcome :)


title: Semi-discrete probability distribution

text: I am somewhat new to statistics, so bear with me. I have a dataset that contains integer values. I created a model for this dataset such that               y(x) = f(x) + c      where `f(x)` is my model and `c` is the error in my model. My `f` function is continuous, and therefore `f(x)` is continuous. I am assuming that `c` is noise, and I am trying to compute the distribution of `c`. The problem is that, while `c` is real-valued, it is not continuous, since it is defined as a discrete value minus a real-value (`c = y(x) - f(x)`, where `y(x)` is an integer). The histogram of `c` looks like an exponential distribution, except there are spaces in between (this site wouldn't let me upload my image for lack of reputation). My goal is to create a synthetic dataset that is essentially `y(x) = f(x) + c_2`, where `c_2` are random values from the same distribution as `c`. What distribution should I use for this? EDIT: Here's some more information. My data is a time-series that are essentially the number of times an action occurs in a given time bin (Poisson- like). So I have data `y(t)`, but I have noticed that I can form pretty good estimates from another variable `x`, such that I can form good models using `y(x) = f(x)`. The model is not perfect, however, hence the error term `c`. My goal is to find the distribution of `c` so that I can create an estimated new dataset with an underlying `f(x)` model with the same distribution of noise as the original dataset.


title: Constructing Random Forests for binary classification by minimizing entropy

text: I'm looking to perform a binary classification using random forests, but I do not quite understand how to minimize the entropy of the data / what tests I should run on the nodes to do so. I'm fairly new to the field, so I'd greatly appreciate any pointers or resources that could lead me in the correct direction. Thanks!


title: Changing sign estimate manually

text: validated readers, Is someone allowed to manually change sign of an estimate (obtained through a OLS), if this is supported by an underlying theory? My first idea was that this is basically a manipulation, which forces the estimates to go in the direction we want to; in fact, the sign of the estimate should directly capture the "right" direction of the effect. However, I am starting to think that it might be justified in some cases. My case is the following. My assumption is that a covariate X has a negative effect on Y, however X is represented by a proxy X1 (which is what eventually is included in the model; X does not appear in the regression). So, this is like saying that X is a function of X1. X1 can either have a negative or a positive effect on X. Thus, the sign of the estimate of B1 is negative (as expected since it proxies the effect of a variable with a negative impact on Y), when the impact of X1 on X is also negative, but it is positive when the impact of X1 on X is positive. In the latter situation, can I put a - in front of B1 and justify in the way I did now (yet more formally)? * * * I should correct myself. We have not changed the sign of the coefficient, but we changed sign of the observed values for this variable. We have experimental data, and we use proxies for competition on the demand side of the labour market (our variable of interest); we know from the previous literature that the competition has a negative effect on the dependent variable. 2 proxies we used are negatively related to Y, and thus also on competition (the proxy decreases in value, the competition increases, Y decreases). The third proxie we use is a rate and moves instead on the same direction of competition (increase in the rate, increase in the competition), therefore it has a positive effect on Y, which makes it a bad proxy. So, I thought we might use the opposite of this variable to be a proxy for competition, so that we keep the negative relationship between the proxy and the dependent variable. The opposite of this rate (by that I mean: 1-rate) is a measure that actually makes sense.


title: How many measurements (10 million molecules each time) are needed to analyze 15 million unique molecues in a 20 million pool?

text: This is my homework: There are 20 million DNA molecules in a library for high throughput DNA sequencing, each sequencing run can generate 10 million reads (i.e., analyze 10 millions of DNA molecules), and there is variation in the selection probabilities across each DNA molecule, which follows a negative binomial distribution (p=0.5, r=0.5). The question is how many sequencing run are required to have 15 million unique DNA molecules analyzed? Can anyone help me out? Thanks a lot!


title: Accelerated Failure Time scale parameter excessively large when data is right censored

text: I am using an accelerated failure time model with the Weibull distribution to predict failure times. My failure times range from 1 - 365, with many (80%) data points that are right censored (no observed failure). Among the observations that did have a failure, the median time to failure was 90 days. When you include the machines that make it all the way through the year, the distribution of failure time is U shaped - lots of mass between 1-90 days, and lots of mass for the machines that last the year. According to survreg I can get the scale and shape parameters for the Weibull distribution from a survreg fit. The scale parameter is equal to `exp(intercept term)` from survreg, which in my case is exp(9.3) ~ 11k. This seems excessively large, for example               hist(rweibull(100000, scale = 11000, shape = 0.5))      shows failure times that are extremely large - the median of that histogram is 5000 days, nearly 14 years! I understand that for most of the data we could not observe a failure, so the estimated time until failure would be somewhat large, but this is a bit unreasonable for prediction. Is there anything I can do with the model or the data, or is this phenomenon because 80% of the data is right censored and we don't know any failure times beyond a year?


title: Building the dataset for Random Forest training procedure

text: I should use the bagging ( _bootstrap aggregating_ ) technique in order to train a random forest classifier. I read here the description of this learning technique, but I have not figured out how I initially organize the dataset. Currently I first load all the positive examples and immediately after the negative ones. Moreover, positive examples are less than half of the negative ones, so by sampling from the dataset uniformly, the probability of obtaining a negative example is greater than that of obtaining a positive example. How should I build the initial dataset? Should I shuffle the initial dataset containing positive and negative examples?


title: Do mean, variance and median exist for a continuous random variable with continuous PDF over the real axis and a well defined CDF?

text: For a continuous random variable with continuous PDF over the real axis and well defined CDF, are the mean, variance, and median always well defined? Mean and variance do not always exist, e.g. for a Cauchy random variable. But what about the median? Since the the median $m$ satisfies $$F_X(m)=\int\limits_{-\infty}^{m} f_X(x)dx=0.5,$$ can I say there must be a value $m$ satisfying the condition since it is continuous random variable with well defined CDF?


title: Discretizing Exponential distribution

text: I am working on a packet generator that can generate packets of 50 different sizes and the packet sizes follow exponential distribution. Given a mean packet size how can I choose the 50 packet sizes from the continuous exponential distribution so that the resulting discrete distribution accurately follows the continuous one.Many thanks in advance.


title: How to analyse three way ANOVA that also has pre-post and one within subjects factor?

text: I'm currently having a problem understanding on how to enter and calculate my data in SPSS. I have a three way design in which there are 2 variables each in   * Experts vs Novice;    * Guided Discovery vs Control;    * Undisguised vs Disguised I've exposed my participants to a pre- and post test but do not understand how to enter the 3rd variable correctly as I've exposed undisguised and disguised together to all my participants. Could someone please clarify?


title: SVM kernel parameters?

text: Hayy all,, Im going to do classification using SVM.   As I understand we have to project our data into higher dimensional by using kernel. And there are 4 common use kernel (linear, RBF, polynomial, and sigmoid).   1. What is the different between those kernel?   2. We need to find the linear separable line / hyperplane to classify our data. And to do that, we have to chose large margin to avoid any overfitting? What is overfitting?   3. As u can see in this image, there is a w. And in many reference, they said w is a line? What is that actually? And what is that lamda??   4. What is SMO mean in SVM?


title: how to solve this probability question

text: I have a question regarding how to solve this problem: There are two predictors A and B. If A is positive, there is 60% chance of raining. If B is positive, there is 60% chance of raining. A and B is independent decision maker. What's the chance of raining when both A and B are positive. How to think about this problem in the bayesian framework? I guess there is not enough condition to solve this problem, but can we build some sort of confidence interval? thanks


title: Group Classification for uni-variate Data R

text: Given that I have some data for prices, and i have several prices groups, say,               Data = [$15,$20,$20,$10,$15,$25,$20,$15]      and Groups:               Below $5         >= $5,  < $10         >= $10, < $15         >= $15, < $20         >= $20, < $25     ...      and I want to have a method to find out which group is the best to assign my data to. Can anyone gives me some suggestion? Thanks


title: Trap 66 in WinBUGS in a hierarchical Bayesian modeling

text: I want to analyze a multilevel multidimensional model in WinBUGS. the model is as below (N=2362 students responding to K=45 items of a test, students are nested within J=116 schools):               model{     #responses     for(i in 1:N){         for(j in 1:K){             logit(p[i,j])<- a1[j]*th[i,1]+a2[j]*th[i,2]-b[j]             y[i,j]~dbern(p[i,j] )         }         th[i,1:2]~dmnorm(mu[sc[i],1:2],tau.p[1:2,1:2])     }     #school level     for(j in 1:J){           mu[j,1:2]~dmnorm(m[j,1:2],tau.s[1:2,1:2])     }              #priors     for(j in 1:J){         m[j,1:2]~dmnorm(m0[1:2],cov[1:2,1:2])     }          tau.p[1:2,1:2]~dwish(cov[1:2,1:2],2)     tau.s[1:2,1:2]~dwish(cov[1:2,1:2],2)     sigma.p[1:2,1:2]<-inverse(tau.p[,])     sigma.s[1:2,1:2]<-inverse(tau.s[,])     s2p<-sum(sigma.p[,])     s2s<-sum(sigma.s[,])     rho<-(s2s)/(s2s+s2p)          a1[1]~dlnorm(0,4)     a2[1]<-0     b[1]~dnorm(0,1)     for(s in 2:K) {         a1[s]~dlnorm(0,4)         a2[s]~dlnorm(0,4)         b[s]~dnorm(0,1)     }         }      I've set these functions as initial values:               ini<-function(){     list(tau.p=matrix(rgamma(4,100,100),2,2),     tau.s=matrix(rgamma(4,100,100),2,2),     th=rmvnorm(N,mean=c(0,0),sigma=diag(2)),     m=rmvnorm(J,mean=c(0,0),sigma=diag(2)),     mu=rmvnorm(J,mean=c(0,0),sigma=diag(2)),     a1=rlnorm(K,0, 0.4),     a2=c(NA,rlnorm(K-1,0, 0.4)),     b=rnorm(45,0,0.5))     }      I use rube package in R to check and run my analysis and everything looks fine. When I run the model I receive "Trap 66 (postcondition violated)" or "undefined real result". I think the problem is from the initials but I have no idea how to solve it. Any idea?


title: Fisher's Exact Test to assess the significance of a difference between false positive rates

text: I have trained two binary classification models on the same data and evaluated them using the same test set. For each model I have calculated a false positive rate (and a count of false positives) at the threshold which produces the maximum sum of sensitivity and specificity. The models are both random forests although they use different feature sets. I am interested statistically evaluating the difference in the two false positive rates and I was wondering if it would be valid to use a Fisher's Exact Test on a 2x2 table containing the counts of false positives and non- false positives for each algorithm? Or more generally on the entire contingency table of classification results? I know that the paper by Dietterich (1998) and several other questions have suggested the use of McNemar's test, but I have not seen any suggestions regarding Fisher's Exact Test. Note: I have previous performed a DeLong test comparing of the area under the ROC curve but I am also interested in the question of whether one classifier is statistically significantly different at a given threshold.


title: Anyone know anything about the Bivariate Granger Causality Test?

text: I am a beginner trying to conduct a granger causality test using the free online calculator http://www.wessa.net/ I have no idea what any of it means to be honest, If anyone could point me in the right direction that would be much appreciated. I am looking for answers to things such as what values to put into the Box-Cox transformation parameter, Degree (d) of non-seasonal differencing, Seasonal Period and so on. I am trying to conduct a causality test on relationship between gross domestic product and foreign direct investment inflows, not sure if this is relevant. Anything will be appreciated, thank you.


title: Apply double non-negative constrain on Sparse PCA

text: In Non-negative Sparse PCA, we apply a non-negative constrain in the coordinate matrix. Here I'm up to apply non-negative constrains on both the basis matrix and the coordinate matrix. I'm wondering if there's any research on this. By the way, I'm using nsprcomp package in R and can I apply this extra constrain using this package?


title: What is the difference between period prevalence and incidence rate?

text: I am confused about the difference between period prevalence and incidence rate. Following are from Wikipedia: > The incidence rate is the number of new cases per population in a given time > period. > > Period prevalence is the proportion of the population with a given disease > or condition over a specific period of time. Thanks!


title: Is it possible to have an estimator that is unbiased and bounded?

text: I have a parameter $\theta$ which lies between $[0,1]$. Let us say that I can run an experiment and obtain $\hat{\theta} = \theta + w$, where $w$ is a standard Gaussian. What I need is an estimate of $\theta$ which is 1) unbiased 2) almost surely bounded. Requirement (2) is crucial for me. The natural think to do is to construct a new estimator setting $\hat{\theta}$ to $1$ if it is above $1$ and to $0$ if it is below $0$. But then the estimator will not be unbiased. So what should I do? Formally, the question is whether there exists a function $f: \mathbb{R} \rightarrow \mathbb{R}$ such that $f(\hat \theta)$ satisfies (1) and (2) above. Moreover, would the situation be different if I drew more than a single sample?


title: Inverse probability weighted estimation of censored longitudinal data in R?

text: I'm looking for a package for estimating the mean of a longitudinal response under monotone dropout using a "state-of-the-art" method based on GEEs, in `R` (AIPW methods, doubly robust, whatever you want to call it). I've been Googling pretty hard and haven't come up with anything useful. So we are all on the same page, I have a sequence of censored, continuous random vectors $\mathbf Y_1, ..., \mathbf Y_n$ where $\mathbf Y_i = (Y_{i1}, ..., Y_{iJ})$ and am interested in estimating the mean response at time $J$. The data is subject to dropout, with dropout time $S_i$. Under MAR (missing at random) I think I can estimate $\mu_J = E(Y_J)$ as $$ \frac 1 n \sum_{i= 1} ^ n \frac{I(S_i = J) Y_{iJ}}{\pi(\mathbf Y_i)} \tag{$*$} $$ where $\pi(\mathbf Y_i) = P(S_i = J | \mathbf Y_i)$ is the probability of making it to the end of the study, given $\mathbf Y_i$. I can assume a model for $\pi(\mathbf Y_i)$ (probably in terms of the hazard function) and under MAR I can estimate it pretty well in principle. So I could model the hazards using the observed data to estimate the weights and then plug them into that equation. This estimator is not state-of-the-art as far as I know - it can have problems with the weights being too large and (I think) it has other problems as well. There's a large literature on estimating equation based methods that solve augmented estimating equations, which is in the same spirit but fixes the problems with $(*)$. The problem is I can't find any software to solve this in `R`. I've seen stuff aimed at MSMs (marginal structure models) and causal inference, but nothing specifically for longitudinal missing data with dropout.


title: Prepare a sampling plan using probability sampling techniques

text: Yesterday I have done my MBA exam, I have faced a question in Probability, I will put the questions please answer anyone knows? > A research agency is planning to conduct a survey regarding "Students > Attitude towards facebook". They want to collect data from engineering > students across Tamilnadu. The sample size is 1000. Prepare a Sampling Plan > using Probability Sampling Techniques.


title: Estimation of covariance and mixing parameter by a bivariate normal-lognormal model

text: I have to create a model which is a mixture of a normal and log-normal distribution. To create it, I need to estimate the 2 covariance matrixes and the mixing parameter (total=7 parameters) by maximizing the log-likelihood function. This maximization has to be performed by the `nlm` routine. As I use relative data, the means are known and equal to 1. I’ve already tried to do it in 1 dimension (with 1 set of relative data) and it works well. However, when I introduce the 2nd set of relative data I get illogical results for the correlation and a lot of warnings messages. To estimates the parameters I defined first the log-likelihood function with the 2 commands `dmvnorm` and `dlnorm.plus`. Then I assign starting values of the parameters and finally I use the nlm routine to estimate the parameters (see script below).               P <- read.ascii.grid("d:/Documents/JOINT_FREQUENCY/grid_E727_P-3000.asc", return.header= FALSE )           V <- read.ascii.grid("d:/Documents/JOINT_FREQUENCY/grid_E727_V-3000.asc", return.header= FALSE )           p <- c(P) # tranform matrix into a vector     v <- c(V)          p<- p[!is.na(p)] # removing NA values          v<- v[!is.na(v)]          p_rel <- p/mean(p) #Transforming the data to relative values          v_rel <- v/mean(v)           PV <- cbind(p_rel, v_rel) # create a matrix of vectors          L <- function(par,p_rel,v_rel) {          return (-sum(log( (1- par[7])*dmvnorm(PV, mean=c(1,1), sigma= matrix(c(par[1], par[1]*par[2]*par[3],par[1]*par[2]*par[3], par[2] ),nrow=2, ncol=2))+     par[7]*dlnorm.rplus(PV, meanlog=c(1,1), varlog= matrix(c(par[4],par[4]*par[5]*par[6],par[4]*par[5]*par[6],par[5]), nrow=2,ncol=2))            )))          }     par.start<- c(0.74, 0.66 ,0.40, 1.4, 1.2, 0.4, 0.5) # log-likelihood estimators          result<-nlm(L,par.start,v_rel=v_rel,p_rel=p_rel, hessian=TRUE, iterlim=200, check.analyticals= TRUE)          Il y a eu 50 avis ou plus (utilisez warnings() pour voir les 50 premiers)          1: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) : production de NaN          2: In sqrt(2 * pi * det(varlog)) : production de NaN          3: In nlm(L, par.start, v_rel = v_rel, p_rel = p_rel, hessian = TRUE,  ... :  NA/Inf           replaced by maximum positive value          4: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) : production de NaN          5: In sqrt(2 * pi * det(varlog)) : production de NaN          6: In nlm(L, par.start, v_rel = v_rel, p_rel = p_rel, hessian = TRUE,  ... :  NA/Inf replaced by maximum positive value      ...... until 50               par.hat <- result$estimate          cat("sigN_p =", par[1],"\n","sigN_v =", par[2],"\n","rhoN =", par[3],"\n","sigLN_p =", par[4],"\n","sigLN_v =", par[5],"\n","rhoLN =", par[6],"\n","mixing parameter =", par[7],"\n")           sigN_p = 0.2919377            sigN_v = 0.4445056            rhoN = 1.737904            sigLN_p = 2.911735            sigLN_v = 2.539405            rhoLN = 0.3580525            mixing parameter = 0.8112917      Does someone know what is wrong in my model or how should I do to find these parameters in 2 dimensions? Thank you very much for taking time to look at my questions.


title: libsvm data format

text: I'm using the libsvm (http://www.csie.ntu.edu.tw/~cjlin/libsvm/) tool for support vector classification. However, I'm confused about the format of the input data. From the README: > The format of training and testing data file is: >   >   >     <label> <index1>:<value1> <index2>:<value2> ... >     . >     . >     . >   > > Each line contains an instance and is ended by a '\n' character. For > classification, `<label>` is an integer indicating the class label (multi- > class is supported). For regression, `<label>` is the target value which can > be any real number. For one-class SVM, it's not used so can be any number. > The pair `<index>:<value>` gives a feature (attribute) value: `<index>` is > an integer starting from 1 and `<value>` is a real number. The only > exception is the precomputed kernel, where `<index>` starts from 0; see the > section of precomputed kernels. Indices must be in ASCENDING order. Labels > in the testing file are only used to calculate accuracy or errors. If they > are unknown, just fill the first column with any numbers. I have the following questions:   1. What is the use of the `<index>`? What purpose does it serve?   2. Is there a correspondence between the same index values of different data instances?   3. What if I miss/skip an index in between ? I ask because the datafile *heart_scale* which is included in the package for libsvm, on line 12, the index starts from 2. Is the `<value>` for index 1 taken as unknown/missing? Note: the tools/checkdata.py tool provided with the package says that the *heart_scale* file is correct.


title: Assessing whether differences between datasets are significant

text: I am unsure which statistical test I should use to test whether two datasets are significantly different. The data looks something like this:               Group   Dataset 1   Dataset 2     A   5   6     B   1   7     C   2   10     D   1   15     E   4   11     F   1   5     G   2   17     H   2   10     I   3   19     J   1   6     K   14  67     L   13  38     M   3   16     N   14  59     O   14  33     P   0   4     Q   10  28     R   0   6     S   4   9     T   3   5     U   3   13     V   0   13      I am hoping that, although one dataset has a larger sum, the proportions between groups mean that they are not significantly different. Does anybody have any helpful advice or suggestions of a suitable test? Thanks in advance


title: Leave-one-out cross validation in selecting predictor

text: I am a newbie here. There are 155 total samples. Five different predictors Xi (i=1,2...5) are used to predict Y, like X1 X2 X3 X4 X5 Y .... The objective is to find the best predictor Xi to predict Y. Leave-one-out cross validation (LOOCV) is used. There would be 154 data in training and 1 data in test in each round. There are total 155 rounds. My question is   1. how to find the best Xi?   2. what is the last regression line? using total 155 data? or choosing one result of 155 rounds?


title: Do you need to report the standard error of the mean when the sample is the population?

text: The other day, my tutor told me that whenever a mean value is reported, the standard error of the mean must be reported alongside. But surely this isn't always the case? For example, let's say we calculate the mean number of goals a football team scores per game across a season. We're not sampling from a population in this case - the number of goals the team scores is the population. Am I right in saying the standard error would not need to be reported in this case?


title: A modelling question about point processes with heavy tails

text: I am trying to model a number of point processes for which I have data. If I choose to model each one using a (different) homogeneous Poisson process and estimate the rate using MLE then for some of the processes the goodness of fit is quite good. However for others it isn't and when I plot a histogram of the interarrival times, I see that they look exponential at first but then have a long tail that the Poisson process doesn't explain. Naively I would like to try a model where the interarrival times are heavy tailed. However I understand this is no longer a Poisson process but rather a renewal process now. Is there some way to do MLE (or similar) to estimate the parameters of a renewal process so I can then do a goodness of fit test to see how plausible the model is?


title: Can Poisson regression be run with a small sample?

text: I am just a beginner. I have a sample of 71. I wish to do a Poisson regression. Is this sufficient for estimating a Poisson model? Should I use only counts as $Y$ variable or should I use percentages (i.e., counts converted into percentages)? For the given sample size up to how many $X$ variables could I use in a single model?


title: Is there a better test than chi-2 to determine plausibility of a distribution?

text: So, I'm trying to determine what's the best distribution to use to analyze my data. Previous authors suggest using a chi-2 test to determine whether the data is Poisson or Gaussian distributed. All approaches are empirical and there's no consensus on how the data should be distributed. It can be either, continuous if you use newer programs, or discrete if you use the standard software for the field; not updated since the early 90's. The measurements can only be be recorded as integers initially. The portion that I'm working on is considered the same data point measured 2075 times. Each of the 2075 measurements is an average of that pixel across around 300 images. The problem comes from how to find an uncertainty for that value. The uncertainty is necessary since later data analysis programs need an estimate. The data is curved because it should be due to the nature of the experiment. There's a lot of junk; a small population of very high values from two erroneous sources, a minor problem with size limits on integer values in the old program that redistributes a very small number of intensities, there is a background and an attempt to remove the background. The biggest problem is the very high intensities because they dramatically affect mean. So we needed to find better ways to target them for removal and a better estimate for uncertainty. Older lit suggests just uses poisson and newer suggests gaussian. The newest stuff wants to use a normal-Pareto. Is there a more powerful test to fit the dataset than a chi-square comparison to a gaussian or poisson fit? Preferably something included in matlab.


title: Reliability estimate for each item (comparing raters to each other)

text: I am working with a 15 item scale of pre-school teaching. Five observers observed 58 teachers. Its a complete crossed matrix i.e, all raters observed all teachers. In an attempt to weed out the weak raters, I would like to calculate reliability estimate for each item, which compares raters to each other. Is there a way I can do that through SPSS or R? danke!


title: What is this distribution?

text: Does anyone know what the following distribution is, please? I can tell it is discrete and $\theta$ is somehow a probability. Thank you! $$\mathbb P_X(x; \theta) = (x-1) \theta^2 (1-\theta)^{x-2}, \ \ x=2,3,\ldots; \ \ 0 < \theta < 1.$$


title: Detecting outliers to achieve better performance for neural network

text: I'm trying to find the ways/algorithm that detects outliers to build better (classification) model by using neural network. According to numbers of research papers, I've found that outliers significantly effects NN modelling. Yet, I couldn't find a method to detect them. Does anyone know the method? (it would be great if you could kindly suggest R package!) Thanks!


title: 2 x 2 mixed anova?

text: 96 participants were randomly assigned to one of the two groups regardless of their motivation level. Participants in both the groups were exposed to both the experimental conditions of gain framed messages and loss framed messages. The only difference was - counterbalancing such that the participants were shown gain frame message in group 1 then its corresponding loss frame message was shown in group 2. So here the within group factor is: framing condition - gain frame and loss frame between group factor means that this should vary across two groups but in my study both participants of high and low motivation were present in both groups and exposed to both gain frame and low frame messages. Participants in two groups were different but exposed to same but counterbalanced design. I am not sure if this is a 2 x 2 mixed anova design with framing condition as within group and motivation as between group factor. Please help.


title: SPSS: Crosstabs multiple comparison

text: I have a table with 11 columns and 2 rows. Column 1 is the control group and I want to compare all the other columns with it. I'm not interested in compraing columns 2 - 11 with each other and I'm concerned that Bonferroni's correction (which seems to be the only option available) is overcompensating by multiplying the p value by 55 or however many comparisons it is. Is there something analogous to Dunnett's so that I can compare just with the control group? Or should I just do 10 chi-squared tests and multiply each p value by 10? Also, when I do use Bonferroni'sin a chi-squared, is it possible to see the exact p values? Thank you


title: Is " independent and identically distributed" an assumption or a fact ?

text: This is in the context of two random variables. A frequent assumption (e.g. of the error term in ANOVA) is of independent and identically distributed random variables. There is a question on this site asking how the assumption can be checked in a given dataset. Is it an assumption or a fact?


title: Is there any procedure to test the significance of meta-analytic estimate of correlation based on a set of sample-correlations?

text: Meta-analysis generates an estimate of an overall correlation. How could we examine that it is a valid estimate.


title: Clustering University Courses using Machine Learning

text: I have a database with 32344 Courses from Swedish universities. A course have the following attributes:               private String code;     private String name;     private String description;     private String localName;     private List<CourseRank> ranks;     private List<Comment> comments;     private List<UserBook> userBooks;     private List<CourseBook> courseBooks;     private University university;     private String teacherName;     private String courseLink;     private String requirementsLink;     private String localCourseLink;     private String localRequirementsLink;     private String level;     private String period;     private String timeOfDay;     private String type;     private int pace;     private int numberOfAssignments;     private int numberOfProjects;     private int numberOfExams;     private double credits;      Most of the courses only have name, university, level, period, pace, type and credits filled in but the goal is that most courses will have everything filled in. Each course have a list with ranks, which basically are course evaluations from students. A rank looks like this:               private int workload;     private int structure;     private int teachers;     private int avarageScore;     private String comment;     private String date;     private int userId;      I want to apply some kind of machine learning to this database. I am thinking of using clustering algorithms to categorize similar courses to each other. What I had in mind was something like this:   1. Let a clustering algorithm compare name and description of a course to show similar hits when a user is searching on the website.   2. Or use all attributes of a course, including course ranks to classify courses into clusters. This could help students find courses they like by suggesting courses that are similar to the one they have marked as favorite or given a good rank. Does this make sense and what algorithms would you recommend? Any other suggestions on what kind of fun machine learning I could apply here? Best regards, Nilsi


title: How to select discrete cosine transform coefficients as a feature vector?

text: I need to use DCT on frames of videos as a feature vector to train a Feed- forward artificial neural network, but the problem is the large number of coefficients (in thousands or so). How can I choose some of these coefficients and still capture the majority of information for my image (frame)?


title: Generate a probability distribution from a histogram with two peaks in matlab

text: I have a histogram with two peaks and I want to generate the corresponding probability distribution. I have used the following MATLAB code:               A=mydata;     M1=max(A);     M2=min(A);     I=(0:100).*(M1-M2)./100+M2;     [n,x]=hist(A,I);     bar(x,n/(1000*0.352))      I found this code frequently to explain how we can find a prob distribution for a histogram of random numbers normally distributed but I don't know whether if this true for a histogram with two peaks and generate a normalised probability distribution. The histogram is shown here: ![The bimodal histogram](http://i.stack.imgur.com/wSxjf.png)


title: How to expand data frame in R

text: I'm having following problem while doing some analysis with R. I have a dataframe like this:               Name | Group | Count     Person 1 | A | 3     Person 2 | A | 1     Person 3 | A | 0     Person 1 | B | 5      Person 2 | B | 0     Person 3 | B | 1     Person 1 | C | 1      and I'd need to "expand" it (not sure if the right term) to be like this:               Person 1 | A     Person 1 | A     Person 1 | A     Person 2 | A     Person 1 | B     Person 1 | B      etc. So it takes value of pair Person 1 and A (in this example, 3) and makes three rows with Person 1 and A and does so for every Person - Group -combination. Can't figure out any good words for searching online.


title: exact definition of rate of convergence

text: Can anyone concisely clarify within the framework of parametric regressions and sample size $N$ what is the rate of convergence? Increasing the exponent $k$ in $N^k$ (or $N^{-k}$) will increase or decrease this rate and why is this considered to be good (or bad)? The reason why I am asking is that when googling "convergence rate" I find both sources that claim the rate to be $N^{0.5}$, while other sources claim it is $N^{-0.5}$. Both can't be right.


title: Zero-inflated negative binomial model for true zeros

text: The `zeroinfl` function in the `pscl` package in R assumes that zeros include both false zeros and true zeros. I have a zero inflated dataset that I'm confident only includes true zeros. Is it appropriate to use `zeroinfl` function on this dataset? Or would the `hurdle` function be more appropriate?


title: Output of logistic model in R

text: I'm trying to interpret the following type of logistic model:               mdl <- glm(c(suc,fail) ~ fac1 + fac2, data=df, family=binomial)      Is the output of `predict(mdl)` the expected odds of success for each data point? Is there a simple way to tabulate the odds for each factor level of the model, rather than all the data points?


title: Relationship between number of replications and the standard deviation of entry means

text: I'm working at a plant breeding company. During some spare time I recently did some exploring of our data (yield trials), mostly just out of curiosity and wanting to learn R better. I started by plotting standard deviation of means created with different number of replications (in the actual data each color would represent a particular seed lot grown in a particular environment). Like this (examples are not actual yield data):   ![Figure 1](http://i.stack.imgur.com/cF3lz.png) Figure 1. Standard deviation of means created from different number of normally distributed numbers. That looks perfectly understandable to me. I then got the idea to plot the standard deviation relative to the 1-rep standard deviation: ![Figure 2](http://i.stack.imgur.com/ASjm6.png) Figure 2. Standard deviation, relative to standard deviation in 1 rep, of means created from different number of normally distributed numbers.2 Obviously this relative standard deviation has a strong relationship with the number of replications. I have been thinking about it a while, but I'm having trouble to formulate the math behind this. So my question is this: **Can someone please explain to me, in mathematical terms, why we see the relationship visualized in the second figure?** (I have a feeling that the answer will leave me feeling utterly stupid for asking. I also realize that this may sound an awful lot like I'm trying to get you to do a homework assignment for me, but I'm really just curious and have access to loads of data but not a lot of training in statistics and mathematics. We have a statistician that I could have tried to ask, but he's on vacation.)


title: Factor analysis: What to do when determinant is almost zero and when KMO for a variable is low?

text: I'm conducting a factor analysis on 40 interval-level variables, and have two immediate concerns:   1. The determinant is `6.608E-006`, which is much lower than the cut-off of `0.00001`. I went back and screened the correlation matrix to find significant, too-high correlations between variables; there's nothing even approaching `0.8`. What now? Do I proceed anyway, just noting the determinant value, or is there something else I should be doing to adjust the data to exclude multicollinearity?   2. KMO is a respectable `.846`, but on examining the anti-correlation matrix, I found five variables whose individual KMO values were less than `0.5`. Now what? Exclude? Ignore, as overall KMO is high enough already?


title: Can I do parallel analysis with any type of exploratory factor analysis/principal component analysis?

text: I wish to perform parallel analysis to determine how many factors I should extract from my maximum likelihood exploratory factor analysis. I have been referred to a program that calculates the eigenvalues for random data using Monte Carlo for principal component analysis. I am not doing principal component analysis, however. I am doing maximum likelihood exploratory factor analysis. I have been told that you can do it for any type of EFA, but I am uncertain. For example, a few of the macros I have seen require that you identify if you are using PCA or PAF. This clearly means that it matters to some extent. Note that I have tried doing parallel analysis for PAF using the engine at http://ires.ku.edu/~smishra/parallelengine.htm. Unfortunately, the results did not make much sense. Every eigenvalue in my data was greater than the eigenvalue from the website, which I believe meant I should extra 100+ factors. Questions:   * Can I use the parallel analysis results designed to be used with PCA to determine how many factors I should extract from my maximum likelihood factor analysis?   * If not (I expect the answer to the above question is "no"), how can I do a parallel analysis to determine how many factors I should extract from my maximum likelihood exploratory factor analysis?


title: Can any probability distribution be written as a Boltzmann distribution?

text: By redefining the energy function, $E(x)$, can any $p(x)$, be written as a boltzmann distribution, ie. $p(x) = \frac{e^{-E(x)}}{Z}$, where Z is the partition function?


title: How to explain usage pattern that follows heavy-tail distribution to non-experts?

text: I need to measure and communicate to non-statisticians the usage pattern of a software component. It seems that based on data the usage follows a nice smooth heavy-tail distribution: lots of users don't use the component at all, few users use it a lot. I don't want to summarise my report using a single mean, since I believe it hides the true nature of the data and non-experts will easily interpret it as it was a normal distribution. On the other hand I can't plot the distribution and start explaining what is power-law, why it doesn't look like normal distribution, why "outliers" are not really outliers but well expected users etc, etc... So I ended up splitting the distribution in 3 parts: users that make some small use of the component, average users, and power users (that use component a lot). For each group I measured mean values (and error of the mean) and in my opinion reporting these few numbers is something more actionable and easier to digest by non-experts. I wonder if there is some better way to perform that kind of analysis for the non experts. Any suggestions more than welcome!


title: How many times should we repeat a K-fold CV?

text: I came across this thread looking at the differences between bootstrapping and cross validation - great answer and references by the way. What I am wondering now is, if I was to perform repeated 10-fold CV say to calculate a classifier's accuracy, how many times _n_ should I repeat it? Does _n_ depend on the number of folds? On the sample size? Is there any rule for this? (In my case, I have samples as big as 5000, and if I choose anything larger than _n = 20_ my computer takes way too long to perform the calculation. )


title: How can I properly relate the dispersion parameter to other estimated variance parameters?

text: In a negative binomial model, could I treat the dispersion parameter as an estimate of the residual variance? In other words, could I do the following with the dispersion $D$: $$ \sigma^2 = \ln(\mu + \frac{\mu^2}{D}) $$ And use $\sigma^2$ along with other variance estimates (e.g., random slope or random intercept variability) to calculate proportion of explained variance of some particular effect? If not, what is the best way to treat the dispersion estimate in relation to other estimated variance components?


title: Why do I have 0 degrees of freedom for my group mean variables in my mixed model? nlme package

text: I am trying to run a mixed effects model in the R package 'nlme'. I have 3 groups and 294 observations, and am using the following model:                model1<-lme(shy ~ int + ind + grpint +grpind, random =~1 | group,                   data = dt, control=list(opt="nlmimb"))      The variables 'grpint' and 'grpind' are vectors of group means for the 'int' and 'ind' variables. (I am running this test to see if individual scores in 'int' or 'ind' still significantly predict group differences in 'shy' over and above group scores.) When looking at my output (below), I am very confused as to why my group variables have 0 df. Did I enter something incorrectly? Is this still acceptable when attempting to test my research question? Thank you so much for sharing your knowledge and expertise!                Linear mixed-effects model fit by REML       Data: dt              AIC      BIC    logLik        93.80127 119.4663 -39.90063           Random effects:      Formula: ~1 | group              (Intercept)  Residual     StdDev:  0.07421543 0.2755107           Fixed effects: shyness ~ interdependent + independent + grpint + grpind                           Value Std.Error  DF    t-value p-value      (Intercept)     2.5677638  37.76813 289  0.0679876  0.9458      interdependent -0.0523904   0.05406 289 -0.9691753  0.3333      independent     0.1623483   0.05329 289  3.0462865  0.0025      grpint          1.5382962  23.23831   0  0.0661965     NaN      grpind         -1.6347691   7.78430   0 -0.2100085     NaN       Correlation:                      (Intr) intrdp indpnd grpint      interdependent  0.000                           independent     0.000 -0.061                    grpint         -0.998 -0.002  0.000             grpind          0.979  0.000 -0.007 -0.991           Standardized Within-Group Residuals:              Min          Q1         Med          Q3         Max       -4.39909336 -0.43557519 -0.09616286  0.48480368  3.05524916            Number of Observations: 294      Number of Groups: 3


title: Desirable and undesirable properties of Latin squares in experiments?

text: A cursory search reveals that Latin squares are fairly extensively used in the design of experiments. During my PhD, I have studied various theoretical properties of Latin squares (from a combinatorics point-of-view), but do not have a deep understanding of what is it about Latin squares that make them particularly well-suited to experimental design. I understand that Latin squares are good at allowing statisticians to efficiently study situations where there are two factors which vary in different "directions". But, I'm also fairly confident there would be many other techniques that could be used. > What is it, in particular, about Latin squares that make them so well suited > for the design of experiments, that other designs do not have? Moreover, there are zillions of Latin squares to choose from, so which Latin square do you choose? I understand that choosing one at random is important, but there would still be some Latin squares that would be less suited to running experiments than others (e.g. the Cayley table of a cyclic group). This raises the following question. > Which properties of Latin squares are desirable and which properties of > Latin squares are undesirable for experimental design?


title: LOESS that allows discontinuities

text: * Is there a modelling technique like LOESS that allows for zero, one, or more discontinuities, where the timing of the discontinuities are not known apriori?   * If a technique exists, is there an existing implementation in R?


title: When does rules based classifier outperforms decision trees?

text: I am trying to identify approximate 3% of the population for some characteristic feature. Standard decision tree or logistic regression gives too many false positives. Is there a chance that rules based classifier can improve performance? I would like to get approx 75% of recall with 95% of precission (i.e. FalsePositives <= 5% of Positives)


title: Test for significance of unexpected non-linear trend

text: How can I test whether there is a significant non-linear effect of a continuous predictor on a continuous DV? Here's my situation. I ran a simple linear regression and found a significant effect of a predictor (predictor2 below). However, when I plot my dv against that predictor, the trend does not look linear - in fact it looks like the dv has two local maxima. When I submitted these findings to a journal, reviewers commented that I should explore the possible statistical significance of this apparent non-linear effect. The effect was not expected and so I do not have a particular non-linear model to test here. I am just looking for an appropriate off-the-shelf technique of first call to use in a situation like this. Ideas? For reference, my data look like this:               > D[1:10,]        sid      posttest      pretest       predictor2      predictor3       predictor4        predictor5     1    2          0.75         0.75        0.8823529      1.00000000        0.8888889                18     2    3          1.00         1.00        0.8666667      0.06666667        1.0000000                16     3    5          1.00         1.00        0.9000000      1.00000000        0.9047619                21     4    6          0.75         0.75        0.3333333      0.93333333        0.9375000                16     5    7          0.50         0.75        0.8750000      0.87500000        0.8333333                18     6    8          0.75         1.00        0.5000000      0.50000000        0.8823529                17     7    9          1.00         0.50        0.8823529      0.58823529        0.4444444                18     8   10          1.00         0.50        0.8750000      0.87500000        0.5882353                17     9   13          0.75         0.75        0.7777778      1.00000000        0.8947368                19     10  15          0.50         0.25        0.8666667      0.46666667        0.8125000                16     ...      The linear regression model and its output look like this:               > lmA <- lm( posttest ~ pretest + predictor2 + predictor3 + predictor4 + predictor5, data=D )     > summary( lmA )          Call:     lm(formula = posttest ~ pretest + predictor2 + predictor3 + predictor4 + predictor5, data = D)          Residuals:          Min       1Q   Median       3Q      Max      -0.68810 -0.11853  0.05159  0.13589  0.36904           Coefficients:                     Estimate  Std. Error   t value  Pr(>|t|)         (Intercept)     0.773645    0.106033     7.296  8.45e-13 ***     pretest         0.292674    0.030905     9.470   < 2e-16 ***     predictor2      0.063541    0.024557     2.588   0.00988 **      predictor3      0.005529    0.026369     0.210   0.83398         predictor4      -0.026421   0.065049    -0.406   0.68474         predictor5      -0.009412   0.005010    -1.879   0.06070 .       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.1801 on 665 degrees of freedom     Multiple R-squared:  0.1298,    Adjusted R-squared:  0.1232      F-statistic: 19.83 on 5 and 665 DF,  p-value: < 2.2e-16      Finally, here is the plot that gave rise to the reviewer comment. Note that binned data with bubbles were used, with bubble diameter indicating number of data rows in each bin, instead of a scatterplot because there are only 5 levels of the dv (0.0, .25, .5, .75, 1.0) so you can't see anything at all with a scatterplot. Reviewers commented that bins could be arbitrarily chosen to make the linear trend look more persuasive, which is true, but binned data still has to be used for the above reason I think. ![Plot of posttest against binned predictor2](http://i.stack.imgur.com/ClB0N.png)


title: Testing for the compatibility of inferences

text: Consider that I have two balances (called 1 and 2). Each of these balances gives a posterior distribution for the weight of the object of the form $m_1 \pm s_1$ (for balance 1) and $m_2 \pm s_2$ (for balance 2). I will then assume that the posterior density are normal with mean $m_i$ and std $s_i$ and that conditionnaly to the true weight, the two posteriors are independent each other. Then consider that I weigh a given object $N$ times with the two methods that gives me two sets $(m_{1,i} \pm s_{1,i})_{i=1,N}$ and $(m_{2,i} \pm s_{2,i})_{i=1,N}$. I would like to test if the obtained inferences are compatible each other at a given "confidence level" i.e. to assess if the inferences provided by two methods (not relying on the same kind of observations but on the same object) are compatible. If yes, one can thrust the results, if not one can be suspicious. How can I proceed?