title: Why direct oblimin rotation results in greater eigen values?

text: I came across this in the Wikipedia page about Factor Analysis. Is that true that _direct oblimin rotation results in greater eigen values_? If that is true, what's the reason behind it and does it generalize to other **oblique** rotations ? (to avoid any confusions consider we only consider PCA as the factor extraction methods)


title: A tool to analyze multidimensional tabular data

text: I have a bunch of data in the following format:               Count|| Fruit  |Fruit Col.|| Veget. |Veg.Col.|| Spice         ---------------------------------------------------------      123 || Apple  | Red      || Carrot | Orange || Cumin          20 || Apple  | Red      || Carrot | Green  || Pepper       ... 20K more combinations...        1 || Orange | Orange   || Potato | White  || Allspice      This is to be read as follows: **_Red apples** in combination with **orange carrots** and **cumin** were encountered **123** times during the experiment._ I need to be able to answer at least following questions about the dataset:   * What are the most popular spices for the most popular vegetables among the most popular fruits, (a) with and (b) without the regard for the color of either?   * Which color of potatoes is encountered most often along with the most popular vegetable?   * What is the most popular red fruit that is encountered with green vegetable? Etc. "Most popular" here assumes answer like:   1. 5% Apple   2. 4% Orange   3. 1% Others Please suggest a tool which will allow me to analyze my dataset in the way I described above. Ideally the tool should _not_ require me to write complicated code to ask the specific question (some coding beforehand is acceptable). A tool that will do what I need without coding would be perfect.


title: Java support for PMML

text: I am new in PMML, Predictive Model Markup Language and I was wondering if there is some kind of Java support (Open Source / professional) for creating/parsing PMML files. Initially I only have in mind the possibility of creating/parsing PMML files programatically from Java environments. I have been "googling" and I have found several possibilities:   * Open source `jpmml` from PMML 3.2.   * JDM's `javax.datamining`, however it seems dead... Does someone have more info?   * Professional library Zementis.   * Use an XML Java library and build DIY parser/writer of PMML files. I appreciate all your opinions. Thanks in advance Oscar


title: How to present dummy variables from linear regression in table?

text: I have carried out this linear regression that includes month coded as a dummy variable:               library(plyr)     set.seed(1)     y <- rnorm(120)     x1 <- c(rep("adult", 60), rep("juvenile", 60))     x2 <- c(rep("male", 60), rep("female", 60))     x3 <- unlist(llply(month.abb, function(x) rep(x, 10)))          summary(lm(y ~ x1 + x2 + x3))          Call:     lm(formula = y ~ x1 + x2 + x3)          Residuals:          Min       1Q   Median       3Q      Max      -2.46354 -0.51524 -0.03981  0.57625  1.95041           Coefficients: (2 not defined because of singularities)                 Estimate Std. Error t value Pr(>|t|)     (Intercept)  0.12073    0.28564   0.423    0.673     x1juvenile   0.00663    0.40396   0.016    0.987     x2male            NA         NA      NA       NA     x3Aug       -0.37510    0.40396  -0.929    0.355     x3Dec       -0.24718    0.40396  -0.612    0.542     x3Feb        0.12812    0.40396   0.317    0.752     x3Jan        0.01147    0.40396   0.028    0.977     x3Jul        0.32385    0.40396   0.802    0.424     x3Jun        0.02273    0.40396   0.056    0.955     x3Mar       -0.25440    0.40396  -0.630    0.530     x3May        0.01341    0.40396   0.033    0.974     x3Nov        0.22012    0.40396   0.545    0.587     x3Oct       -0.01502    0.40396  -0.037    0.970     x3Sep             NA         NA      NA       NA          Residual standard error: 0.9033 on 108 degrees of freedom     Multiple R-squared:  0.04703,   Adjusted R-squared:  -0.05003      F-statistic: 0.4845 on 11 and 108 DF,  p-value: 0.9093      I now want to present the results of this linear regression within a table. Instead of presenting the beta for every month, is there a way to summarise the overall effect of month on `y` within the same table? For example, would if be acceptable to summarise the beta, se, t value and p value of `x3` by using their mean values across months?


title: Multiple Imputation and prediction

text: hopefully someone has an idea regarding the following problem: I did a multiple imputation using the mice package - now just one step is left: I want to generate final and "valid" imputations for my NAs. How can I do that? With multiple imputation I get by default 5 datasets with different values. These I need to calculate a prediction model with statements about the standard errors. But: How do I get a final value for my NAs? Here my code:               imp <- mice(data, m = 5) # Compute 5 imputated data sets     mifit <- with(imp, lm(y ~ x1 + x2 + x3 + c1 + c2)) # Compute out of the 5 data sets from imp 5 linear models with y as dependent variable     pmifit <- pool(mifit) # pool the models to receive final coefficients       If the pmifit would be of class(pmifit) == "lm" I would use the predict function to replace my missings but pmifit is an object with no method in predict. Anyone an advise?


title: Material on Structural Equation Modelling (SEM)

text: I've little background on SEM. Now I want to focus on its theoretical aspects like statistical model, model estimation and model evaluation. The seminal book on this topic is by Bollen: Structural Equations with Latent Variables. Unfortunately, this book is not in my access in my region. I wonder if there are any online class material or articles available freely to fulfil my requirements. Thanks in advance.


title: Why do I get difference results when using pooled ols and random effects model?

text: I was wondering what the difference between pooled OLS and random effects model is? I know that random effects eliminates the time-constant effect. However, when i run my regression with pooled ols and random effects model I get different results. The pooled ols regression looks find where all the variables are significant. When i run the random effects model (with the same dependent and independent variables) many of the variables get insignificant. The sample size is smaller for the random effects model because for some observations I only have one year which I can use in my pooled OLS regression but not in my random effects regression? Could this be the reason why I get different results?


title: How to run two motion plots together in R?

text: I'm doing a real recreation of the annual tree growth from tree-ring measurements, and I have many problems working with graphs and their margins. I have done the two graphs separately and these work fine, but at the moment to put them together, it doesn't work! and if it is possible I would like also save this as GIF file. My main idea is to plot these two graphs together, and if is possible save them as a Gif file. Here my two motion graphs separately:               tree01 <- read.table(url("http://darwin.ec/test/tree01.txt"),header=T, dec =",")          x = tree01$x         y = tree01$y     z = tree01$z          year=tree01$year         TRW = tree01$TRW          # PLOTTING DISC     plot.new()     for(i in 1:125) {       Sys.sleep(0.1)       par(new=T)       plot(x[1],y[1], xlim = c(-1000, 1000), ylim = c(-10000, 10000))       par(new=T)       points(x[i],y[i], cex= z[i], col="black")       legend("topright", legend= year[i], border="") # I tried to draw without border     }          # PLOTTING LINES     plot.new()          for (i in 1:125){       Sys.sleep(0.1)       plot(x=year, y=TRW, type="l", ylim = c(0,5))        points(year[i], TRW[i], col="red", pch=20, font = 20)       abline(v=year[i], col = "grey")     }      That is why a tried to put them together, but in different panels, but unfortunately doesn't work:               library(animation)          # saveGIF({          for(i in 1:125) {       old.par <- par(mfrow=c(2,1))       #disc       Sys.sleep(0.1)       par(new=T)       plot(x[1],y[1], xlim = c(-10000, 10000), ylim = c(-10000, 10000))       par(new=T)       points(x[i],y[i], cex= z[i], col="black")       legend("topright", legend= year[i], border="")       #lines       par(new=T)       plot(x=year, y=TRW, type="l", ylim= c(0,500))       points(year[i], TRW[i], col="red", pch=20, font = 20)       abline(v=year[i], col = "grey")       par(old.par)      }          # }, movie.name = "tree01.gif", interval = 0.1, ani.width = 600, ani.height = 600)      and at the end I would like to have a graph in two panels, but running at the same time like this: ![Annual tree growth from tree-rings](http://i.stack.imgur.com/DPn6h.png) I will be very grateful for your help. Thanks!


title: Is it fair to say p-values tell us nothing about the probability null hypotheses are true?

text: I understand that NHST assumes the null hypothesis is true, and that we can't infer from a low-value that a null hypothesis is necessarily false. However, is it fair to say that the p-value tells us _nothing_ about the conditional probability $Pr(H_0 = T|T_{obs})?$ If it is fair to say that p-values tell us nothing about the probability null hypotheses are true, why does anyone use them?


title: Zero inflated negative binomial with selection

text: I am looking for a Stata (or R/Matlab if there's no Stata) implementation of the model described by Greene (1994) (http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1293115). It is essentially a Heckit for ZINB instead of OLS. Would anyone know if such an implementation exists? Also, a critical formula in this model (3.6, p18) is somehow missing from the Working Paper. If anyone knows of another version that has it, that would be greatly, appreciated.


title: Is that overfitting?

text: Given a document-term matrix $X$, where $$X(d, t) = \textit{occurrences of 't' in 'd'}$$, it's possible to compute it's Truncated Singular Value Decomposition:$$X_k = U_k \Sigma_k V_k^T$$ Then, for a document profile $d$ (i.e. the vector of occurences in that document), it is possible to compute a 'predicted' document profile $\hat{d}$ that, according to the model, is more consistent, as follows: $$\hat{d} = d V_k V_k^T$$ This is the LSA/LSI/SVD method. Now, suppose that instead computing only $V_k$, I compute $C$ matrices $V_{k,c}, c = 0, 1, \dots ,C $, with the same dimensions of $V_k$ but different entries. Then, for a given document profile, I compute $C$ 'predicted' profiles, as: $$\hat{d}_c = d V_{k,c} V_{k,c}^T$$ one for each $c = 0, 1, \dots ,C $. Among these, I choose as final predicted profile, the one which minimize the L-2 norm, w.r.t. the input profile: $$\hat{d} = \arg \min_c \| d-\hat{d}_c\|_2 $$. I'm not totally persuaded by this algorithm; choosing the output with a L-2 norm w.r.t. the input seems to me to be source of overfitting. Can I have to your opinion, suggestion, critiques please? Thank you very much.


title: significance test to use for multiple conversions in an A/B test

text: What is the proper way to measure significance in an A/B test, when a visitor can potentially "convert" more than once? I understand a chi-square test or z-test are most often used for A/B testing, but from what I understand that is only if it is categorical, like if a visitor purchased or not. In my case, customers can purchase multiple times in a test period, and I want to make sure I'm using the right significance test for this. For example:               test A: 800,000 visitors, 50,000 purchases     test B: 45,000 visitors, 2,700 purchases      I saw a similar question here, but it was asking about revenue: AB Testing other factors besides conversion rate. In my case, should I be using a Wilcoxon test also?


title: What number minimises RMSD between itself and randomly generated numbers?

text: Let us assume that we have a random number generator given by its probability density function. Now, we need to use this distribution to generate a number that has the following property. It minimises the root mean square deviation between itself and random numbers generated by the considered distribution. My question is this number has a special name. May be it is just a mean or median?


title: When to include a random effect in a model

text: I'm new to mixed modelling and i'm confused as to whether its appropriate to use a random effect in an analysis I'm doing. Any advice would be appreciated. my study is testing how well a newly developed index of mammal abundance can predict the value of an established but more labour intensive index. i've been measuring these indices in multiple forest patches, with multiple plots in each forest patch. because i'm not directly interested in the effect of forest patches, and because my sample plots are nested within forest patches, ive been using forest patch as a random effect. However, I've got a couple of questions about this: first, i know that random effects allow you to generalise your results across all possible levels of the random factor, not just the ones you sampled. but it seems to me that to make this kind of inference your levels would have to be randomly sampled? My forest patches were not randomly sampled, so can I still use them as a random effect? second, Ive read that you can test whether it is necessary to have a random effect by doing eg a likelihood ratio test to compare models with and without the effect. I've done this, and it suggests that the random effect model does not explain the data as well as a fixed effects only model. my issue with this is that my plots are still nested within forest patches, and so presumably not independent. so, can i use this LRT approach to justify excluding the random effect, or do i still need to include it to account for nestedness? and if i do end up removing the random effect, is there a way to verify that plots within forest patches can be considered independent? Thanks for your help! Jay


title: Hypothesis test on a series of wagers with different probabilities and returns

text: I have a model that predicts the probability of an event occurring, and then based on a given return if the event occurs, calculates the expected value. If the expected value is positive, a wager is placed. If the event occurs, the outcome is the return, and if it does not occur the wager is lost. I have tested the model on historical data, and the results are positive, but not so positive that I feel confident that my average outcome is statistically different from 0. How can I conduct a hypothesis test on my results from testing on the historical data? Below is the distribution of my results, with -1 being a lost wager, 0 being an event where no wager was placed due to a negative expected value, and anything greater than 0 being the return from a successful wager. ![All Outcomes](http://i.stack.imgur.com/ordMA.jpg) ![Successful Wagers](http://i.stack.imgur.com/oNbXG.jpg)


title: Whitening data and testing for correlation

text: Before I start I would like to mention that I do NOT have a mathematical background, so please answer in a ... easy to follow manner. I'm testing 2 sets of stock market data (Shanghai Stock exchange (SSE Composite index) and Hong Kong stock exchange (HANG SENG INDEX)). I want to test for correlation between the two sets of data (I don't care about lead/lag structures). I read some of the posts on this website and from what I understand I first need to "whiten" my data before I can test for correlation. My questions are as following;   1. Do I need to whiten my data? If so, how would I do this? Is this possible in Excel, or do I need special software (which software programs could do this for me?)    2. When my data are whitened, can I just use a Pearson correlation method to test for correlation?


title: Issues in least square estimation

text: I want to perform system identification of an AR (2) system with a white gaussian signal using the LS estimation technique. Assuming the received signal is generated by channel model given by               x(t)= 0.919*x(t-1) - 0.92*x(t-2)+n(t)      where n(t) is white noise. I need to find the parameters by minimizing the least square errors between predicted and experimental values when the receiver section only gets x(t). I also need to find the 95% confidence interval for each parameter. Now MATLAB has its own system identification command `ar()` which gives us the coefficients. It is hard to follow the flow of this function and hence wanted to know how exactly the parameter estimation is performed. **EDIT this is how I did but the coefficients are incorrect.**               x(1)=0.001;     x(2)=0.002;     for i =3 : 200            x(i) =  0.195 *x(i-1) -0.95*x(i-2) ;  %output     end          d= awgn(x,10, 'measured');  %output received signal          xb = toeplitz(d,[d(1),zeros(1,N-1)]) ;% Input matrix               P1 = pinv(xb)*d;      AND when I used `model1 = ar(x, 2, 'ls');` I get the coefficients as `-0.1950 0.9500` But the signs are reversed! And when I used the noise corrupted signal, I got completely different result. Should I not be using the signal d instead of x? What is wrong? Can somebody please provide in detail what is the correct method to find using pinv() and ar() go about it?


title: Questions regarding predict.glmpath()

text: I'm trying to do LASSO in R with the package glmpath. However, I'm not sure if I am using the accompanying prediction function _predict.glmpath()_ correctly. Suppose I fit some regularized binomial regression model like so:               fit <- glmpath(x = data$x, y=data$y, family=binomial)      Then I can use predict.glmpath() to estimate the value of the response variable $y$ at $x$ for varying values of $\lambda$ through               pred <- predict.glmpath(fit, newx = x, mode="lambda", s=seq(0,10,1),type="response")      However, in the help file it can be seen that there is also an option _newy_. How should one interpret the result when calling _predict.glmpath()_ with _newy = some.y_? **[Edit]** An additional question came to mind: The option _type_ can have the following values, according to the help file:                                     description in help file          "response"            the estimated responses are returned     "loglik"              the log-likelihoods are returned     "coefficients"        the coefficients are returned. The coefficients for the initial input variables are returned (rather than the standardized coefficients)     "link"(default)       the linear predictors are returned      However, to which linear predictors and coefficients are they referring to? Surely not those of the original model?


title: Is multicollinearity a concern in nonparametric statistics?

text: Should I be concerned about multicollinearity in nonparametric statistics?


title: Identifying and quantifying variables (predictors) to understand the relationship

text: What type of analysis do I need to understand and/or gain information from a set of data? For example, I have students' data for SATScore, HighSchoolGPA, HighSchoolRank, etc., and FreshmenGPA. Now I want to develop a model which tells me which of the criteria are important for FreshmenGPA and by how much. I want to quantify the importance of these categories (SATScore, HighSchoolGPA, HighSchoolRank, etc.). What options do I have? Which is the most accurate methodology to do such analysis? What I Want to find out from the data: I want is to identify which of these variables(SATScore, HighSchoolGPA, HighSchoolRank, etc.) can help me predict Students FreshmenGPA. I did some reading on what I need to do for this but there is no definite source or a definite answer. From what I roughly understand, it looks like I need to find out which of these variables are helpful in predicting FreshmenGPA, and whether there is an interaction variable(combination of two or more variables) which is a better predictor than any individual variable. How can I do this? Once I know which variables and/or interaction variables are important for my model, I will need to quantify the variable importance, for this I'll need to find out the weights(coefficients) for variables which I found are important for FreshmenGPA. Could anyone show me the steps I need to follow to understand and do this? Also, It looks like there are multiple ways of doing this, Not sure if this is correct, but If so which of these methods give the most accurate model for the data?


title: Average of several categorical random variables

text: Say I have a categorical random variable $X$ over a discrete label space $L$ = {Sky, Road, Tree, Unknown}. In other words each $X_i \in L$. Now I store the parameters of each $X_i$ which are the individual probabilities i.e P(X_i = Sky), P(X_i = Road), ...$. Now say say I have several of these $X_1, X_2,...,X_m$. I want to get a new random variable $X_A$ which denotes the the average of $X_1, X_2,...,X_m$. So currently I am computing the parameters of $X_A$ by simply averaging the parameters of $X_1, X_2,...,X_m$. **My question is that if the above approach is correct?** * * * More Details: Since I store the parameters of each $X_i$ : P(X_i = Sky), P(X_i = Road),.. as a vector $P_i$ of size $|L|$. I am using the following equation to compute the average: $P_A = (\sum_{i=1}^{m} P_i )/m$ and then normalize as $P_A = P_A / (\sum_{l=1}^{|L|} P_A[l] )$ And if I am storing Log Probabilities in each $P_i= [\log P(X_i= Sky), \log P(X_i= Road), ..]$ then I am doing: $P_A = \log \sum_{i=1}^{m} \exp (P_i) - \log(m)$ and normalize as $P_A = P_A - \log \sum_{l=1}^{|L|} \exp(P_A[l])$ where $\log \sum \exp$ is as described here. **One side question** is that if there is more efficient way of computing this when storing log-probabilities (2nd equation)? A grounded application of this averaging can be as following: Say you 3D voxels of some size where each voxel stores a random variable $X_i$ denoting its label Sky, Road, Tree, Unknown. Now you want to move to more coarse representation with bigger voxels (octree) so you need compute the average $X_A$ of all the smaller voxels inside the new voxel.


title: Decision Level Fusion of SVR outputs

text: I have two sets of features predicting the same outputs. But instead of training everything at once, I would like to train them separately and fuse the decisions. In SVM classification, we can take the probability values for the classes which can be used to train another SVM. But in SVR, how can we do this? Any ideas? Thanks :)


title: Residual plot for regression tree: What should it look like?

text: * I realize that decision trees are nonparametric methods   * What should residual vs. actual/fitted look like for a well behaved regression tree?   * My argument would be that since each observation assigned to a terminal node is assigned (as a predicted value) the average of the dependent variable at that terminal node, you would expect the conditional distribution (that is, for each node) to be approximately normal.   * I have attached two plots for my decision tree (validates at 63% on test set, so kind of weak), residuals vs. fitted and residuals vs. actual -Basically, my question: wouldn't a strong regression tree look like a step-function of sorts? ![Residuals vs. actual](http://i.stack.imgur.com/YOF5c.jpg) ![Resdiduals vs predicted](http://i.stack.imgur.com/KYEXe.jpg)


title: Which equation will give me meaningful insight?

text: I have 10 iPads. I am logging the number of times an app is crashing each day for each of these iPads. The number of crashes tends to be skewed towards just a couple of the devices such that taking an average will not give anything meaningful. What type of equation should I use to get some meaningful data related to my app crashing? And / or (more importantly): What should I search for in order to learn about the type of equation(s) that will help me? Example dataset:               {iPad,numCrashes}     [{0,11},{1,7},{2,10},{3,0},{4,0},{5,0},{6,0},{7,0},{8,0},{9,0}]      Crashes per device is 2.8, but only 30% of the devices actually crash.


title: Estimate non-parametric prediction interval

text: Maybe let me first describe the real world application of my problem before going to a more mathematical model. Suppose I have a system that measures my current position up to an unknown error (let's say this would be a GPS-device). At every point that is recorded I also save the accurate position obtained by other means. Now I want to estimate an prediction interval of the recorded error, thus I want a reliable value that say 95 % of all data points are closer than the estimated value to the real position. Which means I basically would like to estimate the probability distribution function of the values $$e_i:=|X_i-P_i|,$$ where $X_i$ is my recorded (random) position and $P_i$ the known reference value. Thus far everything is good. Problems do now arise looking at the sample values of the $e_i$. I see highly correlated data points and a distribution function that is far away from all parametric distributions I know of. Do you know of a robust, non-parametric method of estimating the distribution (or at least one predefined quantile of that distribution) with correlated data? Sample sizes shall not be a problem, I have about 50000 data points available. Of course one way would be to compute the empiric distribution function, but does this still work with correlated data? The limit theorem I know of assumes independent samples to guarantee correct asymptotical behavior. Thanks for your input.


title: Clustering high dimensional data (p > n) in R

text: I have a situation where we have a number of quantitative features / variables (p) than the number of samples (n). My objective is to classify these samples into groups (may be hierarchical). I can see a good discussion about this in this Q/A post, here at CV. I am aware of discussion about clustering based on high-dimensional data in wiki and its needs. Here is data example for workout:               set.seed(123)     # matrix of X variable      xmat <- matrix(sample(-1:1, 2000000, replace = TRUE), ncol = 10000)     colnames(xmat) <- paste ("M", 1:10000, sep ="")     rownames(xmat) <- paste("sample", 1:200, sep = "")      Here are my questions:   1. What would be best approach ?    2. I am interested to find implementation codes for a suitable method (may be Subspace clustering or Projected clustering or Correlation clustering or Hybrid approaches) for my case.


title: Iterative Addition of Variables to Model Based on P Value

text: Suppose I have 64 columns that I have chosen out of 500+ columns based on the fact that they have the highest pairwise correlation (is this a good way?). I take 16 of these columns and run a simple multi-variate regression. I find that 4 of these columns have P values less than .05. Can I go back, do this for the next set of 16 (mutually exclusive from the first set), check the P values, and add the successful columns to the model without issue? Moreover, is it okay to only pick the four mutually exclusive sets of 16 columns, or would I have to actually pick 64 choose 16 times and find commonalities between each set? I have a feeling there is something wrong with the way I am doing the above, so I hope someone can shed some light on the proper way to go about this. As well, please let me know if there is anything I can do to improve this question.


title: Appropriate Statistical Test for multiple groups

text: I have a question regarding which test I should run on some data. I have small groups of people who answered questions. The answers to these questions have been classified as correct or incorrect. The N of the total number of questions asked is fairly large, however, the N of the total number of groups is quite small. I'd like to analyse the accuracy based on group size. My initial thought was an ANOVA would be appropriate as I would divide the groups into size 1-2, 3-4, and 5+ and use overall group accuracy for each instance. However, as I said, the N is small. A colleague suggested that this is incorrect and a Chi Squared would be more appropriate where the N in question would be the instances of all the correct/incorrect responses instead of the overall group accuracy. Is either one of these or another test more appropriate?


title: lmmlasso does not work for p>n

text: I am using the following example from the R manual:               # (2) Use the lmmlasso on a small simulated data set     set.seed(54)     N    <- 20 # number of groups     p    <- 6 # number of covariates (including intercept)     q    <- 2 # number of random effect covariates     ni   <- rep(6,N) # observations per group     ntot <- sum(ni) # total number of observations     grp  <- factor(rep(1:N,each=ni)) # grouping variable     beta <- c(1,2,4,3,0,0,0) # fixed-effects coefficients     x    <- cbind(1,matrix(rnorm(ntot*p),nrow=ntot)) # design matrix     bi1  <- rep(rnorm(N,0,3),each=ni) # Psi=diag(3,2)     bi2  <- rep(rnorm(N,0,2),each=ni)     bi   <- rbind(bi1,bi2)     z    <- x[,1:2,drop=FALSE]     y    <- numeric(ntot)     for(k in 1:ntot) y[k] <- x[k,]%*%beta + t(z[k,])%*%bi[,grp[k]] + rnorm(1)     # correct random effects structure     fit2 <- lmmlasso(x=x,y=y,z=z,grp=grp,lambda=10,pdMat="pdDiag")     summary(fit2)      When I am trying to increase the value of `p` to be more than `ntot`, then I get an error of `Redundant covariance parameters` (after I have increased the lambda considerably). I am not sure if I am missing something here. Any help will be highly appreciated. My code is the following:               # lmmlasso used on a simulated dataset     set.seed(54)     require(lmmlasso)     N    <- 30           # number of groups     p    <- 500          # number of covariates (including intercept)     q    <- 1            # number of random effect covariates     ni   <- rep(6,N)     # observations per group     ntot <- sum(ni)      # total number of observations          grp  <- factor(rep(1:N,each=ni)) # grouping variable          beta <- c(1,2,4,3,3,rep(0,496)) # fixed-effects true coefficients     x    <- cbind(1,matrix(rnorm(ntot*(p)),nrow=ntot)) # X design matrix     z    <- x[,1,drop=FALSE]#Z design matri is a subset of X design matrix          bi1  <- (rep(rnorm(N,0,sqrt(0.56)),each=ni))      bi   = rbind(bi1)          y    <- numeric(ntot)     for(k in 1:ntot) y[k] <- x[k,]%*%beta + t(z[k,])%*%bi[,grp[k]] + rnorm(1,0,sqrt(0.25))          # correct random effects structure     fit2 <- lmmlasso(x=x,y=y,z=z,grp=grp,lambda=1200,pdMat="pdDiag")     summary(fit2)


title: How to tell which recipe rating is better?

text: Comparing ratings for recipes I have often wondered is, say, 50 reviews of average 4.5 better than 31 reviews of average 4.8? Orbiting these considerations has been some idea of statistical significance. All uncertainty threatens to be dispelled today since I decided to ask here. No doubt there is a clear answer to the question, is it possible to determine simply from the number of ratings (and their average, or even their distribution), which recipe is better? And if so, what is a quick formula to mentally evaluate this? **Edit** : To provide some more clarity, a definition of 'better':   1. All reviews are rated on the real interval $[0,5]$. A recipe is modelled as having some real, independent, "objective" "score" (on the same interval) that is estimated, with error, by the ratings of the recipe reviewers. It is to this "essential" score we refer when we assert one recipe is better in comparison to another. We attempt to access this essential score through models of it, such as through the average ratings, and number of ratings of that recipe.   2. If the ratings (of two recipes $F$ and $G$ for the same enreciped menu item, $M$) are on the same site, and have the same number of reviews, then barring the inclusion of any additional information, recipe $F$ is better than $G$ if $avg(ratings(F)) > avg(ratings(G))$.   3. If the ratings are on the same site, and have different numbers of reviews, **I don't know.**   4. If the ratings are on different sites, then, barring the inclusion of any additional information, as per 2 and 3.   5. If additional information is available, such as, for example (and readers may choose to delve as deeply as they like in this direction):      * the reliability of site $A$ versus site $B$. This is left to the reader's discretion, and it is imagined a probability value giving the "confidence of belief in the reliability of" site $A$'s ratings, that is, a measure of site $A$'s rating's fidelity to the underlying "true" value of the recipe's rating. So, if site $A$ provides very high ratings to recipe of an item $M$, but is generally thought to be unreliable signal of betterness, while site $B$ provides somewhat lower average, with much greater reliability, to a different recipe of the same item $M$, one could make a reasonable case that site B's recipe is better than site $A$'s. Any metric of site reliability and associated mechanism of computing an ordering on recipe's must reach the same conclusion.      * the distribution of reliability of the raters, that is, a measure of how closely a rater's score for any recipe, say $score(rater, recipe)$, models the recipe's underlying "objective" score.  In terms of what's required to answer this question, 3 is necessary and sufficient, and 5 is unnecessary and (alone) insufficient, although welcome.


title: Specifying lag in `dlnm` when passing arguments to `crossbasis`

text: I am using the `dlnm` package to build a finite distribute lag linear model. I intend on testing the model-fit based on various lag levels to assess which lag is suitable. Needless to mention I will apply some domain knowledge to make a good call. I am using these two available resources to carry out this rather complex exercise: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3191524/ http://cran.r-project.org/web/packages/dlnm/vignettes/dlnmOverview.pdf I first created a Matrix to include non-linear effects of both my predictors using this: `ImpressionsA.x<-onebasis(locanatmodelset$ImpressionsA.x, fun="poly", degree=2)` `ImpressionsA.y<-onebasis(locanatmodelset$ImpressionsA.y, fun="poly", degree=2)` Since my non-linear relationship is tending to be quadratic, I am using 2nd degree polynomials. Now I need to factor in the lags for each variable by calling `crossbasis`. This is where I am confused. The function as per the package documentation is of the format: `cb <- crossbasis(chicagoNMMAPS$temp,lag=30,argvar=list("thr",thr.value=c(10,20)), arglag=list(knots=c(1,4,12)))` my doubts are : 1) in `lag=30` is 30 of the same units as my response variable or will it be in the units of my predictor variable? In my case, my lag is specified in days. I want to specify 5 days as lag before fitting the model. How should I pass the argument? 2) Since I have already created my basis matrix for predictor variables using `onebasis`, how should I pass the arguments `argvar` and `arglag`? 3) I also wanted to extract the lagged values of my predictor variables (ImpressionsA.x and ImpressionsA.y) in their original units. The matrix itself is not helpful. It transforms everything into negative values on some other scale all together.


title: Difference powerlaw, lognormal and stretched exponential (Weibull) function

text: I am currently fitting above mentioned functions to my data and I can observe, that both lognormal and Weibull are better fits than a power-law. In the literature it is often suggested that it is hard to distinguish these functions, but I can in my data. I do not completely understand the differences between them though. What might be an explanation for the underlying process that generates my data that I can infer from these results?


title: Is Cohen's $d$ a better measure of effect-size than a simple mean difference for a meta-analysis?

text: Cohen's $d$ measures the difference between the mean for the experimental group and the control group divided by standard deviation. The Cohen's $d$ statistic produces something that does not add any value to the observed effect-size i.e. the difference between the mean for the experimental group and the mean of control group.Should we use Cohen's $d$ or simple difference produced by a controlled experiment for a meta-analysis?


title: the relationship between statistics and probability

text: > **Possible Duplicate:**   >  What's the difference between probability and statistics? Today I just came back from a student seminar in our department. And I heard that for the Pareto distribution, it is kind of difficult to estimate the parameter. There are hundreds of estimations of the parameter and they differ a lot, which means there is no way to estimate it. Thus, I am just wondering what is statistics? Statistics could not estimate such a parameter? Statistics relies what to get the estimation? Not Probability properties of this distribution? Statistics just relies on a mass data set and how could it come up a statistic to estimate the parameter of interest? We use mean to estimate the expectation just because of the central limit theorem, which is totally probability. So what is statistics?


title: Generate variables with fixed Rsquare and coefficients

text: I have a random variable `X`.               set.seed(1)     X <- rnorm(10,0,30)      I would like to generate three other variables `(y1, y2, y3)` that have a pre- specified R^2 of `.50`, `.36` and `.26`. i.e: `cor(x,cbind(y1,y2,y3))^2=0.5,0.36,0.26` and also the coefficients for a linear model where y1 y2 and y3 are regressed on X return: `0.4, 0.25 and 0.1`, i.e.: `lm(X~y1+y2+y3) -> 0.4*y1 + 0.25*y2 + 0.15*y3` Is this possible to do starting from X?


title: Automating statistical correlation between "texts" and "data"

text: I am collecting textual data surrounding press releases, blog posts, reviews, etc of certain companies' products and performance. Specifically, I am looking to see if there are correlations between certain _types_ and/or _sources_ of such "textual" content with market valuations of the companies' stock symbols. Such _apparent_ correlations can be found by the human mind fairly quickly - but that is not scalable. How can I go about automating such analysis of disparate sources?


title: Is it typical for group means of composite z-scores to have means with values of the same magnitude but opposite signs?

text: I developed a composite score by converting four items to z-scores, then I summed up the z-scores to form a composite measure. However, when I've calculated means to compare finding for two different groups on a couple datasets using this composite score, the mean scores are the same distance from zero. For example,   * For one dataset, M1= 0.76 and M2= -0.76    * For anotherdataset M1= 1.34 and M2= -1.34.  Is this typical finding for this type of composite score? It seems odd to me.


title: Interpolation of influenza data that conserves weekly mean

text: **Edit** I have found a paper describing exactly the procedure I need. The only difference is that the paper interpolates monthly mean data to daily, while preserving the monthly means. I have trouble to implement the approach in `R`. Any hints are appreciated. **Original** For each week, I have the following count data (one value per week):   * Number of doctors' consultations   * Number of cases of influenza My goal is to obtain _daily_ data by interpolation (I thought of linear or truncated splines). The important thing is that I want to **conserve the weekly mean,** i.e. the mean of the daily interpolated data should equal the recorded value of this week. In addition, the interpolation should be smooth. One problem that could arise is that a certain week has less than 7 days (e.g. at the beginning or end of a year). I would be grateful for advice on this matter. Thanks a lot. Here's a sample data set for the year 1995 ( **updated** ):               structure(list(daily.ts = structure(c(9131, 9132, 9133, 9134,      9135, 9136, 9137, 9138, 9139, 9140, 9141, 9142, 9143, 9144, 9145,      9146, 9147, 9148, 9149, 9150, 9151, 9152, 9153, 9154, 9155, 9156,      9157, 9158, 9159, 9160, 9161, 9162, 9163, 9164, 9165, 9166, 9167,      9168, 9169, 9170, 9171, 9172, 9173, 9174, 9175, 9176, 9177, 9178,      9179, 9180, 9181, 9182, 9183, 9184, 9185, 9186, 9187, 9188, 9189,      9190, 9191, 9192, 9193, 9194, 9195, 9196, 9197, 9198, 9199, 9200,      9201, 9202, 9203, 9204, 9205, 9206, 9207, 9208, 9209, 9210, 9211,      9212, 9213, 9214, 9215, 9216, 9217, 9218, 9219, 9220, 9221, 9222,      9223, 9224, 9225, 9226, 9227, 9228, 9229, 9230, 9231, 9232, 9233,      9234, 9235, 9236, 9237, 9238, 9239, 9240, 9241, 9242, 9243, 9244,      9245, 9246, 9247, 9248, 9249, 9250, 9251, 9252, 9253, 9254, 9255,      9256, 9257, 9258, 9259, 9260, 9261, 9262, 9263, 9264, 9265, 9266,      9267, 9268, 9269, 9270, 9271, 9272, 9273, 9274, 9275, 9276, 9277,      9278, 9279, 9280, 9281, 9282, 9283, 9284, 9285, 9286, 9287, 9288,      9289, 9290, 9291, 9292, 9293, 9294, 9295, 9296, 9297, 9298, 9299,      9300, 9301, 9302, 9303, 9304, 9305, 9306, 9307, 9308, 9309, 9310,      9311, 9312, 9313, 9314, 9315, 9316, 9317, 9318, 9319, 9320, 9321,      9322, 9323, 9324, 9325, 9326, 9327, 9328, 9329, 9330, 9331, 9332,      9333, 9334, 9335, 9336, 9337, 9338, 9339, 9340, 9341, 9342, 9343,      9344, 9345, 9346, 9347, 9348, 9349, 9350, 9351, 9352, 9353, 9354,      9355, 9356, 9357, 9358, 9359, 9360, 9361, 9362, 9363, 9364, 9365,      9366, 9367, 9368, 9369, 9370, 9371, 9372, 9373, 9374, 9375, 9376,      9377, 9378, 9379, 9380, 9381, 9382, 9383, 9384, 9385, 9386, 9387,      9388, 9389, 9390, 9391, 9392, 9393, 9394, 9395, 9396, 9397, 9398,      9399, 9400, 9401, 9402, 9403, 9404, 9405, 9406, 9407, 9408, 9409,      9410, 9411, 9412, 9413, 9414, 9415, 9416, 9417, 9418, 9419, 9420,      9421, 9422, 9423, 9424, 9425, 9426, 9427, 9428, 9429, 9430, 9431,      9432, 9433, 9434, 9435, 9436, 9437, 9438, 9439, 9440, 9441, 9442,      9443, 9444, 9445, 9446, 9447, 9448, 9449, 9450, 9451, 9452, 9453,      9454, 9455, 9456, 9457, 9458, 9459, 9460, 9461, 9462, 9463, 9464,      9465, 9466, 9467, 9468, 9469, 9470, 9471, 9472, 9473, 9474, 9475,      9476, 9477, 9478, 9479, 9480, 9481, 9482, 9483, 9484, 9485, 9486,      9487, 9488, 9489, 9490, 9491, 9492, 9493, 9494, 9495), class = "Date"),          wdayno = c(0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L,          5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L,          6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L,          0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L,          1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L,          2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L,          3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L,          4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L,          5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L,          6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L,          0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L,          1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L,          2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L,          3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L,          4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L,          5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L,          6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L,          0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L,          1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L,          2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L,          3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L,          4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L,          5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L,          6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L,          0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L), month = c(1, 1, 1, 1, 1,          1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,          1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,          2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3,          3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,          3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,          4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,          4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,          5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6,          6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,          6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,          7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8,          8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,          8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9,          9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,          9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,          10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,          10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,          11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,          11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,          12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,          12, 12, 12, 12), year = c(1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,          1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995), yearday = 0:364,          no.influ.cases = c(NA, NA, NA, 168L, NA, NA, NA, NA, NA,          NA, 199L, NA, NA, NA, NA, NA, NA, 214L, NA, NA, NA, NA, NA,          NA, 230L, NA, NA, NA, NA, NA, NA, 267L, NA, NA, NA, NA, NA,          NA, 373L, NA, NA, NA, NA, NA, NA, 387L, NA, NA, NA, NA, NA,          NA, 443L, NA, NA, NA, NA, NA, NA, 579L, NA, NA, NA, NA, NA,          NA, 821L, NA, NA, NA, NA, NA, NA, 1229L, NA, NA, NA, NA,          NA, NA, 1014L, NA, NA, NA, NA, NA, NA, 831L, NA, NA, NA,          NA, NA, NA, 648L, NA, NA, NA, NA, NA, NA, 257L, NA, NA, NA,          NA, NA, NA, 203L, NA, NA, NA, NA, NA, NA, 137L, NA, NA, NA,          NA, NA, NA, 78L, NA, NA, NA, NA, NA, NA, 82L, NA, NA, NA,          NA, NA, NA, 69L, NA, NA, NA, NA, NA, NA, 45L, NA, NA, NA,          NA, NA, NA, 51L, NA, NA, NA, NA, NA, NA, 45L, NA, NA, NA,          NA, NA, NA, 63L, NA, NA, NA, NA, NA, NA, 55L, NA, NA, NA,          NA, NA, NA, 54L, NA, NA, NA, NA, NA, NA, 52L, NA, NA, NA,          NA, NA, NA, 27L, NA, NA, NA, NA, NA, NA, 24L, NA, NA, NA,          NA, NA, NA, 12L, NA, NA, NA, NA, NA, NA, 10L, NA, NA, NA,          NA, NA, NA, 22L, NA, NA, NA, NA, NA, NA, 42L, NA, NA, NA,          NA, NA, NA, 32L, NA, NA, NA, NA, NA, NA, 52L, NA, NA, NA,          NA, NA, NA, 82L, NA, NA, NA, NA, NA, NA, 95L, NA, NA, NA,          NA, NA, NA, 91L, NA, NA, NA, NA, NA, NA, 104L, NA, NA, NA,          NA, NA, NA, 143L, NA, NA, NA, NA, NA, NA, 114L, NA, NA, NA,          NA, NA, NA, 100L, NA, NA, NA, NA, NA, NA, 83L, NA, NA, NA,          NA, NA, NA, 113L, NA, NA, NA, NA, NA, NA, 145L, NA, NA, NA,          NA, NA, NA, 175L, NA, NA, NA, NA, NA, NA, 222L, NA, NA, NA,          NA, NA, NA, 258L, NA, NA, NA, NA, NA, NA, 384L, NA, NA, NA,          NA, NA, NA, 755L, NA, NA, NA, NA, NA, NA, 976L, NA, NA, NA,          NA, NA, NA, 879L, NA, NA, NA, NA), no.consultations = c(NA,          NA, NA, 15093L, NA, NA, NA, NA, NA, NA, 20336L, NA, NA, NA,          NA, NA, NA, 20777L, NA, NA, NA, NA, NA, NA, 21108L, NA, NA,          NA, NA, NA, NA, 20967L, NA, NA, NA, NA, NA, NA, 20753L, NA,          NA, NA, NA, NA, NA, 18782L, NA, NA, NA, NA, NA, NA, 19778L,          NA, NA, NA, NA, NA, NA, 19223L, NA, NA, NA, NA, NA, NA, 21188L,          NA, NA, NA, NA, NA, NA, 22172L, NA, NA, NA, NA, NA, NA, 21965L,          NA, NA, NA, NA, NA, NA, 21768L, NA, NA, NA, NA, NA, NA, 21277L,          NA, NA, NA, NA, NA, NA, 16383L, NA, NA, NA, NA, NA, NA, 15337L,          NA, NA, NA, NA, NA, NA, 19179L, NA, NA, NA, NA, NA, NA, 18705L,          NA, NA, NA, NA, NA, NA, 19623L, NA, NA, NA, NA, NA, NA, 19363L,          NA, NA, NA, NA, NA, NA, 16257L, NA, NA, NA, NA, NA, NA, 19219L,          NA, NA, NA, NA, NA, NA, 17048L, NA, NA, NA, NA, NA, NA, 19231L,          NA, NA, NA, NA, NA, NA, 20023L, NA, NA, NA, NA, NA, NA, 19331L,          NA, NA, NA, NA, NA, NA, 18995L, NA, NA, NA, NA, NA, NA, 16571L,          NA, NA, NA, NA, NA, NA, 15010L, NA, NA, NA, NA, NA, NA, 13714L,          NA, NA, NA, NA, NA, NA, 10451L, NA, NA, NA, NA, NA, NA, 14216L,          NA, NA, NA, NA, NA, NA, 16800L, NA, NA, NA, NA, NA, NA, 18305L,          NA, NA, NA, NA, NA, NA, 18911L, NA, NA, NA, NA, NA, NA, 17812L,          NA, NA, NA, NA, NA, NA, 18665L, NA, NA, NA, NA, NA, NA, 18977L,          NA, NA, NA, NA, NA, NA, 19512L, NA, NA, NA, NA, NA, NA, 17424L,          NA, NA, NA, NA, NA, NA, 14464L, NA, NA, NA, NA, NA, NA, 16383L,          NA, NA, NA, NA, NA, NA, 19916L, NA, NA, NA, NA, NA, NA, 18255L,          NA, NA, NA, NA, NA, NA, 20113L, NA, NA, NA, NA, NA, NA, 20084L,          NA, NA, NA, NA, NA, NA, 20196L, NA, NA, NA, NA, NA, NA, 20184L,          NA, NA, NA, NA, NA, NA, 20261L, NA, NA, NA, NA, NA, NA, 22246L,          NA, NA, NA, NA, NA, NA, 23030L, NA, NA, NA, NA, NA, NA, 10487L,          NA, NA, NA, NA)), .Names = c("daily.ts", "wdayno", "month",      "year", "yearday", "no.influ.cases", "no.consultations"), row.names = c(NA,      -365L), class = "data.frame")


title: Testing of sensitivity for different algorithms

text: I have some queries I would like to discuss, I am doing some analysis for sensitivity. I am not from the field of statistics. I work with RNA-Seq data. The thing is I have found some genes which are differentially expressed by 4 different methods at different replicates. Out of which one of the method is our own developed algorithm. I want to now do the test of positives with my data. So I take all the genes given by each method and take the union out of it. So this gives me the universe of all the genes from all the method. Now I want to see at different pair replicate models how many positives are there for each of the methods. I want to do a curve where I can show at each replicate how number of positives are there for each method and alternatively the number of false positives at different replicate from each method. Can anyone tell me if my explanation is correct or not. For the positives I can simply take in account the percentage shared at each replicate for different methods calculating from the union set of all methods and then plot the curve and then for the false positives what should I do? It should be the remaining percentage right that are not shared in the union? I would like to have some understanding here.


title: Multiple regression or multivariate regression

text: Are there any difference in beta coefficients when doing several multiple regression analysis as compared to doing multivariate regression?


title: Interpretation of $\theta$ in negative binomial regression

text: First off, a very similar question has been asked before. But the answers to this question did not explain what high/low values of theta mean. Here's my crack at trying to figure out what high/low values of theta mean. So please don't close this question! Let's assume you've made two models: a negative binomial regression (NB) and a zero-inflated negative binomial regression (ZINB). The NB regression has a theta of 0.5 and the ZINB regression has a theta of 2. As I understand it, the higher theta in the ZINB regression indicates that more variance in the residuals has been accounted for, and therefore the negative binomial distribution that the model assumes has a more slender shape. Is this correct? Can anybody provide a more precise definition of the theta value, but without using equations? I also quickly sketched a visualisation of my understanding. The residuals in the NB are more spread out, meaning the theta is smaller and the shape of the negative binomial distributions are more fat. The residuals in the ZINB are less spread out, meaning the theta is larger and the shape of the negative binomial distributions are more slender. ![enter image description here](http://i.stack.imgur.com/Fc5ez.jpg)


title: How to calculate the variance of vectors for clustering?

text: I'm interested in various methods of measuring dispersion of vectors mainly for use in cluster analysis. I can think of three methods:   1. Find the mean vector (centroid), then calculate the variance of the distances of all vectors to this mean vector. It is possible that the set of vectors could all be different, but have the same distance to the mean vector. In this case this would not appear to be a great measure, though this situation may be unlikely in practice. It seems that the Davies-Bouldin cluster quality measure uses this to measure intra-cluster quality.    2. Use the mean pairwise distance between vectors. I've seen this used to measure both intra and inter cluster quality. This would seem to require some sort of distance matrix. The implementation may be difficult if one tries to add or remove vectors and update on the distance matrix on the fly.    3. Calculate the population variance for each component of the vectors. This would result in a vector containing the population variance for each component. Then take the sum of the components in this vector.  My questions:   * Any thoughts on these measures?    * Any other good measures?    * Also does anyone know a one pass algorithm for calculating #1 and #2.  I know how to compute #3 with a numerically stable one pass algorithm. Essentially every time I add or remove a vector from a cluster I would like the measure of quality to update automatically. I've had some luck with this for certain measures.


title: Use of data from ROC curve

text: In order to find an optimal time for initiation of treatment post surgery (oncologic patients) I created a ROC curve with death defined as event. The AUC was not significant. However, I decided to use the time that provided the best specificity and sensitivity in cox regression as a binary variable. The Cox parameter was significant. I was wondering whether I am committing a sin using data from ROC curve with p-value > 0.05 and by using the data obtained from exploratory analysis to analyse the same training data. I would be very interested in hearing feedback and suggestions (maybe for validation?) Thanks in advance, Liran


title: Likelihood based model selection

text: Let's say I got a set of models $M = \\{M_1, M_2, \dots M_n\\}$. Now say I got some data $x$ and I would like to know, which model represents the data best. I know how to calculate the likelihood $L(\theta | x)$, with $\theta$ being the parameters of any of those models. I realize that the likelihood value of one model alone won't tell me anything useful. But what I can do is compare them to each other. Now I know which of the given models is the most likely. **But** : I would also like to know, how likely it is none of the models represents a model well enough? That is, I'm interested in a statistical sound way to tell, that I should create a new model for that data. Any pointers on how I could calculate this?


title: Want to evaluate how uniform my data is?

text: I have a feature vector of size nx3 where n = no. of data points and no. of features = 3.   So the mean is a 3x1 vector. The data actually is a set of pixels of an image segment(not a whole image). and the 3 features are R,G,B for the time being(, although it can be anything). My goal: To define a measure that will evaluate how uniform or homogeneous is my data (or image segment). I want to calculate this score for all such feature vectors (or all such image segments). Then, I will apply a threshold to decide if each of my feature vector is uniform or non-uniform. P.S: I do not want to do any kind of clustering among my data points. My concern : When I find variance (`var` command in matlab), I get a 3x1 vector, each element for each feature.   I instead want to find one scalar value, as in, the average euclidean distance of all data points from the mean point. something like this: `(1/n) * sum(sqrt( (x1^2 - m1^2) + (x2^2 - m2^2) + (x3^2 - m3^2)))` `[x1 x2 x3]` is ith sample. `[m1 m2 m3]` is the mean vector. (1) Is this right way to measure the uniformity? (2) What is it called? should I call is mean squared error(MSE)? (3) Is there any better way I should evaluate the uniformity of my data points? Thanks :)


title: Reporting operative effect in paired t test

text: I just want to make sure I have something clear in my head. When I calculate the effect size for a paired samples t-test after obtaining a significant result, I simply take the mean of the differences divided by the standard deviation of the differences to get an effect size d. Do I then need to take d and divide it by the square root of 1-r, where r is the correlation between pairs and r is estimated from the sample pairs? I am confused because dividing by the square root of 1-r supposedly gives me an "operative effect size" and I'm not really sure if the operative effect size is what I should be reporting in my analyses. For example, in this report I am working on, I need to know if there was an effect size of 2 SD. So when I calculate my effect size, should I be dividing by square root 1-r? I don't think so, I think I need to report the actual detected effect size and not the operative effect, but I would love a second opinion. Thanks!


title: F-test for equality of regression coefficients

text: I have been trying to look for a reference on the theory behind using the F-test to test for the equality of regression coefficients. The closest I could find is this, which is for general linear restrictions: http://www.mattblackwell.org/files/teaching/ftests.pdf. Google seems to have failed me. Could you point me in the right direction for a theoretical reference on using F-tests to test for equality of regression coefficients? Thanks a lot in advance!


title: Best subset selection

text: My statistical learning text claims that for best subset selection, 2^p total models must be fit through regression if for p covariates, we fit p choose k models at each k, k = 1,...,p. I interpret this mathematically as 2^p = p choose 1 + ... + p choose k. Why 2^p models?


title: Testing survival against frequency of some event

text: I have the following data for 200 cases/subjects:   * time to death over a 15 year period, $t$. A datum with a value of 180 months means that the subject did not die over the 15 year period.   * the frequency of three particular types of event associated with each subject $a_{i}$, $b_{i}$ and $c_{i}$, Apart from the obvious survival statistics and KM curves, how can I test the association between $t$ and each of $a_{i}$, $b_{i}$ and $c_{i}$. I am hypothesizing that a profile of $t$ will be distinct for each of $a_{i}$, $b_{i}$ and $c_{i}$. One way I was thinking of testing this hypothesis was to simply do Pearson's correlation between $t$ and each of $a_{i}$, $b_{i}$ and $c_{i}$. I could then visualize this with a scatter plot of $t$ _vs_ each of $a_{i}$, $b_{i}$ and $c_{i}$. Any suggestions on a more sophisticated analysis including variations of KM analysis? Code examples with R appreciated.


title: Multivariate Data into VAR model in the vars package

text: I have been trying to use the vars package in r for my multivariate time series data. I am having some problems. Do I have to do something to my multivariate data before using the vars package? Can someone help me with examples. Thank you. Sam


title: How do I conduct a $\chi^2$ test on two different surveys?

text: I conducted a pre and post survey that consisted of five questions with answers being either agree or disagree. The survey was to assess if a persons perception of Pit Bulls was improved after they interacted with it. The first or pre survey was before the individual interacted with the dog. After they did, the respondents immediately were given the same survey (the post survey). I now have the results of the two surveys and am uncertain how I am to conduct the $\chi^2$ test to see if there was a significant improvement in the perception of the dog? Specifically I do not know if I am to assess each question individually or if the numbers are compiled for all questions.


title: Getting fixed-effect only predictions from mixed model on new data in R

text: I would like to construct predictions for a mixed model (logistic via glmer) on a new data set using only the fixed effects, holding the random effects to 0. But I am having trouble setting up the model matrix to be able to calculate them. Since the mer class doesn't have a predict method, and since I want to omit the random effects for predictions on the new data set, I think I need to construct a model matrix for the fixed effects of the same structure used in the original model, but using the new data. Then multiply by the fixed effect coefficients in the model. The fixed effect portion of my model formula contains factors and interaction terms between numeric fixed effects, so it's a little more complicated than just extracting the fixed variables from the matrix. e.g. I need to ensure the factor contrast expansion is the same as the original, interaction terms are properly listed, etc. So my question is: what is the more straightforward general approach for constructing a new model matrix that mimics the structure of the original model matrix used in creating the model? I've tried model.matrix(my.model, data=newdata) but that seems to return the original model matrix, not one based on newdata. Sample code:               library(lme4)          cake2 <- head(cake) # cake2 is "new" data frame for future predictions          # recipe is a fixed effect factor, temp is fixed effect numeric, replicate is random effect     m <- lmer(angle ~ temp + recipe + (1 | replicate), data=cake)     summary(m)          nrow(cake2)         # but new data frame has 6 rows     nrow(cake)          # original data frame has 270 rows          # attempt to make new model matrix using different data frame     mod.mat.cake2 <- model.matrix(m, data=cake2)     nrow(mod.mat.cake2) # 270 rows, same as orig data frame      I tried other methods like extracting the terms from the formula and building a new formula from that, but it seemed overly convoluted, and brittle in handling factors and interaction terms. How can I get mod.mat.cake2 to be a fixed effect model matrix based on the formula in m, but using values from cake2? Or is there an easier way to go about getting fixed-effect only predictions from an lmer model? All help is appreciated. Thank you.


title: Overlapping time series: is there any better way to visualize them?

text: I have this time series dataset: ![enter image description here](http://i.stack.imgur.com/eXXhd.png) The graph shows trend lines for 7 stock prices. They are very close and overlapping, but you will be able to get an idea that trend lines are layered (i.e. brown on the top and red/orange at the bottom, though far from conspicuous). Is there any way to better visualize this data? like transforming $y$-axis to another scaling, mapping the whole thing onto cylinder/cone etc.? I tried with moving average, but the improvement is not so good. NOTE: This is not an ML/DM problem. I am looking for a better/alternative/suitable visualization technique, that's all.


title: Model loyality effect of returning customers

text: I have customers which book an offer about once a year. I'd like to determine if there is a "loyality effect" - i.e. a customer from this year is more likely to book next year again. I can see that after booking this year, a customer will most likely book either about one year later, or next January. Moreover, just by chance an arbitrary customer will book with me (with a general demand peak in January or Midyear). Can you suggest an approach to split and quantify these effects from booking dates I know for my customers? (therefore I know: customer1:[date1, date2, ...], customer2:[date3, ...], ...) I'd like to separate the general random demand from the loyality bonus. Any ideas?


title: R: random sampling for multivariate normal and log-normal distributions

text: I want to generate random monthly _(m)_ temperature ( _T_ ) and Precipitation ( _P_ ) data considering that both variables are intercorrelated ( _rTP[m]_ ) The tricky thing is that my random variables that have specific quantitative properties: temperatures are normally distributed, while precipitations follow a log-normal distribution and should be log-transformed **mvrnorm** of the package MASS could be used.               mT=c(1,2,4,7,10,15,17,18,17,10,5,1)     mP=c(3.9,3.7,3.9,4.1,4.5,4.7,4.8,4.8,4.4,4.1,4.2,3.9) #log-transformed     sdT=c(1,1,1,1,1,1,1,1,1,1,1,1)     sdP=c(0.7,0.8,0.7,0.6,0.4,0.4,0.4,0.5,0.6,1,0.8,0.6)  #log-transformed     rTP=c(0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4)     covTP=rTP*(sdT*sdP)          simP=NULL     for (m in 1:12)     {     out=mvrnorm(500, mu = c(mT[m],mP[m]), Sigma = matrix(c(sdT[m]*sdT[m],covTP[m],covTP[m],sdP[m]*sdP[m]), ncol = 2), empirical = TRUE)     simP[m]=mean(exp(out[,2])-1)     }      In this case I generate two random time-series that are inter-correlated which is great. However, simulated precipitations ( _simP_ ) are on average higher than the observed one ( _mP_ )               plot(exp(mP)-1, type="l", lwd=2, ylim=c(0,250)); points(simP, type="l", lwd=2, lty=2)      I could use **rlnorm** or **rlnorm.rplus** to consider than precipitations are log-transformed, but then I have troubles with temperatures that are normally distributed. My question is: How can I create random sampling for variables that have specific quantitative properties (log-normal and normal distributions)? Thanks !


title: Logistic Regression, SVM or NN?

text: Just attended Andrew Ng’s online course on ML and although I’ve understood the methods I seem to be missing the intuition on where to apply them in terms of classification problems. What are the pros and cons of each method and in which situation should I chose one over the others? From the course I got the impression that NN is a more complex algorithm and therefore capable of identification more complex, no linear, patterns but might be computationally expensive to train. What is the advantage of SVM over LR?


title: Using association rules for classification (unusual results in parameter tuning)

text: I am using association rules for a classification problem. I am trying to predict tags for short pieces of text. What I am trying to do now is tune the confidence and support parameters in my association rules so that I am only using the rules that produce the best classification performance in an out-of- sample test set. I am getting some strange results. The best confidence seems to be around 55% (not so strange) but the best support is 1. The more I lowered support, the better the classification performance! This seems to contradict some of the literature that I have read which advises to use a support that is about 10% of the test set size. My test set size is of several million. I have checked my code several times and I am pretty sure that it is correct. The performance metric that I am using is mean F1-score. Can anyone think of a reason why this is happening?


title: Clarification on interpreting confidence intervals?

text: My current understanding of the notion "confidence interval with confidence level $1 - \alpha$" is that if we tried to calculate the confidence interval many times (each time with a fresh sample), it would contain the correct parameter $1 - \alpha$ of the time. Though I realize that this is not the same as "probability that the true parameter lies in this interval", there's something I want to clarify. **[Major Update]** Before we calculate a 95% confidence interval, there is a 95% probability that the interval we calculate will cover the true parameter. After we've calculated the confidence interval and obtained a particular interval $[a,b]$, we can no longer say this. We can't even make some sort of non-frequentist argument that we're 95% sure the true parameter will lie in $[a,b]$; for if we could, it would contradict counterexamples such as this one: What, precisely, is a confidence interval? I don't want to make this a debate about the philosophy of probability; instead, I'm looking for a precise, mathematical explanation of the how and why seeing the particular interval $[a,b]$ changes (or doesn't change) the 95% probability we had before seeing that interval. If you argue that "after seeing the interval, the notion of probability no longer makes sense", then fine, let's work in an interpretation of probability in which it _does_ make sense. _More precisely:_ Suppose we program a computer to calculate a 95% confidence interval. The computer does some number crunching, calculates an interval, and refuses to show me the interval until I enter a password. Before I've entered the password and seen the interval (but after the computer has already calculated it), what's the probability that the interval will contain the true parameter? It's 95%, and _this part is not up for debate_ : this is the interpretation of probability that I'm interested in for this particular question (I realize there are major philosophical issues that I'm suppressing, and this is intentional). But as soon as I type in the password and make the computer show me the interval it calculated, the probability (that the interval contains the true parameter) could change. Any claim that this probability never changes would contradict the counterexample above. In this counterexample, the probability could change from 50% to 100%, but...   * Are there any examples where the probability changes to something other than 100% or 0% (EDIT: and if so, what are they)?   * Are there any examples where the probability doesn't change after seeing the particular interval $[a,b]$ (i.e. the probability that the true parameter lies in $[a,b]$ is still 95%)?   * How (and why) does the probability change in general after seeing the computer spit out $[a,b]$? **[Edit]** Thanks for all the great answers and helpful discussions!


title: Pairwise comparisons for a regression with sandwich estimates (in R)

text: **The question in short** I run a robust regression in R and made a boxplot with only one of the predictor variable grouping the response variable. On this boxplot I'd like to add some information about the statistical model. What information (and how to display it (it is not an issue of programming)) would you suggest me to provide? * * * **The developed question** I have a several predictors: two non-ordinal predictor and one continuous predictor (below coded in R)               set.seed(81)     pred1 = rep(c('Car', 'Bike', 'Train', 'Airplane'), 6)     pred2 = rep(c('High', 'Low', 'Middle'), 8)     pred3 = rnorm(24)     resp = c(rnorm(12, sd = 1), rnorm(12, sd = 5))      `resp` is the response variable. I ran a robuste regression with sandwich estimates:               require(sandwich)     require(lmtest)               m = aov(resp ~ pred1 + pred2)     coeftest(m, sandwich)          t test of coefficients:                      Estimate Std. Error t value Pr(>|t|)       (Intercept) -0.49642    0.73911 -0.6716  0.51034       pred1Bike    1.55917    1.16568  1.3376  0.19769       pred1Car     1.23873    1.24080  0.9983  0.33135       pred1Train   2.50882    0.91468  2.7428  0.01338 *     pred2Low     0.11613    1.00540  0.1155  0.90932       pred2Middle  0.51476    0.90924  0.5661  0.57829       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      And I boxploted the groups for `pred1`:               require(ggplot2)     ggplot(data.frame(pred1, resp), aes(x=pred1, y=resp)) + geom_boxplot()      ![enter image description here](http://i.stack.imgur.com/8PnBDm.png) On this plot I'd like to add some letters in order to indicate the groups which are statistically similar (p.value < 0.05) as discussed here. Something like this:               ggplot(data.frame(pred1, resp), aes(x=pred1, y=resp)) + geom_boxplot() + annotate('text', x=1:4, y=6, label=c('a','b','a','b'), size = 8, color='red')      ![enter image description here](http://i.stack.imgur.com/wBvF3m.png) **My question is:** How can I find these p.values for the pairwise comparisons with my robust regression? I can do what follows where `m` is a simple aov model:               TukeyHSD(m)      But the following doesn't work:               TukeyHSD(coeftest(m, sandwich))      I might missunderstand what are these pairwise comparisons and what the results I currently have mean! Please let me know if you feel this! **The aim of my question is for me to understand what is the best way to display the results of my statistical model on a boxplot.** Note: the variables `pred2` and `pred3` are used to withdraw some part of variance that I don't want to be accounted to the effect of `pred1` (as `pred1`, `pred2`and `pred3` are correlated in my case). Therefore, I guess it is better not to run simple pairwise t-test in order to get the p.values I'd like to add at the top of each boxplot.


title: Panel data forecasting from Arellano-Bond GMM estimation

text: I want to come up with predictions of final energy demand per capita (fe) for a panel of countries. Explanatory variables are GDP per capita (gdp) and population density (pop) -- all variables are expressed in natural logarithms. I use _R_ and the _plm_ package for estimation. My model is based on the following equation and I use the _pgmm_ function of _plm_ , which does an Arellano-Bond estimation: $\text{log(fe)}_{it} = \delta \, \text{log(fe)}_{i,t-1} + \beta_1 \,\text{log(gdp)}_{it} + \beta_2 \,\text{log(gdp)}_{it}^2 + \beta_3 \,\text{log(pop)}_{it} + \beta_5 \, \text{log(year)}$ Now I want to forecast future values of final energy demand using scenarios for GDP per capita and population density and a logarithmic (sic!) time trend. If I understand the procedure correctly, this should work like this: $\text{fe}_{it} = exp(\delta \, \text{log(fe)}_{i,t-1} + \beta_1 \,\text{log(gdp)}_{it} + \beta_2 \,\text{log(gdp)}_{it}^2 + \beta_3 \,\text{log(pop)}_{it} + \beta_5 \, \text{log(year)}) + exp(s^2/2)$ (where $s^2$ is the mean squared error of the residuals.). I know simply plug in my scenario values and the final energy demand of the previous year, the numbers that I get are not even close to what I expected but rather off by several (dozens) orders of magnitude. For convenience, I leave out the term $exp(s^2/2)$, but since this is positive, it would only make my problem even larger. What am I doing wrong?


title: Correlation coefficients for ordered data: Kendall's Tau vs Polychoric vs Spearman's rho

text: It seems like for managing with ordered measurements researchers **usually** deal with Polychoric Correlation. (For example, for making matrix before doing Factor Analysis.) Why so? Kendall Tau Rank Correlation Coefficient and Spearman's rank correlation coefficient are also suitable for ordered data. Any points 'pro' and 'contra' for these correlation coefficients are welcome.


title: Using Univariate ANOVA with non-normally distributed data

text: If my data are non-normally distributed and I'm conducting a 2x2 ANOVA, what can I do to correct for this problem so I can report the main effect and interaction output appropriately? Only one finding is significant (one of the main effects). I've read that bootstrapping cannot be applied to a univariate ANOVA in SPSS...I also tried this out in SPSS myself and could not obtain any bootstrapped output. **More information:** I used a 2 x 2 (Age [6 years, 7 years] × Educational Classification [Group A, Group B]) univariate ANOVA to explore the hypothesis that test performance (percentage of items correct) would increase with age for children in Group B but would not differ significantly across age for students in Group A. I transformed the percentages into arcsine values to meet the assumption of a continuous dependent variable. My sample is small, so the sample size for each cell created by the 2 x 2 are as follows: Group A 6-year-olds = 10, Group A 7-year-olds = 13, Group B 6-year-olds = 20, and Group B 7-year-olds = 14. I know that scores are non-normally distributed across Group A and across Group B due to most students obtaining high scores on the test (and I ran tests of normality and looked at the Q-Q plots). Group A, which is also the smaller group of the 2, has larger variance. But should I run normality tests for all 4 groups created by the 2 x 2? I found a Main Effect for Educational Classification but not for Age, and there was no interaction. I just want to know if my findings are valid given the non-normality, or if I should find a way to fix this problem in SPSS (e.g., bootstrapping).


title: Ranking randomly distributed variables based on decreasingness

text: For some research I am conducting, I have run into the following problem, which I must solve repeatedly: My data consists of sequences of observations $\\{O_1,O_2,\ldots,O_n\\}$, each with a different normally distributed error. The sequences are variable in length, ranging from 5 to 44 observations. **I need a ranking system that expresses how much a given sequence of observations is _decreasing._** * * * Here's what I've come up with so far: Since everything is normally distributed, it's easy to compute $P(O_1 > \cdots > O_n )$ or $\prod_{i=1}^{n-1}P(O_i > O_{i+1})$ via Monte Carlo. So my first naïve way of tackling this problem is to try these two scores. There are problems with this approach. In my data, the observations may correspond to an exponentially decreasing quantity, or follow a power law. In these asymptotically convergent cases, for large $n$ these scores could be very low. What's more, a short sequence with wide error bars could easily get a higher score. So a good ranking system should punish short, ambiguous sequences and forgive long, asymptotically convergent sequences. Thinking of the ambiguous case, one can show that if all of the observation are i.i.d, then regardless of the distribution:   * $\displaystyle P(O_1 > O_2 > \cdots > O_n) = \frac{1}{n!}$ and    * $\displaystyle\prod_{i=1}^{n-1}P(O_i > O_{i+1}) = \frac{1}{2^{n-1}}$ From these observations, I was inspired to come up with the following rankings:   * $P(O_1 > \cdots > O_n ) \times {n!}$   * $\displaystyle\prod_{i=1}^{n-1}P(O_i > O_{i+1}) \times 2^{n-1}$ I can't say anything about the statistical properties of these rankings, as they are my own invention and I can't find any papers on this problem. Also, I have yet to implement anything, since I only encountered this problem very recently. I will update this question with my developments. Input and criticisms are welcome! Thanks everyone in advance!


title: Difference between panel data & mixed model

text: I would like to know the difference between panel data analysis & mixed model analysis. To my knowledge, both panel data & mixed models use fixed & random effects. If so, why do they have different names? Or are they synonymous? I've read the following post, which describes the definition of fixed, random & mixed effect, but doesn't exactly answer my question: What is the difference between fixed effect, random effect and mixed effect models? I would be grateful if somebody could refer me to a brief (about 200 page) reference on mixed model analysis.Just to add, I would prefer Mix modelling reference irrespective of software treatment. Mainly theoretical explanation of Mix modelling.


title: Two way ANOVA or ANCOVA

text: I need to show that BMI of adults is neither age nor gender specific. Should I do ANCOVA or randomly make age groups and do two way ANOVA?


title: Which test to use to check if a possible confounder impacts a 0 / 1 result?

text: I've been given a task with the following question:   1. Investigate whether or not the type or treatment (0 or 1) has an impact on the result (0 or 1)    2. The same as 1), but keeping in mind another variable which is either 0,1 or 2 may be a confounder for the result.  How would I check whether or not it has an impact? For question one, I think I can use a chi-squared test, is this correct? And for question two, I'm not sure how to include the confounder in the test.


title: Determining smoothing parameter in HP filter for hourly data

text: I'm trying to determine a smoothing parameter for the Hodrick-Prescott filter. I've seen that there are papers on the topic but they are far too advanced for my comprehension. If I have a data set, $X$, what are the steps to take? My data is hourly. Using MATLAB I know the smoothing parameter values for monthly data etc., but how can I compute it for hourly data?


title: Specifying the number of clusters in nearest neighbor clustering

text: I've got a distance matrix between examples. I want to cluster them into _m_ clusters with a nearest neighbor algorithm which works like this:               1. Set i = 1 and k = 1. Assign example x_1 to cluster C_1.     2. Set i = i + 1. Find nearest neighbour of x_i        among the patterns already assigned to clusters.        Let d_n  denote the distance from x_i to its nearest neighbour.        Suppose the nearest neighbour is in cluster n.     3. If d_n less than or equal to t then assign x_i to C_n where t is the         threshold specified by the user. Otherwise set k = k+1 and assign x_i  to a              new cluster C_k.      How could I adapt this algorithm so I could specify how many clusters I want? Is anybody aware of an existing R implementation of nearest neighbour clustering?


title: Limitation on group size with non-parametric permutation tests

text: I would like to know if there is any limitation in the size of the groups when you are implementing non-parametric permutation tests. My case belongs to the realm of medical imaging. I am comparing one patient image versus a group of control images. For each image (including the patient's), I calculate a maximum T-statistic. This is my 'first level' metric. In my second level design I use my first level metric to compare the patient against the distribution of the first level metric in the group of controls. I obtain a T-score. Then I permute the labels (patient, control) of my images and obtain again a T-score for that permutation. I run all permutations where I have one image in the patient group and the rest of the images (30) in the control group. If I have 1 patient and 30 controls, this will yield 31 permutations with a configuration of 1 patient and 30 controls . After building the permutation distribution with this procedure, I notice that the score for my initial configuration (the one I am testing) does not surpass the 0.05 threshold obtained from the permutation distribution. This is kind of expected, since I don't have enough permutations to have a low enough p-value. However, (and this is a big however), when I ran ALL possible permutations ($2^{31}$) assigning the label patient or control to each subject, I notice that my initial configuration statistic surpasses the 0.05 threshold. Can I do this, though? Can I run all of the possible permutations when my original design is 1 versus many? Is there something that I am missing? Or does my procedure seem reasonable?


title: PhD vs MS in biostatistics

text: What are the benefits of getting a PhD vs a masters in biostatistics (besides being a professor)? Will the first job will be the same for both levels? Does the PhD in biostatistics really make a difference in industry?


The empirical Bayesian method is a new concept to me. It raises my interest, because it offers a different philosophical and methodological perspective of statistical analysis. From my limited knowledge, I am not sure if it is useful/powerful, or inferior to alternative methods, and if it is worth learning (and to what degree) compared to learning other useful topics in Bayesian analysis? How do you think about it? What are pros and cons of empirical Bayesian methods?


title: Simple experimental design - should I counterbalance?

text: I'm designing a pretty simple experiment that goes like this. Participants will be shown a series of stimuli and after viewing each one they will answer a few questions where they will make judgments about the stimulus - all Likert items. There are two kinds of stimuli. Probably obvious, but the hypothesis is that there will be a difference between answers for A vs B stimuli. There will be 30 or so stimuli, with an equal number of A and B stimuli. All participants will see all the stimuli (within-subjects). I'm wondering if there would be a benefit to counterbalancing the order in which they receive the items, vs just showing everyone the same randomized sequence of stimuli (which is easier to setup). If there's a better method I need to consider, I'd be interested in hearing about it. I also looked into blocking designs, but this is so simple that I don't think those apply here. I'm planning to analyze with t-tests or Mann- Whitney-Wilcoxon.


title: Regression with error in covariates

text: I am looking for some advice for a colleague who is dealing with regression models for which it is know, that the continuous covariate of interest $X_1$ was measured with error. More precisely, we know only know that the true $X_1$ lies within an interval $(X_L, X_R)$. My first idea was to illustrate/describe the effect of the measurment error by fitting three models of the form \begin{align*} \text{Model 1:} \quad Y & = \beta_0 + \beta_1X_M + \beta_2 X_2 +... \\\ \text{Model 2:} \quad Y & = \beta_0 + \beta_1X_L + \beta_2 X_2 +... \\\ \text{Model 3:} \quad Y & = \beta_0 + \beta_1X_R + \beta_2 X_2 +... \\\ \end{align*} where $X_M$ is the midpoint of the interval. What other regression methods are available that take into account that the true value of $X_1 \in (X_L, X_R)$? Specifically, when using linear regression models, logistic regression models, or generalized linear mixed models. There is some information on the Wikipedia-Site about Errors-in-variables models (http://en.wikipedia.org/wiki/Errors-in-variables_models), but I'm not sure if those methods can be applied in this case.


title: How to separate out the regression effect vs treatment effect without a control group?

text: I'm looking at a dataset that has pre-post test measurements on users' stress, depression and anxiety levels collected from a website's online health assessment. On average, the healthier participants at baseline got worse over time, and the sicker participants at baseline got much better, and the middle group gets a little better. There's definitely a regression effect going on here, but also a treatment effect too. As this data was collected based on website usage, there isn't really a control group (all of the "post" measurements come from people that have used the online program). There are probably ways that I could synthesize a control group using the people who I can guess didn't make much use out of the treatment (based on number of logins or length of time between logins), but is there a way to separate out the treatment effect from the regression effect when you can't use difference-in-difference techniques using a control group or anything like that? Thanks!


title: Poisson distribution vs multiplying probabilities

text: I am a TA for a stats course for engineers, and I had a really good question from a student today, which I don't know the answer to. We were going through the following word problem: "4 computers run continuously for the Toronto Stock Exchange. The probability of a computer failure in a day is estimated at 5%. Assuming differing computers fail independently, what is the probability that all 4 computers fail in a day?" Since the sampling takes place over an interval, the way I would approach this is using the Poisson distribution, with the average number of computers failing on a day $\equiv\lambda = 0.05$. If four computers fail, then $k = 4$. Thus, \begin{align*} P(k; \lambda) &= \frac{\lambda^{k} e^{-\lambda}}{k!} \\\ P(k=4; \lambda = 0.05) &= \frac{0.05^{4} e^{-0.05}}{4!} \\\ & = 2.477\times 10^{-7} \end{align*} However, a student asked why it would not be appropriate to just multiply the probability of each computer failing. Since the probability of each computer failing each day $\equiv p = 0.05$, and since each computer failure is independent, he argued that, \begin{align*} P(k=4) &= p^4 \\\ &= 0.05^4 = 6.25\times 10^{-6} \end{align*} Which one of these approaches is wrong given the question? And why? What underlying assumption of the wrong approach is violated by the question? Thank you for your help. UPDATE: I left out some information in the problem the first time this was posted, and I apologize.


title: Inference for between group differences (with non-linear time component, outliers and mixed effects)

text: **TLDR** : How can I perform inference for the between group differences in a possibly logistic growth with time in the presence of outliers, unequal measurement times and frequency, bounded measurements and possible random effects on individual and per study level? I am attempting to analyse a dataset where measurements for individuals were made at different time points. Measurements start low at time 0 and follow (very roughly) a logistic growth pattern with time. I am trying to establish if there are differences between two groups of individuals. The analysis is complicated by the following factors:   * The effect of time is non-linear, so either a non-linear logistic regression (biologically plausible, but not particularly well fitting) or a non-parametric regression seem appropriate   * There are massive outliers, so regression using the sum of squared residuals seems off the table. Quantile regression seems appropriate.   * Random effects may be appropriate on a per individual and per study level. Mixed effects models seems appropriate.   * Measurement times, number of available measurements and end of monitoring differ between individuals. Survival analysis techniques seem appropriate. Possibly also applying weights equal to 1 / number of observations for individual.   * Measurements are bounded below at 0 and while there is no obvious boundary above, arbitrarily high measurements seem biologically implausible. However, quite a few individuals have some measurements of zero (partly due to the measurement accuracy of the device).   * A few models I tried so far failed to fit, usually with an unhelpful error related to the numerical procedure. This leads me to believe that I will need a reasonably robust method able to deal with this somewhat ugly dataset.   * Finally, I want to produce inference of the form "group 1 has faster growth than group 2" or "group 1 has a higher asymptotic level than group 2". What I have tried so far (all in R) - I was aware that most of the below are not particularly appropriate for the dataset, but I wanted to see which models could actually be fitted without numerical errors:   * Non-parametric regression using crs in the crs package. Nicely produces a curve reasonably close to logistic growth for most of the time period with some strange behavious toward the end of the monitoring period (where there are fewer observations). Using individuals as fixed effects reveals some outliers. Using the variable of interest as fixed effects shows some difference. However, I am not sure if there is any way to assess fits and do inference on a model this complex.   * Non-linear mixed effects regression using nlme in package nlme and SSlogis. Gradually building up the model with update() works reasonably well. Getting too complex with the fixed effects or the random effects leads to convergence failure. Since there is no reason to believe that residuals are normally distributed, it seems a bad idea to pursue this further. Edit: I have recently become aware that it is possible to specify autocorrelated residuals in nlme. However, at the moment it seems I cannot even get fixed weights to work. Advice on the correct syntax is welcome.   * Non-linear mixed effects regression using nlmer in package LME4 and a custom likelihood for the logistic growth model. Works fairly well, but standard errors on the fixed effects get massive, probably due to the outliers. I also have the slight suspicion that some of the models fail to fit without error, as I sometimes get tiny random effects (about 10^10 smaller than with slightly simpler models). Since there is no reason to believe that residuals are normally distributed, it seems a bad idea to pursue this further.   * Non-linear quantile regression using nlrq in package quantreg and SSlogis. Fits reliably and quickly, but percentile lines intersect. This means that an area containing 90% of the data is not fully contained in an area containing 95% of the data.   * Non-parametric quantile regression using the LMS method with package VGAM. Even trivial models failed with obscure errors using this dataset. I believe the number of zeros in the dataset and / or the large range of the data while also getting close to zero may be an issue.   * To complete this list, I should probably also mention the lqmm package for Linear Quantile Mixed Models, which I have not used yet. While the package cannot use non-linear models as far as I know, transforming the time variable may produce something reasonably close. I would appreciate feedback if these or any other method might be used to produce reasonably robust inference in this scenario. Maybe regression is not needed at all and another, possibly simpler method is sufficient. I'd be happy to provide an example dataset, if required, but think this question might also be of interest beyond the current dataset.


title: Kaplan Meier curve with different time lengths across groups

text: I'm not too familiar with survival analysis, so this may be a somewhat basic question. I am interested in the time between two events [$A$ - medication] and [$B$ - time of first treatment], the latter which can be equal to $A$. Specifically, I want to know if this varies across years. The subjects in question have at least 5 years of data, and I have data from 2001-2011. Because of this, tests on summary statistics won't do (subjects in 2011 have the capacity to have up to 11 years while subjects in 2006 only have the capacity to have 6). So instead, I'm considering the value $A-B$ to be censored if $B$ falls within one year of the subject entering the dataset, and then constructing a KM curve. However, the time length covered by each group will still end up different, like so: ![Un-truncated & censored](http://i.stack.imgur.com/zjKVx.png) By visual inspection, the survival curves are different, though only the 2011 vs any other group comparisons are significant by the logrank test. Alternatively, I can truncate time differences greater than 5 years (considering them censored): ![Truncated & censored](http://i.stack.imgur.com/NJ6Bo.png) This results in the exact same curve for years <= 5 as expected, however, the logrank test returns more significant differences (2010 vs any other group). I'm not sure which of these setups is more appropriate.


title: Sample variance Fisher consistency

text: I am reading conflicting references regarding the Fisher consistency of the sample variance. $$s^2 = \frac{\sum_i^n (x_i - \bar{x})^2}{n}$$ Could anyone explain me how to proof whether $s^2$ is Fisher consistent or not? The discussion page on Wikipedia was inconclusive as well.


title: Help in understanding maximum likelihood estimation

text: I am having trouble in understanding link how the 5.4 has come from and the line **" Since the distribution has to be equal to 1 at $r_u$, it must be divided by $\tilde{\phi} \exp (-d \tau K_{2}) r_u^{\nu}$"** what does "it" mean and where from is that distribution coming? After, the distribution how is Equation (5.4) realized? I have found no references, saying this distribution is applied in Takens formula. Please help.


title: Comparing standard deviations of two dependent samples

text: As far as I know you can do a F-test ($F = s_1^2/s_2^2$) or a chi-squared test ($\chi^2 = (n-1)(s_1^2/s_2^2$) for testing if the standard deviations of two independent samples are different. But does this also hold for dependent samples?