title: Is randomization reliable with small samples?

text: Jerome Cornfield has written: > One of the finest fruits of the Fisherian revolution was the idea of > randomization, and statisticians who agree on few other things have at least > agreed on this. But despite this agreement and despite the widespread use of > randomized allocation procedures in clinical and in other forms of > experimentation, its logical status, i.e., the exact function it performs, > is still obscure. Cornfield, Jerome (1976). "Recent Methodological Contributions to Clinical Trials". American Journal of Epidemiology 104 (4): 408–421. Throughout this site and in a variety of literature I consistently see confident claims about the powers of randomization. Strong terminology such as "it **eliminates** the issue of confounding variables" are common. See here, for example. However, many times experiments are run with small samples (3-10 samples per group) for practical/ethical reasons. This is very common in preclinical research using animals and cell cultures and the researchers commonly report p values in support of their conclusions. This got me wondering, how good is randomization at balancing confounds. For this plot I modeled a situation comparing treatment and control groups with one confound that could take on two values with 50/50 chance (eg type1/type2, male/female). It shows the distribution of "% Unbalanced" (Difference in # of type1 between treatment and control samples divided by sample size) for studies of a variety of small sample sizes. The red lines and right side axes show the ecdf. **Probability of various degrees of balance under randomization for small sample sizes:** ![enter image description here](http://i.stack.imgur.com/H3n20.png) Two things are clear from this plot (unless I messed up somewhere). 1) The probability of getting exactly balanced samples decreases as sample size is increased. 2) The probability of getting a very unbalanced sample decreases as sample size increases. 3) In the case of n=3 for both groups, there is a 3% chance of getting a completely unbalanced set of groups (all type1 in the control, all type2 in the treatment). N=3 is common for molecular biology experiments (eg measure mRNA with PCR, or proteins with western blot) When I examined the n=3 case further, I observed strange behaviour of the p values under these conditions. The left side shows the overall distribution of pvalues calculating using t-tests under conditions of different means for the type2 subgroup. The mean for type1 was 0, and sd=1 for both groups. The right panels show the corresponding false positive rates for nominal "significance cutoffs" from .05 to.0001. **Distribution of p-values for n=3 with two subgroups and different means of the second subgroup when compared via t test (10000 monte carlo runs):** ![enter image description here](http://i.stack.imgur.com/PBxSw.png) **Here are the results for n=4 for both groups:** ![enter image description here](http://i.stack.imgur.com/yAxux.png) **For n=5 for both groups:** ![enter image description here](http://i.stack.imgur.com/MbM4w.png) **For n=10 for both groups:** ![enter image description here](http://i.stack.imgur.com/9Dt0R.png) As can be seen from the charts above there appears to be an interaction between sample size and difference between subgroups that results in a variety of p-value distributions under the null hypothesis that are not uniform. So can we conclude that p-values are not reliable for properly randomized and controlled experiments with small sample size? **R code for first plot**               require(gtools)          #pdf("sim.pdf")     par(mfrow=c(4,2))     for(n in c(3,4,5,6,7,8,9,10)){       #n<-3       p<-permutations(2, n, repeats.allowed=T)            #a<-p[-which(duplicated(rowSums(p))==T),]       #b<-p[-which(duplicated(rowSums(p))==T),]            a<-p       b<-p            cnts=matrix(nrow=nrow(a))       for(i in 1:nrow(a)){         cnts[i]<-length(which(a[i,]==1))       }                 d=matrix(nrow=nrow(cnts)^2)       c<-1       for(j in 1:nrow(cnts)){         for(i in 1:nrow(cnts)){           d[c]<-cnts[j]-cnts[i]           c<-c+1         }       }       d<-100*abs(d)/n            perc<-round(100*length(which(d<=50))/length(d),2)            hist(d, freq=F, col="Grey", breaks=seq(0,100,by=1), xlab="% Unbalanced",            ylim=c(0,.4), main=c(paste("n=",n))       )       axis(side=4, at=seq(0,.4,by=.4*.25),labels=seq(0,1,,by=.25), pos=101)       segments(0,seq(0,.4,by=.1),100,seq(0,.4,by=.1))       lines(seq(1,100,by=1),.4*cumsum(hist(d, plot=F, breaks=seq(0,100,by=1))$density),             col="Red", lwd=2)          }      **R code for plots 2-5**               for(samp.size in c(6,8,10,20)){       dev.new()       par(mfrow=c(4,2))       for(mean2 in c(2,3,10,100)){         p.out=matrix(nrow=10000)              for(i in 1:10000){                d=NULL           #samp.size<-20           for(n in 1:samp.size){             s<-rbinom(1,1,.5)             if(s==1){               d<-rbind(d,rnorm(1,0,1))             }else{               d<-rbind(d,rnorm(1,mean2,1))             }           }                p<-t.test(d[1:(samp.size/2)],d[(1+ samp.size/2):samp.size], var.equal=T)$p.value                p.out[i]<-p         }                   hist(p.out, main=c(paste("Sample Size=",samp.size/2),                            paste( "% <0.05 =", round(100*length(which(p.out<0.05))/length(p.out),2)),                            paste("Mean2=",mean2)         ), breaks=seq(0,1,by=.05), col="Grey", freq=F         )              out=NULL         alpha<-.05         while(alpha >.0001){                out<-rbind(out,cbind(alpha,length(which(p.out<alpha))/length(p.out)))           alpha<-alpha-.0001         }              par(mar=c(5.1,4.1,1.1,2.1))         plot(out, ylim=c(0,max(.05,out[,2])),              xlab="Nominal alpha", ylab="False Postive Rate"         )         par(mar=c(5.1,4.1,4.1,2.1))       }          }     #dev.off()


title: Differences between clustering and segmentation

text: I have read about piecewise aggregate approximation (PAA) mining time series data, sliding window, top down and bottom up approaches for time series segmentation but these are applicable to single dimension time series. What are the techniques for multi dimensional time series segmentation? Can Gaussian mixture model or K means clustering be used for segmentation? If so then what is the difference between segmentation and clustering? I want to apply segmentation to motion data represented by joint angle trajectories so that each segment corresponds to an action. How do I do it? Paper1 applies Gaussian Mixture model (GMM) for clustering but that would imply losing time information. Paper2 also uses GMM for clustering and segmenting the multi dimensional data where each pixel is represented by a pdf. My objective is how to segment the multi dimensional time series so that the temporal information is retained and each subsequence can be used to represent an action and then classification can be used. There is a tool box Symbolic Aggregate Approximation which uses PAA to segment the data nad then labels each segment so obtained. I am not able to find a new technique for segmenting the data. Also, I am a bit confused when authors apply GMM for segmenting since it is a model based clustering algorithm. Therefore, the questions are :   1. What is the difference between segmenting and clustering   2. How do I segment motion time series data so that the temporal information is retained.   3. What are the algorithms for doing so - What are the techniques for multi dimensional time series segmentation?? Please provide links or ideas thank you.


title: How do I set up a regression analysis and verbally describe what I'm doing to find the following?

text: I'm looking to see if there is first: \- a main effect of the gender of the interviewer \- a main effect of the gender of the interviewee And then if significant: \- the main effect of each gender interaction/dyad: 4 separate main effects - male/male, male/female, female/female, female/male (bowles) \- the main effect of the interaction between the negotiation strategy used and 4 interaction/dyad makeups (16 interactions) I also need to control for demographic variables. These are my hypotheses: 1st Hypothesis: Effect of gender interaction 2nd Hypothesis - Interviewer men will treat women differently (Bowles) 3rd Hypothesis - women will receive lower salary outcomes The outcome variable is increased salary. Does that make sense? Let me know if you need more info.


title: Find the Distribution

text: Given $ Y_1, Y_2..Y_n$ are iid from a distribution with pmf,   $f(y) = a^{2}$ for $y=0$, $f(y) = 2a(1-a)$ for $y=1$ , $f(y) = (1-a)^{2}$ for $y=2$, where $0<a<1$. For large n, calculate the approximate distribution of a) $\sqrt {\bar{Y}}$ - Solution to part(a) posted as answer(awaiting confirmation) b) $\sqrt n ({\bar{Y}-\mu)}+\bar Y^2$ , where $\mu=E(Y_1)$ Could you please verify my solution for part (b) : By CLT $\sqrt n ({\bar{Y}-\mu)} \rightarrow N(0,\sigma^2)$ (convergence in probability) For $\bar Y^2$, applying delta method, $\bar Y^2 \rightarrow N(\mu^2,\frac{4\mu^2\sigma^2}{n^2})$ (converges in distribution) **{EDIT} - Can I say :** $\bar Y^2 \rightarrow \mu^2$ in probability where $\sigma^2 = Var Y$ and $Var \bar Y^2 = \sigma^2/n$ **Can I apply slusky theorem** , as one distribution converges in probability and other in distribution: By Slutsky theorem , $\sqrt n ({\bar{Y}-\mu)} + \bar Y^2 \rightarrow [\mu^2 + N(\mu^2,\frac{4\mu^2\sigma^2}{n^2})]$ Thanks!


title: What does error refer to in linear regression notation?

text: I regularly see linear regression models written in this notation: $y = a + \beta X + error$ I've never really pinned down what $error$ actually refers to. In the linear regression plotted below, does $error$ refer to the grey band, therefore the confidence interval for the prediction line? ![enter image description here](http://i.stack.imgur.com/RaVXN.jpg)


title: Is it possible to interpret regression output omitting the intercept?

text: Is it possible to interpret regression output when omitting the intercept? Can this omission be justified?


title: Identifying the population and samples in a study

text: I apologise for the silliness of this question. I ran an e-learning experiment using a class of undergrads. Participation was voluntary, so only half of the class participated. I know that my sample is the people who participated. Is the population all of the students in that class, or is the population all undergrads in similar course (or whatever I'm trying to generalise the study to)? If it's the latter, then what do I refer to the students in the class as? The cohort?


title: Question about Confidence Interval

text: You have graduated from DePaul with a degree in HCI. Congrats! Your new boss tells you to survey 25 of the company's customers, asking them to rate on a 5-point scale their satisfaction with the website. He wants the margin of error to be 0.25. You know from previous experiments that the standard deviation is 1.0. How large will your confidence level be? 70% 80% 90% 99% 99.9%


title: Using Principal Component Analysis to determine variable interaction

text: I'm characterizing a communication channel by measuring bit error rate (ber) as a function of 15 channel settings: ber = Func(settings[1:15]). The dataset is about 100,000 samples. I wonder if I can use PCA to explore interaction of settings and their effect on the bit error rate. Specifically:   1. How to prepare data for PCA; do I run PCA on {ber, settings[1:15]}, or only on settings[1:15] ? I've seen PCA used on "symmetric" data, such as images, text, etc.   2. What do projections of settings[1:15] on principal components tell me in this case   3. What do outlier datapoints in principal component space tell me.


title: Calculate variance of survival estimate generated from 'baseline' of coxph

text: **Short version:** calculate variance of estimator of survival for a new subject (new covariates) where the coefficients for the covariates have been calculated from existing data. Sorry this is so long, not sure quite where I'm going wrong. Feel free to edit if would be easier to have equations then code. This is really a follow-up question to this earlier one. I'm trying to replicate Klein & Moeschberger: Survival Analysis (2nd ed.), example 2, pg. 284. I've got the predictor for survival, but not the variance of the estimator. If this is already implemented in `package:survival` I'd be grateful, but I'd prefer to calculate it manually. The textbook appears to follow the same logic as this paper from 1979 so it should be possible to calculate this without having the above textbook to hand. First part is adapted from the answer to the earlier question by @ocram: **This part seems OK**               library(survival)     data(larynx, package="KMsurv")     ### make data.frame containing unique event times     ### and no. events at those times     t1 <- as.data.frame(table(larynx[larynx$delta==1, "time"]))     names(t1) <- c("t","e")     ### change factor to numeric     t1$t <- as.numeric(levels(t1$t))     ### fit intercept model     c2 <- coxph(Surv(time, delta) ~ 1, data=larynx)     ### use this to generate baseline hazard     ### i.e. all covariates = 0     b1 <- basehaz(c1, centered=FALSE)     ### only keep rows where events occurred     b1 <- b1[ b1[, "time"] %in% t1$t, ]      ### get dummy variables (for factors)     m1 <- model.matrix(c1)     ### store results (used later)     Wk <- matrix(nrow=nrow(t1), ncol=ncol(m1))     Q3 <- rep(NA, ncol(m1))     ### new covariates     Z0 <- c(60, 1, 0, 0)      Now try to find $Q1$, $Q2$, $Q3$ as per textbook. I have put this in a `for` loop to make it easier to read. Note this spans multiple code blocks to allow appropriate equations:               for (i in 1:nrow(t1)){     ### one row per time with at least one event      Now $b$ is the vector of estimators of the coefficients ( $p=4, k=1,2,3,4$ in this case) and $Z$ is the matrix of parameters from the `model.matrix`. So: \begin{equation} W(t_i;b) = \sum_{j \in R(t_i)} \exp \left( \sum_{k=1}^p b_k Z_{jb} \right) \end{equation}                   t1$W[i] <- sum(  exp( rowSums( sweep( m1[larynx$time >= t1$t[i], ],                                             MARGIN=2,                                             as.matrix(c1$coef), "*" )                                  ) )  )      (`sweep()` is used to facilitate matrix * vector multiplication as illustrated in this question). Then, with $e_i$ = no. events per time period, the estimated hazard at that time is: \begin{equation} \hat{h}_0(t) = \frac{e_i}{W(t_i;b)} \end{equation} and the cumulative hazard is: \begin{equation} \hat{H}_0(t) = \sum_{t <= t_i} \hat{h}_0 \end{equation} i.e.                     t1$H[i] <- t1$e[i] / t1$W[i]           t1$CH <- cumsum(t1$H)      The baseline survival function (i.e. all covariates = 0) is \begin{equation} \hat{S}_0(t) = \exp{ -\hat{H}_0(t) } \end{equation} and the estimate of survival for a subject with covariates $\textbf{Z0}$ is \begin{equation} \hat{S}(t|Z=\textbf{Z0}) = \hat{S}_0(t)^{ \exp (b \textbf{Z0} ) } \end{equation} i.e.                     t1$S0 <- exp(-t1$CH)           t1$S1 <- t1$S0^(exp(sum(c1$coef*Z0)))      Next, $Q1$ which is \begin{equation} Q1(t) = \sum_{t <= t_i} \frac{e_i}{W(t_i; b)^2} \end{equation} or                    t1$Q1[i] <- t1$e[i]/(t1$W[i]^2)      **Problem near here I think** Now, the tricky part. I need the estimator $W_k(t_i;b)$ where \begin{equation} W_k(t_i;b) = \sum_{j \in R(t_i)} Z_{jk} \exp(b\textbf{Z0}) \end{equation} which (I think) is                     Wk[i, ] <- colSums(  sweep(m1[larynx$time >= t1$t[i], ],                                MARGIN=2,                                exp(c1$coef * Z0),                                "*") )      This is used to generate $Q3(t, \textbf{Z0})_k$, a vector of length $k$ \begin{equation} Q3(t, \textbf{Z0})_k = \sum_{t <= t_i} \left( \frac{W_k(t_i;b)}{W(t_i;b)} - \textbf{Z0}_k \right) \left( \frac{e_i}{W(t_i;b)} \right) \end{equation} like so                    ### still inside first 'i' loop here          for (k in 1:ncol(m1)){              Q3[k] <- sum ( (Wk[t1$t <= t1$t[i], k] / t1$W[i] - Z0[k] )  * t1$H[i] )          }      then, using the covariance matrix $\hat{V}$ which (I think) is `c1$var`: \begin{equation} Q2(t, \textbf{Z0}) = Q3(t; \textbf{Z0})^T \hat{V} Q3(t; \textbf{Z0}) \end{equation} i.e.                   Q2 <- t(Q3) %*% c1$var %*% Q3      and use this to generate the variance \begin{equation} \hat{V} = [\hat{S}(t|Z=\textbf{Z0}]^2 [Q1(t) + Q2(t;\textbf{Z0})] \end{equation} and the standard error $SE$ \begin{equation} SE = \frac{ \sqrt{\hat{V}}}{\sqrt{n}} \end{equation}                   t1$var <- t1$S1^2 * (t1$Q11 + Q2)             t1$SE <- sqrt(t1$var)/sqrt(nrow(m1))     } # first 'i' loop ends here      Now the result of all this should be that at $t=5$ `t1$S1` should be around $0.66$ which we can see from `print(t1)` or from:               Z0df <- data.frame(age=60, stage=2)     cbind(survfit(c1, newdata=Z0df)$time, survfit(c1, newdata=Z0df)$surv)      but the $SE$ should be around $0.106$. For another example: the $SE$ with `Z0 <- c(60,0,1,0)` should be $0.095$.


title: Interpreting the graph of the following confidence interval for the proportion of binomial in R

text: #Recursion of Binomial            bino <- function(n,p){            px <- (1-p)^n            X <- 0            Fx <- px            U <- runif(1)            while(Fx<U){               X<- X+1               px <- ((n-X+1)*p*px)/(X*(1-p))               Fx <- Fx+ px            }            return(X)         }              #inputs          N <- 10000         n <- 10         p <- 0.7              set.seed(100)         X <- rep(0,N)              for(i in 1:N) X[i]=bino(n,p)              phat <- rep(0,n+1)              for(i in 0:n)phat[i+1]=(sum(X==i))/N              phat.CI <- 1.96*sqrt(phat*(1-phat)/N)                   # plot output              plot(0:n, dbinom(0:n, n, p), type="h", xlab="x", ylab="p(x)")         points(0:n, phat+phat.CI, pch=3)         points(0:n, phat-phat.CI, pch=3)      _observing the graph it is seen that the interval for the proportion when X=10 is wider than that of others. What does it indicate?_ After running the code here is the graph: ![enter image description here](http://i.stack.imgur.com/J9cCE.png)


title: binary data with no correlation but with p-values equal to real data

text: For doing parallel analysis in my real data to determine the number of its dimension, i need to simulate binary data matrix with n(number of case) by k(number of items). the p-values of these simulated data must be equal to the corresponding p-values of items in sample real data with n cases and k items. but unlike the real data that there is correlation between items, in simulated data the correlation of items must be zero.   i want to generate simulate data by the following syntax. it is just an example. let me know the correctness of my procedure:               library(MASS)     set.seed(99)     n<-100 #NUMBER OF PERSONES     K<-5 #NUMBER OF ITEMS     p1<-rep(0.2,n)# p-values of first item is 0.2     p2<-rep(0.8,n)     p3<-rep(0.1,n)     p4<-rep(0.2,n)     p5<-rep(0.3,n)     p<-cbind(p1,p2,p3,p4,p5)     X<-matrix(rbinom((n*K),1,p),n,K)#datafile     write.table(X,file="MATH.txt",row.names=FALSE,col.names=FALSE)


title: Statistic not sensitive to coefficients

text: I'm looking for some kind of function I can compute on a bunch of numbers. The idea is that if I increase _all_ numbers by some coefficient (say 2) the value of this should not change. e.g. f(2,4,6,9)= f(4,8,12,18) * * * Edit: I know I can just divide all the numbers by the largest and then compute any other statistic. I'm more looking for some function that doesn't care about the actual numbers but the relative differences betwen them. The problem with just dividing everything by the largest is that it would make this _extremely_ sensitive to outliers or bad data


title: What's the distribution of these data?

text: I got the data, and plot the distribution of the data, and use the qqnorm function, but is seems doesn't follow a normal distribution, so which distribution should I use to discribe the data? Empirical cumulative distribution function ![enter image description here](http://i.stack.imgur.com/xnQpe.png) ![enter image description here](http://i.stack.imgur.com/oAG0d.png)


title: Issues while implementing stochastic gradient descent

text: Actually, I am trying to learn a model using stochastic gradient descent. I am using just a small subset of my entire data at a time to make each iteration. As far as I know, stochastic gradient descent can fluctuate but still moves towards the maximum point. I am outputting the value of the objective function and I could see that it's fluctuating between 10 and 43 and not getting above it. My assumption is even though it fluctuates it slowly keeps on approaching towards the maximum value. It's not guaranteed that it will reach the global optimum and then terminate. It can overshoot and then undershoot and then keep on fluctuating there. But in my case it is not showing an increasing trend. At one time the value of objective function is 43. The very next time it is 11 and then 35,27. It hasn't gone above 40. Any suggestions? What could have gone wrong?


title: Exact p-value for any nonparametric test

text: After having found the R-package `coin`, I have a question. For example, to get an exact p-value, I could use the function `wilcox.exact` in the `exactRankTests` R package. But if I do other nonparametric test like the `Kruskall-Wallis` \- is it possible to get the exact p-value then? Do You know if there is a R-package in which I can get the exact p-value of any nonparametric test because in the one above I have just found how to get the exact p-value only for the Wilcoxon?


title: How does JAGS deal with state space models?

text: I am trying to use JAGS to deal with the following multivariate state space model. $Y_t=X_t\theta_t+\epsilon_t$ $\theta_t=\theta_{t-1}+\nu_t$ JAGS code is neat but JAGS is running too slow when I am monitoring several hundred variables in the model. I decide to use forward filtering backward sampling algorithm (FFBS) to write my own code. I am wondering how JAGS deals with state space model and whether FFBS is more effcient than JAGS.


title: How to get the P-value from a known Friedman's chi-squared value

text: Provided I have the value of the Friedman's chi-squared statistic, is there any way to get the corresponding P-value without copying the whole table of critical values in my Python program?


title: Differing definitions of Matrix Condition Number

text: When describing the condition number of a correlation matrix I have seen is described as the ratio of the singular values (from singluar value decomposition). I have also seen it described as the ratio of the eigenvalues (from PCA analysis). Can both be correct (under say different norms)? If so which one would be considered standard?


title: Matching file-names with movies from imdb

text: I have 10 millions of filenames and I need to determine whether filename corresponds to some movie or video game. For this I initially thought to compare filenames using Levenshtein distance with imdb list ftp://ftp.funet.fi/pub/mirrors/ftp.imdb.com/pub/ but it turns out that direct compassion requires $O(nm)$ where n is number of filenames and m is number of items in imdb database, both of which in millions range. I think there should be a way to organize imdb list in some tree like suffix structure to speed up comparisons to something like $O(n \log m)$. Any other advices and ideas are very much appreciated, too.


title: What regression model to use when independent variables are percentages to predict % outcome?

text: "Independent" variables: time spent (% at work, % sleeping, % exercising), body mass composition (% fat, % muscle, % bone) Dependent variable: Smoker (1) or Non-Smoker (0) What kind of regression model should I use when subsets of the "independent" variables are percentages and are therefore not completely independent of each other?


title: R: Which distribution to use with gbm for gamma distributed data?

text: When I use GLMs I can use the option `family="Gamma"` for analysing data consisting of positive real numbers. Also package `gbm` provides a large number of distributions to choose from, but there's none that matches the gamma distribution. Which distribution should I choose?


title: Fraction confidence intervals for small sample size

text: What is 95% confidence interval for the fraction of 5 successes out of 12 trials? Is it possible to compute confidence intervals for a sample of this size?


title: Prerequisites for AIC/BIC model comparison

text: I have a question about model selection when using AIC/BIC. So, if two model structures are totally different, can I still directly apply AIC and BIC? Also, for a hierarchical model, how to compute the total number of parameters when computing AIC and BIC? Thanks.


title: Extract Random Effects from plm model

text: The `ranef` function provides a way to extract the random effects coefficients from an `lme` model. **What is the function to extract random effects for a`plm` model?** Find below R code from the JSS paper (http://www.jstatsoft.org/v27/i02/) discussing panel data econmetrics with the `plm` pacakge. The same model is estimated by using `lme` and `plm`. However, a function to extract random effects in not mentioned for the latter.               require(plm)     require(nlme)          data("Grunfeld", package="plm")          # Estimate the same model with plm and lme     reGLS <- plm(inv~value+capital, data=Grunfeld, model="random")     reML <- lme(inv~value+capital, random=~1|firm, data=Grunfeld)          # Check if coefficients are the same     coef(reGLS)     summary(reML)$coef$fixed          # Extract random effects for lme model     ranef(reML)


title: Factor analysis - CATPCA combined with conventional PCA

text: I have some concerns regarding factor analysis and especially about combining the factor analysis for an ordinal scale (categorical data) - CATPCA with conventional PCA. Basically, I need to enter my set of variables in CATPCA, take the transformed variables from CATPCA and introduce them in conventional PCA (because, for example, I need to rotate the factor to obtain a reasonable solution). My concerns are:   1. After introducing the transformed variables in conventional PCA, is it okay just to follow the steps for conventional PCA and not taking in consideration that the initial variables have been transformed (standardized, from what I have observed) in CATPCA? I mean, can I go through evaluating MSA, communalities, factor loadings and so on, as in a regular case?   2. After obtaining the factors, if I want to assess the reliability, is it okay to take into consideration the initial variables (not the transformed ones in CATPCA which I also use in conventional PCA), assuming that I would like to generate the summated scales (for each factor) based on the initial variables and further use the summated scale for different constructs for correlation and regression?   3. Is it okay to calculate summated scales if I am undertaking an exploratory study and some of the factors that result are somewhat different from what I have found in the literature? (I have read that summated scales are not of much use if the solution obtained from the PCA is not supported by literature, for example if I have one or two factors that combine some variables in a different form than from I have found in the literature, thus I also name them differently)   4. Does it have importance in conventional PCA which of the criteria: communality and loading problems is taken into consideration for deleting variables in the conventional PCA? (should I first delete the variables that have a low communality or the variables that have cross-loadings?)   5. If I have a sample of over 350 respondents, is it acceptable to consider loadings over 0.50 as being practically significant or should I consider loadings over 0.30 as being significant (which would lead, for example, in the case of a variable that loads with 0.60 on a factor and respectively 0.30 on another factor to be considered as having a cross-loading)?   6. If I have a certain number of constructs and I would like to use the same methodology for each construct, is it okay to use CATPCA + conventional PCA even for the constructs for which I do not need to rotate the solution of the ones that are unidimensional? (I want to use a combination of CATPCA+PCA for a number of reasons, including to rotate the solution, but some of the constructs I am analyzing are unidimensional or, even I they are multidimensional, the solution may be okay without rotation)   7. Is it accepted to replace missing values with the mode in CATPCA while replacing the missing values with the mean (for the initial variables, not the transformed ones) for summated scales, correlation and regression? (This would also imply that I would use the initial variables and not the transformed ones for calculating summated scales, for correlation and for regression)   8. For summated scales: can I use all the variables that load on a factor or is it necessary to use only the variables with high loading (over 0.70, for example)?   9. Can I calculate a latent variable (construct) after the formula: Sum of (summated scale for factor i * variance explained by factor i), i taking values from 1 to n (n being the number of factors)? Thank you!


title: How to re-construct a matrix from SVD

text: I have a Audio time-series, to which I'm trying to detect the most significant parts of the signal, i.e. the voiced parts and forget the unvoiced parts. $$ T = [0, 0, 1, 1, .....n] $$ I then compute the Spectral result of the time series, which, gives me a dimensional frequency representation of the input signal. From this, I would like the determine the most significant elements of this, thus only leaving me with the voiced parts. I know this is possible using Principle component analysis and Independent component analysis, however, I'm trying to compute the following using SVD. To try and understand the theory and application, I constructed a simple matrix: $$ S = \begin{bmatrix} 1&0 \\\ 0&1 \\\ 1&1 \\\ 0&2 \\\ 1&1 \end{bmatrix} $$ In Python I have the following:               u, s, v = linalg.svd(S);      This therefore gives me three matrices, which, I understand that the following: "s" -> this matrix is a diagonal matrix which is ordered from singular value to lowest In order to re-construct the matrix, pick out the largest values of "s" and their corresponding columns in "u" and "v" Thus meaning that the "s" corresponds to the magnitude which relates to the columns in "u" and "v". But, I don't know how such can be re-constructed (using Python) and how it actually works..? My question mainly relates to taking a 2-dimensional vector and transforming this into another 2-Dimentsional vector containing only the most significant values. Any help would be greatly appreciated


title: Discrete Time Markov Chain - Inventory

text: Let $D_n$ be the demand for an item at a store on day $n$. Suppose that $D_n$ is a sequence of independent and identically distributed (i.i.d.) random variables with probability mass function: $p_k = P(D_n = k)$ Suppose the store follows the following inventory management policy. If inventory at the end of the $n$th day (i.e., after satisfying the demands during the day) is above $s$ the store does nothing. If it is less than $s$, then the store orders enough amount to bring inventory at the beginning of the next day up to $S$. Assume the inventory delivery is instantaneous. Let $X_n$ be the amount of inventory in the store at the beginning of the $n$th day. Model $X_n$ as a DTMC. Not sure how to deal with the probability mass function. I guess it would mean its equal but not sure how to put that in any other way rather than $p_k$. Nor do I understand how to deal with S.


title: Principal component analysis, psychometrics, and number of items required

text: I saw that questions regarding sample sizes and PCA are asked here quite frequently. However, I was not able to find exactly the information I need. I plan to do a psychometric study and conduct a principal component analysis on the data. The data I want to collect will be three-dimensional, e.g. participants × ratings × items. In many studies in the field, the data are aggregated in such a way that averages are taken across one of the dimensions (usually participants) so that in the end a ratings × items matrix remains. On this matrix, a two-way PCA is conducted, usually finding about two or three components that sufficiently describe the data. If I have, let's say, 50 to 100 participants and 10 or 11 ratings, how many items do I need to include? The components I want to find are supposed to summarize the ratings from the (aggregated) ratings × items matrix. Currently, I have seven items that the participants would rate on the 10 or 11 rating scales, the matrix averaged across participants would thus be 11 × 7\. Is this sufficient, or do I need to include more items in order for the PCA on the (aggregated) ratings × items matrix to be reliable? Please point me to the relevant threads here or elsewhere where I could read up more on the issue. Edit: If the low number of items (currently seven) is a problem, would this be alleviated if I calculated the correlations between the ratings for each individual first (so for each participant, I would have a ratings × items matrix and would calculate the correlation matrix of the ratings), and then averaged the correlations of the ratings over all participants? This technique was used in a different paper, and the authors said this would increase their sample size from `number_of_items` to `number_of_participants`. Is this true, and would this help me? Thanks again!


title: How to add two dependent random variables?

text: I know, I can't use convolution. I have two random variables A and B and they're dependent. I need Distributive function of A+B


title: Model Estimation

text: How can I estimate a model where Y ~ X1 + X2 + X3 X1 ~ Z1 + Z2 X2 ~ Z1 + Z2 X3 ~ Z1 + Z2 where Y may/may not be correlated with Z1 and Z2. Is there a R procedure that I could use to estimate the model. Thanks a lot for your help.


title: How do sampling error, measurement error and specification error affect regression coefficients

text: Can anyone answer the question or direct me to the proper resources - ESPECIALLY for sampling error effect on the coefficients? How would the coefficients be affected if the independent variables came from one data set while the dependent variable came from a different data set?


title: Find the amount of variation due to another covariate

text: I'm trying to explain a binary outcome (cardiovascular disease) with a categorical predictor (gene level, coded as 0, 1, or 2 depending on the number of risk alleles present). I'm trying to determine if the effect of the gene on cardiovascular disease is due to its known effect on an additional continuous covariate (blood glucose level), or whether it has a separate mechanism of action on cardiovascular disease. I'm not trying to find that other method of action, just to see if the effect of the gene on cardiovascular disease is due to its known effect on blood glucose, or whether it is independent of this. How, in your opinions, would one go about this? I'm trying to find some new tests that I haven't heard of before. Here's what I've done/thought of:   1. I've included an interaction term in the model to test whether the the expression of the gene is modified by the covariate. I got a significant P for the interaction term, so I know that the effect of the gene on the response is modified by the covariate.   2. I've looked at the Nagelkerke's pseudo R2 for a model with only the gene, only the covariate, and both of them together. I've found that the R2 of the combined models is less than the sum of the individual models (only by about 0.9%). From this I thought that since the gene explained ~ 2% of variation alone, but only added an additional ~1.1% of variation when taked onto the covariate only model, that part of the explained variance of the gene is actually due to its effect on the covariate. i.e. R2Gene only = ~0.02; R2Covariate only = ~0.2; R2Gene and coariate = ~0.211. So alone it explained ~ 2% variation, but when the covariate was included in the model, it only added 1.1% of variation. Is it legitimate to say that some portion of the variation explained by the gene is actually due to the variation in covariate? So my question is: a) are these things that I've done kosher, and b) how would I go about finding whether the variation explained by the gene is actually explained by the gene's effect on a covariate. Thank you for your time, any and all opinions are valuable and helpful.


title: Rare event logistic regression

text: I have a situation in which that there are 4 cases with a condition per 1000 people surveyed (binary). I am interested in performing a logistic regression to determine the risk factors (binary). Is it possible to perform such an analysis with SPSS?


title: Obtaining a representative sample when there is no sampling frame?

text: if there is no plausible sampling frame available for surveying, how will one obtain a representative sample of size n. In the textbooks they use SRS, stratified, cluster, and systematic but these all imply a sampling frame is available like a telephone book etc. So say for example, there is a website selling old cars, but they don't allow you to list them all, but you can search by make or model. How would you in this case get a representative sample of size n? Or if we were trying to get a representative sample of size n to survey the number of library books that are borrowed but we don't have access to the entire catalogue? But representative sample in the sense we can use it to make inferences about the population.


title: What is the difference between different kind of cross validation methods?

text: Now I am using MATLAB, and found one of its function named `crossvalind`, which provides lots of methods for cross validation, including some that I wasn't able to find on this site.               [TRAIN,TEST] = crossvalind('LeaveMOut',N,M), where M is an integer,     returns logical index vectors for cross-validation of N observations by     randomly selecting M of the observations to hold out for the evaluation     set. M defaults to 1 when omitted. Using LeaveMOut cross-validation     within a loop does not guarantee disjointed evaluation sets. Use K-fold     instead.      It seems that this kind of cross-validation can be used in a loop, but it cannot guarantee disjointed evaluation. I think that this kind of cross- validation has the advantage that its partition could be updated in every loop, rather than being fixed before the loop as in k-fold cross-validation. This means that more 'dynamics' could be tested by this kind of cross validation. But I am not sure about this leave-m-out cross-validation: Could anyone states the difference between leave-m-out cross validation and k-fold validation, what are their advantages or disadvantages, and when should we choose one rather than the other? * * * **To@a.desantos** In k-fold cross-validation, from the set of loops one sample can be in the test set for only one time, but in leave-M-out cross-validation one sample could appear in more than one test sets. **Edit 1** Some new issues inspired by @a.desantos. What m and n should I consider? Considering (by convention) that 10-fold cross-validation is nice, I should choose m as 4 when n would be 40; the combination of C(40,4) is 91,390. Should I generate all 91,390 different combinations, then run 91,390 tests, and average the performance of the specific learning algorithm? **To@Frank Harrell** **Edit 2** you mean the train-test routine will be runned for 10 times since each time the data is split into ten disjoint parts in a 10 fold cross validation, this will generate one average test error E1, and this error is not reliable, then i need to run this 10-loop-test-routine for about 5-10 times and average all E1s to get another test error named E2, this E2 stands for the ultimate average error, and is more reliable than E1, right?


title: non-normal data for two-way ANOVA, which transformation to choose?

text: I need to perform a two-way ANOVA on my data. My data is from a non-normal population. Apparently there is no two or three factor test for non-normal populations. I realized I need to transform my data, but I'm unsure about which transformation to perform on my data, I don't know which is the most appropriate. I don't know what is the criteria to choose one from the transformation list of possibilities?


title: What is a good resource on table design?

text: I've seen various theoretical treatments of graphics, such as the Grammar of Graphics. But I have seen nothing equivalent with regards to tables. Over the while I have developed an informal model of good practice in table design. However, I'd like to be able to provide a good reference to students. The APA Style Manual has a few tips on table design, but it is only a starting point. Question: What is a good resource that provides theoretical and practical advice on the presentation of numeric results in tables? **UPDATE:** It would be particularly useful to have a good free online resource. Note: I'm not sure if this should be community wiki. I feel as if there might be a correct answer.


title: Univariate Normal Converted to Multivariate Normal: Covariance Derivation

text: I am reading the paper available at this link: https://drive.google.com/file/d/0B2_rKFnvrjMARnU1QjB4anR3RDA/edit?usp=sharing I am having trouble understanding section 5.1 (page 2741). Essentially it says the following: $\theta_{ABi} \sim \mathrm{N}(\mu_{AB}, \tau^2)$ $\theta_{ACi} \sim \mathrm{N}(\mu_{AC}, \tau^2)$ $\theta_{BCi} \sim \mathrm{N}(\mu_{BC}, \tau^2)$ $\mu_{BC} = \mu_{AC}-\mu_{AB}$ implies $ \begin{pmatrix}\theta_{ABi} \\\ \theta_{ACi}\end{pmatrix} \sim \mathrm{N} \left(\begin{pmatrix}\mu_{ABi} \\\ \mu_{ACi}\end{pmatrix}, \begin{pmatrix}\tau^2 & \tau^2/2 \\\ \tau^2 /2 & \tau^2\end{pmatrix} \right)$ I do not understand how $\mathrm{Cov} \left[ \theta_{ABi}, \theta_{ACi} \right] = \tau^2 /2$ ? Could someone please explain this?


title: Significant lineal model (p value) but no significant parameters

text: I´m interested to know if it´s possible (and what it means) to obtain a significant linear model regression (p value=0.008) with coefficients fitted are not significant (p>0.05)?


title: Whether to use linear regression or not

text: Effects of growth hormone (GH) replacement with recombinant human GH on bone and mineral metabolism were studied in 36 GH-deficient children. Several outcomes, including serum ionized calcium levels, were assessed at pretherapy (zero months) and 1, 3, 6, 9 and 12 months after the beginning of therapy. Simple linear regression was used to determine if changes in ionized calcium levels could be predicted from length of therapy. Each patient’s ionized calcium level was employed in the analysis. **Is the use of simple linear regression appropriate?** **a.** Yes **b.** No **c.** The investigators should have used logistic regression instead


title: How to actually plot a sample tree from randomForest::getTree()?

text: Anyone got library or code suggestions on how to actually **plot a couple of sample trees** from:               getTree(rfobj, k, labelVar=TRUE)      (Yes I know you're not supposed to do this operationally, RF is a blackbox, etc etc. I want to visually sanity-check a tree to see if any variables are behaving counterintuitively, need tweaking/combining/discretization/transformation, check how well my encoded factors are working, etc.) (Prior question without decent answer: 'How to make Random Forests more interpretable?'. Also 'Obtaining knowledge from a random forest' I actually want to **plot a sample tree**. So don't argue with me about that, already. I'm not asking about varImpPlot(Variable Importance Plot) or partialPlot or MDSPlot, or these other plots, I already have those, but they're not a substitute for seeing a sample tree. Yes I can visually inspect the output of `getTree(...,labelVar=TRUE)`) (I guess a `plot.rf.tree()` contribution would be very-well-received.)


title: How does confidence interval vary with sample size?

text: I am trying to convert state-level poll results for 2012 U.S. presidential election into state-level probabilities of each candidate winning. To achieve that, I am using the logistic function, which conveniently goes to 0 or 1 as poll results deviate from 50%, and also makes the probability of each candidate winning a state equal to the probability of the other candidate losing that state. As a reminder, the logistic function is: $ P(x) = \frac{1}{e^{-x} + 1} $ (NOTE: To match the domains of poll results ($p \in [0, 1]$) and logistic function ($x \in [-\infty, \infty]$), I use the tangent function: $x = tan(10*p/\pi - \pi/2)$, which also conveniently sends $p=50\%$ to $x=0$, the center of the logistic function. I'm not kidding! But don't worry about this unless you have a better suggestion off the top of your head.) My question is how to incorporate "sample size" into the logistic function, $P(x)$, so that its sharpness around center reflects what percentage of total population (or likely voters) has been polled. I am thinking that if I could write the confidence interval as a function of sample size, I could divide $x$ by the length of that interval and thus adjust the sharpness of the logistic function.


title: Timeseries regression

text: I'm following an undergraduate course on timeseries using OxMetrics and wanted to reproduce som results in R Estimating an ARMA(3,3) model:               arima(temp,c(3,0,3))      using stats package and also TSA package results in following error:               Error in arima(temp, c(3, 0, 3)) : non-stationary AR part from CSS      Tried to use ML-method:               arima(temp,c(3,0,3),method="ML")      no error so I get estimates, but not the estimates found in OxMetrics. Then I tried installing the package FitARMA and used:               FitARMA(temp,c(3,0,3),TRUE,FALSE)      Which results in the same estimates for the coefficients as when estimated in OxMetrics using exact MaxLikelyhood except for the constantterm. Then I tried changing fourth argument of FitARMA from FALSE to TRUE:               FitARMA(temp,c(3,0,3),TRUE,TRUE)      Estimating mean by max.likelyhood and then getting the same estimate for the constantterm but now the estimates for the coefficients change to approximately the same as when using arima in R. I've already consulted four books on R and I can't really figure out why this is happening. Both OxMetrics and R use BFGS for one. Being quite new to timeseries I'm quite lost...


title: Compare Regression Coefficients

text: I want to test if two regression coefficients from two separate regressions are significantly different from each other. In particular I want to test from the model $$ \begin{aligned} x1_{(t)} &= \mu1 + a*x1_{(t-1)} + b*x2_{(t-1)} + e1 \\\ x2_{(t)} &= \mu2 + c*x1_{(t-1)} + d*x2_{(t-1)} + e2 \end{aligned} $$ if $H0: b=c$ against $H1: b>c$. Can I do this with a simple t-test or do I require another test statistic? * * * Just wanted to give a brief update. In my reference paper, they claim that they can test H0:b=c against H1:b>c in the model stated above with a simple Z-test. Does this mean that my test statistic is simply (b-c)/se(b)?


title: Changepoint Analysis in R for Time Series Data

text: I'm new to R and statistics in general. My background is in Software Development, and I've tasked myself to solve a problem in a 3rd party application. This application is a Game, and I'd like to measure the existence of a performance problem. I measure performance in Frames Per Second (FPS), and for some reason under high load the system exhibits periodic drops (hesitation) in performance. I want to measure the frequency, severity of this drop in performance. The idea is to measure whether changes in the codebase/testing-conditions has changed the base & hesitation performance. I acknowledge this is similar to time-series analysis of data like stocks, temperature, etc. However the resources out there that I would like to apply hasnt quite sunk in yet. So I'm looking for set a responses for me to review and hopefully teach myself. I've been able to load data into R Studio and create graphs. I've followed a few other answers to create some smoothed lines, but it didnt quite work for my dataset. http://imgur.com/a/wnLQs \- An album of the data graphs. Any help would be appreciated.


title: What types of data analysis do not count as statistics?

text: When does data analysis cease to be statistics ? Are the following examples all applications of statistics ?: computer vision, face recognition, compressed sensing, lossy data compression, signal processing.


title: Comparing the result of a study which has unequal group sizes

text: I have conducted two user studies and in my studies I didn't have control over the group sizes. In each study users were put in groups and they were asked to perform some group activities. Here is the way the groups were arranged:   * **Study 1 (Control Groups):** In this study 20 people participated. The participants were in 2 groups of 4 people, and 4 groups of 3 people.   * **Study 2 (Modified Groups):** In this study 25 people participated. The participants were in 5 groups of 4 people, and 3 groups of 3 people. I have gathered two kinds of data during my study:   * Individual data. For example the number of ideas that each person put forward.   * Group related data. For example the time group spent to do something. I am looking for a method which helps me combine the results for each study and then compare the results. I was wondering if anyone could help me with this ? More information about the task and data: The task is Brainstorming ideas about a specific topic. The topic is the same between all of the groups. Users work in a software that I wrote. Test groups see a different version of the software. I am measuring the time it took the group to generate 10 ideas, and I want to see if software version had an effect on it or not. But the main problem is having groups of 3 and groups of 4 users.


title: Maximum power test

text: **Question:** In a count that was held a year ago in a reserve , there were 15 tigers counted. 7 females and 8 males, which were marked with numbers from 1 to 15. a. By recent reports, there was fear that the number of tigers in the reserve had dropped under 15. In order to check that, it was decided upon going out and recording the numbers of 4 tigers in the reserve, which will appear randomly. If the maximum number of these 4 records is k or lesser, then the conclusion is that the number of tigers really decreased.   1. if the count was done with return. what is the value of K, if the required p-value is 0.1?   2. If the count was done without return , calculate the p value of the test when k=9. b. In a random sample of 4 tigers, without return , there was 1 male and 3 females found. Now we need to check if the first count (a year ago) had a mistake, and actually there are 15 tigers in the reserve but among them 8 females and only 7 males. 1\. find a maximum power test based on X - the number of females in a 4 tiger sample. 2\. find the p-value of that test, and make a conclusion for a=0.1 3\. find the probability of a type II error, and the power of the test. **My thoughts** I did part a, and wasn't sure about it. about part b- didn't really know how to do it- is doing the ratio of likelihood functions a good way to start? I don't know what is the parameter for that though... in A.1. i said that $X_i\sim U(1,15), S(X)=max(X_i)$. $P(S(x)\le k)=0.1$ so $0.1=\frac {k-1+1}{15-1+1}=\frac k{15}$ so $k=1.5$ -> $k=1 $(k has to be a whole number) A.2. $P(S(x)\le9)=P(X_1\le9)P(X_2\le8)P(X_3\le7)P(X_4\le6)=\frac 6{65}$


title: Aside from regression coefficients, what are commonly used approaches to measure one variable's "sensitivity" to another variable?

text: I am about to embark on a lengthy study for a work project. At the core of it is the need for a statistical tool to help quantify "sensitivity." I am the lead software developer on the task, but I have a pretty strong background in statistics (at least compared with many other software developers). However, my stats background never touched on this idea of sensitivity, especially not outside of an OLS approach. I have been reading the Wikipedia article on Sensitivity Analysis today, and believe the section on regression analysis touches closely to the problem I am having. Let me explain: The approach that management wants to take is more or less precisely the "garbage-can" approach from Achen's paper "Let's put garbage-can regressions and garbage-can probits where they belong.". The goal is to make a nice software front-end that gives our researchers a simple interface to select data elements from a bunch of database tables, apply time series lags and/or simple pre-treatments and data cleaning, and then "dump" them in to a regression, with all of the coefficients and summary stats similarly "dumped" out to a big table o' coefficients at the end. I have been diligent in making this tool and I'm pretty proud of how easy it is to use. Our parameter scripts require very little to be written beyond a simple string such as "DepVar ~ alpha + beta_1*Var1 + beta_2*Var2 + ..." and some simple definitions about what "DepVar" or "Var1", etc., means in the string. But I also can't help feeling pretty bad that "productionizing" this kind of thing is _super_ misleading. The system is entirely built on ordinary least squares, with no ability to tweak assumptions about heteroskedasticity, error distributions, etc. You cannot swap in other regression techniques for OLS. It's all OLS, all the time, and this was a conscious choice by management. The output of the system will be a very large set of time-series coefficients on an individual-by-individual basis. So, in down-stream parts of our analysis pipeline, other programs can load a big time series of, say, "Individual X" past sensitivities to some variable. These are used for risk modeling and outlier treatments among other things. I can't help but feel like this is a very bad way to quantify the sensitivities. Any time there are violations of the OLS assumptions, and especially when these relationships are non-linear, the coefficients can be extremely misleading, even if they have good-looking Frequentist summary stats. I've raised this issue with my team many times, but the managers in charge seem extremely hell-bent on having simple-to-manipulate statistics. They also seem overly concerned about "interpretability" -- in the sense that they believe these coefficients from OLS can be directly read off as normalized loadings of one variable on another variable. I agree that can sometimes be true, but only if you're really sure that the linear model isn't badly misspecified. And having looked at the data myself, I do believe the linear model is often badly misspecified. I have quantified that and repeatedly shown this well-documented fact, along with the Achen paper, to my managers, but they just will not budge. Lately I have been wondering: how do other approaches tackle this problem? In a machine learning context, like fitting a SVM model or a decision tree, how do practitioners quantify the sensitivity between the different variables? How does this work in a setting where there may be non-linearities?


title: Measure of dispersion for dataset with zero mean

text: I'm testing a machine learning algorithm which, until now used coefficient of variation in order to relative measure dispersion among different attributes in a given dataset. Now that I need to test it thoroughly on larger dataset, it's critical to do something about the possibility of a zero mean, which would immediately breaks down coefficient of variation (cv), and hence the algorithm. Please note that the distributions are not essentially Gaussian, and are often skewed when they are. Also, the measure of dispersion needs to be dimensionless. I searched the web a bit to find ways to 'normalize' cv in some manner, or to even find an alternative to cv. I tried inter-quartile range(iqr), but there was a significant drop in performance when I migrate from cv to iqr. Is there any other dispersion measure that returns similar values to that of cv? If not, is there any principled way to 'normalize' coefficient of variance?


title: Linear Mixed Effects Models

text: I have commonly heard that LME models are more sound in the analysis of accuracy data (i.e., in psychology experiments), in that they can work with binomial and other non-normal distributions that traditional approaches (e.g., ANOVA) can't. What is the mathematical basis of LME models that allow them to incorporate these other distributions, and what are some not-overly-technical papers describing this?


title: Distribution with upper and lower boundaries

text: Is there any continuous distribution that has an upper level of 1 and a lower level of 0? The reason I am asking this is that I want to know the distribution of my data which shows the roll over ratio of customer deposit which might only be between 0 (means all the fund is withdrawn) and 1 (means all the fund is rolled over).


title: What to do when a Mann-Whitney U assumption is violated?

text: I'm using a M-W U test to analyse some Likert scale results I have, as my data is ordinal. One of the assumptions that keeps on appearing on references is the following: "both distributions must be the same shape (i.e., the distribution of scores for both categories of the independent variable must have the same shape). If I fail this assumption, and therefore am unable to use M-W U, what can I use instead? EDIT: I cannot add images, but will try to explain what my data looks like. I have two groups based on my IV ("Yes" and "No"), charted on a boxplot with y-axis ranging from 1 - 5. For "Yes", I have a box extending from 4 to 5, with no variance. For "No", the box ranges from 3.5 to 4.5, with variance from 2.0 to 3.5, and 4.5 to 5.0. I assume this means my shapes are "different". Does this mean I can't use MWU? If not, can you make suggestions as to what I can use instead (SPSS)? EDIT2: Sample data:               "Yes": 5 5 5 5 5 5 4 4 4 4 5 5 4 4     "No":  4 5 5 5 3 2 4 3 4 4 4 4      I was using MWU to detect any significant differences (apparently, "yes") - but now am not sure if that's the right way to do it.


title: Equality in linear regression

text: I am looking for some hints to prove the following equality: $$ y^{\top}y - y^{\top}X(X^{\top}X)^{-1}X^{\top}y = \dfrac{\det(L^{\top}L)}{\det(X^{\top}X)}, $$ where $y$ is a $n\times 1$ vector, $X$ is a $n\times m$ matrix and $L=(X,y)$ (concatenation with $y$ as the last column of $X$). **Upddate:** Thank you all, I (cross-)posted the question in math stackexchange, where it received more attention (at least the attention I needed). Feel free to close this question.


title: Is there something like an absolute scale type?

text: Some measures seem to be of a different type to others, depending on what kind of statements are meaningful. Different scale types try to capture that difference. Known interval scale types are:   1. Nominal   2. Ordinal   3. Interval   4. Ratio The great thing is, that each scale type builds on the previous one, adding additional properties. I came across a fifth scale type, which comes after the ratio scale type. The problem is, I only found it in the book _Software Metrics: A Rigorous & Practical Approach_. also explained on these slides. > Absolute scale measurement is just counting > >   * The attribute must always be of the form of ‘number of occurrences of x > in the entity’ >   * number of failures observed during integration testing >   * number of students in this class >   * Only one possible measurement mapping (the actual count) >   * All arithmetic is meaningful > A previous definition was, that the absolute scale was simply on a fixed range of values, but the more I seek for it I come across this more simple definition: it's limited to counting. Any thoughts on this? I feel a bit insecure here.


title: Can a variable be a fixed effect and random effect variable at the same time?

text: I am building a mixed model, where the same 4 groups are tested repeatedly 5 weeks on the measure M. I want to test the effect of group on the measure M. Can I build a model as following (let's say with R):                  lmer(M ~ Time + Group + (1|Group), data= mydata)      where group is both a fixed and random effect at the same time? Imposing two effects of a group seems contradictory conceptually. But what if I want to test the group effect?


title: Visual explanation of forecasting algorithm

text: I am trying to develop the concept to give an intuitive explanation of how my forecasting algorithm works. The forecasting algorithm uses historic data and some current-day data to forecast the remaining work for the day as it is a workforce planning algorithm. It does not contain any big datasets and only a few input parameters such as number of incoming documents, number of workers etc. The forecasting is just for one day. I am not sure what visualization technique I can use in such a scenario. I have tried to use a line graph with mouse-over but the users are not able to understand the shape of the graph. They want to understand how the rules of the forecasting algorithm affect the graph shape. My question is: what visualization technique can I use to make it more intuitive, so that I can explain how different rules affect the graph shape and also make them interact with the technique?


title: Spearmans and variable units

text: Both my variables for correlation (Spearmans given small n of 17) relate to time (seconds for x, and minutes and secs for y). Should I standardize the units before entering into analysis. I get different answers if I try both ways, which is correct?


title: How to examine interactions between factor and covariate in a mixed effects model?

text: I have 2 factors `A` and `B` (5×3) and one covariate `X` in a within subject design. Here's how I specify my overall model:               lme.out = lme(y~ A*B*X, random=~1|Subject, data=mydata)      My interpretation is that I am looking at a graph `y~x`, where the slope changes due to the covariate, and the line shifts up or down based on the different levels of `A` and `B` (changes in intercepts). What I want to find out is: if I were to fix factor `A` (take any of the levels), then looking at the lines (`y~x`), what is the effect of `B`? Does the levels of `B` shift the line up or down (intercepts) or does it alter the slope of the line (`X`). Should I be running some sort of contrasts analysis? But I am not sure how contrasts work between factors and covariates. One way I could think of is to take the subsets of data corresponding to different levels of A and create models such as: `lme(y~ B+X, random=~1|Subject, data=mydata[which(mydata$A = A1,])`. This way I could compare the resulting intercepts and slopes across these models. Can anyone tell me if what I am doing make sense? Suggestions of any kind would be greatly appreciated!


title: Detect outliers in mixture of Gaussians

text: I have a ton of univariate samples ($x_i \in \mathbb{R}^+$). I'd like an automated method to check for outliers and identify the outliers, if any are present. A reasonable model for the distribution of the non-outliers is a mixture of Gaussians. The number of Gaussians in the mixture and their parameters are not known a priori. Can you suggest a simple method for identifying outliers? Do you have any recommendations? It'd be nice if it were simple to code up in Python. Something quick and dirty -- say, easy to understand, easy to implement, and pretty effective-- beats something complex but optimal. For example, I'm a bit reluctant to wade into something fancy based upon expectation maximization. Example parameters: I might have 10,000 samples or so. The distribution of non-outliers might be a mixture of 2 Gaussians; or I might have a mixture of a few hundred Gaussians. **Update:** People have asked how anything could possibly be an outlier, given these assumptions. (Presumably, the unstated concern is that this problem may be unsolvable: if every data set is always explainable by some mixture model, then there's no basis to ever identify anything as an outlier.) That's a fair question, so let me try to respond. In my application domain, I can reasonably assume that there will be dozens of samples from each component Gaussian. e.g., I might have 40,000 samples from a mixture of 100 Gaussians, where each Gaussian component has a probability no lower than 0.001 (so it is almost guaranteed that I have at least 10 samples from each Gaussian). I realize I didn't state this assumption earlier, and I apologize for that. However, with this additional assumption, I believe the problem is solvable. There exist examples of data sets where one or more points can be considered outliers (they cannot reasonably be explained by any mixture model). For example, consider a data set that has a single isolated point that is very far from all others: if it's far enough away, it can't be explained by the Gaussian mixture model and thus can be recognized as an outlier. In conclusion, I believe that the problem is well-defined and is solvable (given the additional assumption stated here): there do exist example situations where some points can reasonably be identified as outliers. Note that I'm not trying to propose a special or unusual definition of outlier. I am happy to use the standard notion of outlier (e.g., a point that cannot reasonably be explained as having been generated by the hypothesized process, because it is too unlikely to have been generated by that process).


title: Determining largest contributor in a group

text: I dont know much about statistics so bear with me. Lets say I have a set of 1000 workers. I want to figure out who the hardest worker is but I can only measure the amount of work getting done in groups of 1-100 over an hours worth of work. Assuming each worker always does around the same amount of work, over a large number of trials and combinations can I rank my workers by who works hardest? note: this is just a metaphor so dont worry about actually running the tests, just assume I already have a large set of data. **Edit:** when I say "Assuming each worker always does around the same amount of work" I mean each individual does around the same amount of work on a day to day basis. So joey will do around 100 units of work each day and greg will do around 50. The problem is I can only observe the number of units of work done by the group. **More edits:** In regards to the number of workers working at once and the frequency of them working. There could be any number of workers working at the same time. Some workers will probably end up working a lot more than others, that is to say, we can assume some workers will be working almost 90% of the time and others almost never. I know that makes it difficult but I will have a very large dataset so hopefully that makes it a little bit easier. **For each hour we know which workers are working and how much work got done. From that information I want to find out who is doing the most work.** if the data were in JSON format it would look something like this...               [       {         "work_done": 12345,         "Workers": [ "andy", "bob", "cameron", "david" ]       },       {         "work_done": 432,         "Workers": [ "steve", "joe", "andy"]       },       {         "work_done": 59042,         "Workers": [ "bob", "aaron", "michelle", "scott", "henry" ]       },       ...     ]


title: Linear model with log-transformed response vs. generalized linear model with log link

text: In this paper titled "CHOOSING AMONG GENERALIZED LINEAR MODELS APPLIED TO MEDICAL DATA" the authors write: > In a generalized linear model, the mean is transformed, by the link > function, instead of transforming the response itself. The two methods of > transformation can lead to quite different results; for example, **the mean > of log-transformed responses is not the same as the logarithm of the mean > response**. In general, the former cannot easily be transformed to a mean > response. Thus, transforming the mean often allows the results to be more > easily interpreted, especially in that mean parameters remain on the same > scale as the measured responses. It appears they advise the fitting of a generalized linear model (GLM) with log link instead of a linear model (LM) with log-transformed response. I do not grasp the advantages of this approach, and it seems quite unusual to me. My response variable looks log-normally distributed. I get similar results in terms of the coefficients and their standard errors with either approach. Still I wonder: If a variable has a log-normal distribution, isn't **the mean of the log-transformed variable** preferable over **the log of the mean untransformed variable** , as the mean is the natural summary of a normal distribution, and the log-transformed variable is normally distributed, whereas the variable itself is not?


title: Standard deviation of the sum of regression coefficients

text: I'm doing OLS estimation with an independent variable lagged as t-1, t-2, t-3, and t-4 (four beta coefficients). I would like to have the sum of these coefficients for interpreting the net impact of these variables on the endogenous one (for example the net influence of the US's GDP growth lagged in t-1, t-2, t-3, and t-4 on the Thai GDP growth). How can I compute the standard deviation of this new coefficient for computing a t-test? I have found this formula but i'm not sure... $$ SE_{b_{new}} = \sqrt{SE_1^2 + SE_2^2+2Cov(b_1,b_2)} $$ (see this link Adding coefficients to obtain interaction effects - what to do with SEs?) But how do I calculate the covariance of two (four in my case) coefficients? Does that make any sense? Thank you very much. * * * Ok, I need the covariance matrix of coefficient to retrieve the cov(b1,b2) cov(b1,b3)... After this, i should be able to compute the standard errors of the sum of coefficients with the formula above...


title: Random walk with drift are differences white noise?

text: If I have a random walk without drift the differences form a white noise process. But what happens if I incorporate a drift $d$? Does this still hold true? I'm not sure because with the drift term d: $$ y_{t} = d + y_{t-1} + \epsilon_{t} $$ $$ y_{t} - y_{t-1} = d + \epsilon_{t} $$ The right hand side of equation 2 can't be a white noise process (for instance its expected value would be $d \neq 0$. Any help would be appreciated!


title: Assigning weights

text: A Company has 4000 employees. Employee1 works 3 hours on monday, 1 hour on tuesday, 0 hours on wednesday and thursday, 5 hours on friday and 6 hours on saturday. It seems he worked more hours on saturday. But a random employee works 0 hours on monday to thursday and 5 hours on friday and 6 hours on saturday. Therefor employee1 actually worked more hours on monday compared to a given random employee. Is there any way to statistically assign employee one worked more on Monday compared all other days ?


title: How to define initial probabilities for HMM?

text: HI This is first time I was reading about HMM, however I have read so many articles on web, but two things where I am confused are:   1. How to decide number of Hidden States (although HMM says we don't need to know, we just have to make a guess, even for making guess what should be the best criteria)   2. Once define hidden states let say 5, then how to define initial probabilities for each hidden state and the transitional probabilities among each other... I would appreciate if someone can give me example... But please don't give me example of weather system.


title: Rao-Blackwellization of sequential Monte Carlo filters

text: In the seminal paper "Rao-Blackwellised Particle Filtering for Dynamic Bayesian Networks" by A. Doucet et. al. a sequential monte carlo filter (particle filter) is proposed, which makes use of a linear substructure $x^L_k$ in a markov process $x_k = (x^L_k, x^N_k)$. By marginalization of this linear structure, the filter can be split into two parts: a non-linear part which uses a particle filter, and one linear part which can be handled by a Kalman filter (conditioned on the non-linear part $x^N_k$). I understand the marginalization part (and sometimes the described filter is also called marginalized filter). My intuition why it is called a Rao- Blackwellized Particle Filter (RBPF) is that the Gaussian parameters are a sufficient statistics for the underlying linear process, and following from the theorem of Rao-Blackwell an estimator conditioned on these parameters performs at least as good as the sampling estimator. The Rao-Blackwell estimator is defined as $E(\delta(X)|T(X)) = \delta_1(X)$. In this context I would guess that $\delta(X)$ is the monte carlo estimator, $\delta_1(X)$ the RBPF, and $T(X)$ the gaussian parametrization. My problem is that I don't see where this is actually applied in the paper. So why is this called a Rao-Blackwellized Particle Filter, and where does the Rao-Blackwellization actually happen?


title: How to write log-likelihood for beta-geometric with optim() in R?

text: I am trying to estimate the parameters for a shifted beta-geometric distribution to model user churn, as shown in this paper. The log-likelihood function is described there and solved via Excel, and I am attempting to do the paramter estimation in R. The function is below. $$ LL(\alpha,\beta \space|\space data)=\Sigma \space n_t \ln[\frac{B(\alpha+1,\beta+t-1)}{B(\alpha,\beta)}]+(n-\Sigma \space n_t)\ln[\frac{B(\alpha,\beta+t)}{B(\alpha,\beta)}] $$ where $n_t$ is how many people were lost in a period, and $n$ is the number of people in the original cohort. I have $n_t$ in the vector `lost` and precalculate $ n -\Sigma n_t$ in `active[tmax]` My attempt to do so in R is as follows:               sbg.lik<-function(params,t,n,active,lost,tmax) {     a<-params[1]     b<-params[2]     logl<-sum(lost*log(beta(a+1,b+t-1)/beta(a,b)))+active[tmax]*log(beta(a,b+t)/beta(a,b))     return(-logl)     }          optim(c(1,1),sbg.lik,t=1:7,n=1000,active=active,lost=lost,tmax=7)      This is unfortunately not working, with the error `Error in optim(c(1, 1), sbg.lik, t = 1:7, n = 1000, active = active, lost = lost, : objective function in optim evaluates to length 7 not 1` Any ideas? Am i approaching the log likelihood maximation correctly? Any advice would be greatly appreciated!


title: Can this be modelled?

text: _Context:_ USA economy _Background:_ Its generally accepted that the growth of e commerce has certain curbing effects on the CPI-inflation. Because search costs are much lower online, people always go for the lowest prices and this makes online retailers fight to price their products lower than everyone else. Plus, there are not much fixed costs in online retail which helps reduce overall prices. Now with the rise of e- commerce's relevance and ease of use, people are switching to it for grocery shopping and electronics shopping, etc and therefore this is seen to have a curbing effect on inflation CPI, as prices are driven lower. _Consider the two processes:_ Let one process be the CPI inflation   Let another process be the growth of e commerce in US(maybe in terms of %of GDP or some other metric) Is there a way to do a time series analysis(or any other type of analysis if you feel time series is not the way to go) to test if this inflation curbing effect of e commerce exists? Will this test be meaningful? This is for a project so I'd like some input on whether its going to give me interesting results or if this analysis is too general to yield anything useful. Thank you!


title: Parametric modelling with accelerated failed-time models

text: I am a clinician who is trying to fit a model that predicts the absolute risk of individuals developing cardiovascular events at the end of follow up period. I came across accelerated failed time models that should enable me to achieve the above purpose. But I don't know how to interpret the output of the model in R. Can anyone help? How do I generate output that tells me about the absolute predicted risk of individuals??                     surfit <- survreg(Surv(time, delta) ~ as.factor(stage) + age, dist='Weibull')           summary(surfit)                                           Value Std.Error    z     p             (Intercept) 3.5288  0.9041     3.903 9.50e-05        as.factor(stage)2 -0.1477 0.4076    -0.362 7.17e-01       as.factor(stage)3 -0.5866 0.3199    -1.833 6.68e-02       as.factor(stage)4 -1.5441 0.3633    -4.251 2.13e-05                     age -0.0175 0.0128    -1.367 1.72e-01              Log(scale) -0.1223 0.1225    -0.999 3.18e-01


title: Mixed effect model validation

text: I fit a mixed effect model to one set of data and now I want to check its performance when applied to another data set. I want to plot a figure of "predicted value" and "observed value" and calculate $R^2$ to evaluate model performance. Is there any code or command for this in R software?


title: SciKit SGD Regressor RBF Kernel Approximation

text: I am using scikit-learn and would like to run SVR with RBF kernel. My dataset was quite large so from reading other posts, I was advised to use SGD regression and RBF approximation. Interestingly I get worse results from using SGD with RBF compared to just SGD by itself. I'd imagine this would be due to wrong parameter values. I tried looping through gamma and n_components for RBF sampler and tried a number of parameters for SGD regessor to no avail. I also output the training and cross validation error, and both are roughly the same, so I believe its a bias issue.


title: non normal subsample within normal sample: parametric or nonparametric?

text: Potentially simple question here but I can't work out the answer, so... I have created a single question questionnaire with five possible yes/no answers and then calculated the percentages for each answer (N=34) with the distribution being normal as per a Kolmogorov-Smirnov test. I then split this sample into two subsamples, male and female, with the aim of comparing these data with the population but whilst the female data is still normally distributed, for one of the five answers for the males the K-S statistic is now significant (0.003). Do I still calculate the means (parametric) for males and females and compare these with the population mean values or do I calculate medians for the (otherwise normally distributed) population and compare them with the medians for the males and females? Or is there an entirely different approach? All advice most appreciated!


title: How to analyse if there is a link between answers on Likert scale and gender/age?

text: I am trying to analyse the results of a questionnaire with SPSS. I have a number of questions with Likert scale and want to see if there is a link between the age and gender towards the answers given whether they strongly agree/agree/…. I tried Chi square (no idea how to analyse the outcome) but I was told that chi square does not work with nominal data. Is that true? What would be the best method?


title: Significance cutoff point for multiple testing

text: I have some results from a genome wide association study, which basically fits $n$ linear models for one trait $y$ and $n$ genotypes $x_1, x_2, \cdots, x_n$: \begin{align*} y &\sim b_{10} + b_{11} x_1 \\\ y &\sim b_{20} + b_{21} x_2 \\\ &\vdots \\\ y &\sim b_{n0} + b_{n1} x_n \\\ \end{align*} after which we get $n$ p values: $p_1, p_2, \cdots, p_n$. Normally $n$ is in the order of $10^6$, and the $\alpha$ level is usually set at $5 \times 10^{-8}$ in the literature, roughly following the Bonferroni correction. Now based on these p values we have spotted an interesting region for which we have built a series of different models and conducted around 10000 more tests, the problem is how do we determine the significance cutoff point this time? Someone has suggested to set it at $5 \times 10^{-8} / 10^4 = 5 \times 10^{-12}$, but apparently this is not quite reasonable, the 10000 tests are only a small number compared to the millions of test conducted previously.


title: What is the proper way to summarize class attendance and GPA?

text: A friend's assignment question is "Does class attendance affect GPA?" She has 3 survey questions:   * How many courses do you have a week?   * How many classes did you miss a week?   * What is your GPA? She is now told to analyse her data: her assignment question told her to analyse by combining two of the questions, but in her case, I told her she would need to use all three. I was thinking of doing (percentage of class missed, (class missed/class taken a w/) compared to GPA. She thinks she can just do classes missed compared to GPA. I told her this wouldn't work because someone could attend 3/3 of their classes, while someone else can attend 3/5 of their classes, and it wouldn't be well represented. Could someone give some insight on the best method to represent the data to graph?


title: Variance of a difference in marginal proportions in a three-way contingency table

text: Let a multivariate distribution be given by $P(Y,S_1,S_2)$, where all three variables are discrete, $Y$ is multivalued, $S_1=(0,1)$ and $S_2=(0,1)$, respectively, and all may be dependent. Define the quantity $X$ for difference in means on categry $j$ of $Y$ as $$X=E(Y=j|S_1=1)-E(Y=j|S_2=1).$$ To give some background to this statistic, assume $S_1$ and $S_2$ are response indicators that indicate whether $Y$ is observed or not observed, e.g. in two indepdent surveys. $S_1$ and $S_2$ may be dependent, because the probability of response in two different surveys should be assumed correlated within the same individual. $X$ then provides the population difference in $Y$ between response groups to both surveys. I estimate $X$ by $$\bar{x_1}=\frac{1}{\sum _{i=1}^{n}S_{1i}}\sum _{i=1}^{n} I(Y_i=j)S_{1i}-\frac{1}{\sum _{i=1}^{n}S_{2i}}\sum _{i=1}^{n} I(Y_i=j)S_{2i},$$ where $I$ is the indicator function and the sum is taken over a sample of size $n$, where $i=1,..,n$. This is the differene in marginal (sample) proportions in category $j$ of a three-way contingency table made up by $P(Y,S_1,S_2)$, let's say $p_{1j}$ and $p_{2j}$. The goal is to derive the standard error for the difference between both quantities. We have: $$\text{Var}(p_{1j}-p_{2j})=\text{Var}(p_{1j})+\text{Var}(p_{2j})-2\text{Cov}(p_{1j},p_{2j})$$ I know how to estimate the variance parts of the formula (i.e., $p_{1j}(1-p_{1j}/n_1$)), but I am clueless about the covariance. The covariance of a multinomial variable is $-n\pi_k\pi_l$ where $\pi$ is the population proportion and $k$ and $l$ two categories. But I do not think this formula applies here. Thank you for advice on where to move from here.


title: Selection bias correction and a multinomial logit

text: I have a data set for a number of people making 2 decisions - where to live; and how many hours to work. For every observation with a non-zero amount of work, there is an observed wage. I've assigned each person a unique 'choice' representing the combined observed outcomes. There are 8 work outcomes (0 hours, 8 hours ... 56 hours) and 5 location outcomes (Region A, B, C, D, E), giving a total of 5x8 (40) possible choices. I want to estimate BOTH the utility function driving the decision, and the wage for the people currently not working. So of the 40 possible choices, 5 have no observed wage. I can see that if I knew all the wages, I could use a mutlinomial logit model (or random parameters logit model) to estimate the utility function (using package mlogit). I also know that if I was just trying to correct the selection bias in the wage estimation, I could use the Heckman procedure, estimate a logit labour supply and then get the unbiased wage equation (using sampleSelection). I'm interested in combining both, and I was wondering if there was a existing package that enabled the combination of the two?


title: missing value patterns

text: I am doing some data preparation with Python using Pandas and I am working with a dataset that has about 80 variables with missing values and I want to capture any patterns of missingness to cut down on the amount of missing value indicators I have but I am having trouble finding any good strategies for doing this. here is an example of what I am getting at:               MISS_1 MISS_2 MISS_3 MISS_4 MVP1 MVP2 MVP3 MISS_STR          0      0      1      0    0    0    1     0010          1      0      1      0    1    0    0     1010          1      1      1      1    0    1    0     1111          1      1      1      1    0    1    0     1110          1      0      1      0    1    0    0     1010          0      0      1      0    0    0    1     0010          0      0      1      0    0    0    1     0010          0      0      1      0    0    0    1     0010          1      0      1      0    1    0    0     1010      One thing I tried was to create a string variable that concatenated all the missing value indicator variables like the 'MISS_STR', unfortunately the number of unique values in this variable was around 2000. Also, I know I am missing patterns that may be useful because if a pattern exists between Variable_1 and Variable_2 and Variable_3 and Variable_4 are missing completely at random, then concatenating them as strings will not capture the pattern between Variable_1 and Variable_2. Is there a better way of doing this? A good clustering algorithm perhaps? Whatever I do, I need to be able to do the same thing to my scoring data.


I have a quick question: How do you make a $10 \times 2$ $\chi^2$ table in SAS?


title: Clustering of multivariate data

text: Please I am about to cluster some data based which have about 15 different columns all of which are numbers(Some categorical while some are measurements) also some of my values are missing in some columns . Please can you give me pointer on how to go about it. I have previously explored the clustering with weka but I am not sure about the way weka implements so I am going the R route. What I know : I already know about Principal components analysis at least in theory. But is this necessary in all clustering of multiple columns . It will go a long way if anyone could provide me a link to a tutorial on this because Quick-R has for just 2 variables. A sample of my dataset is listed below               1,64,9,30,33,2,3,1,6,1,5,-3.62,-3.71,-2.73,1     2,61,4,30,33,2,3,2,7,4,4,-3.62,-3.71,-2.00,1     3,49,4,18,21,2,3,2,8,17,18,-3.68,-3.88,-2.00,1     4,40,4,10,12,2,2,2,24,20,23,-3.32,-3.42,-2.00,1     5,43,9,10,12,2,2,1,2,1,29,-3.12,-3.19,-2.73,1     6,52,9,16,19,2,3,2,35,34,35,-3.33,-3.26,-2.95,1     7,46,4,15,18,2,3,2,8,40,42,-3.59,-3.50,-2.00,1     8,40,4,10,12,2,2,2,24,20,46,-2.45,-2.69,-2.00,1


title: Expected standard deviation for a sample from a uniform distribution?

text: I've been trying to find information on the sampling distribution of the standard deviation for uniform distributions and have been having a heck of a time figuring out the expected value for the standard deviation of a sample. Hopefully someone can help out! Basically, I'm trying to find the expected standard deviation given sample size n. Given a uniform, continuous distribution with population range (b-a), I know that the population standard deviation would be (b-a)/sqrt(12); however, what would be the expected standard deviation for a sample size n drawn from that continuous uniform distribution? My hope is that given a Uniform Distribution with range W (or "PopRange" in the example below), I would like to calculate an expected standard deviation for a given sample size. Since the mean sample standard deviation changes as a function of the sample size, I'm hoping to quantify it. Here's some MATLAB code illustrating what I am hoping to achieve.               % Number of samples for the brute-force estimation     totalSamples = 10000;          % The population range (b-a)     PopRange = 30;          for n = 2:100         % Draw samples of sample size n         theSamples = rand(n,totalSamples)*PopRange;              % Calculate the mean of the standard deviations and ranges.  This         % should be  equivalent to the means of the sampling distributions         % of the standard deviation and range         uniformSampleSD(n) = mean(std(theSamples,[],1));         uniformSampleRange(n) = mean(max(theSamples,[],1)-min(theSamples,[],1));              % Calculate the expected standard deviation, note that this is just the         % population standard deviation at the moment and does not incorporate         % the sample size         uniformSampleESD(n) = PopRange/sqrt(12);         % SEE: https://en.wikipedia.org/wiki/Prediction_interval#Non-parametric_methods         % Calculate the expected sample range         uniformSampleERange(n) = (n-1)/(n+1)*PopRange;     end          % Note the red curve underestimating the population standard deviation for     % smaller sample sizes? This is fine, but I would like to find a numerical     % estimate for it     figure     plot(uniformSampleSD,'r')     hold on     plot(uniformSampleESD,'b')      Any help or pointers would be greatly appreciated! * * * Edit: Perhaps I mis-worded what I am looking for? I am looking for the mean of the sampling distribution of the standard deviation for a continuous distribution. If I draw 100000 sets from a uniform distribution [0,1] of sample size n, the distribution of standard deviations will not have a mean of sqrt((1-0)^2 / 12) = 0.289. Empirically, the mean of the sampling distribution of the sd will be:               n = 2, sd_mean = 0.235     n = 3, sd_mean = 0.263     n = 4, sd_mean = 0.273     n = 5, sd_mean = 0.278     n = 6, sd_mean = 0.280     n = 7, sd_mean = 0.282     n = 8, sd_mean = 0.283     n = 9, sd_mean = 0.284     n = 10, sd_mean = 0.285      Note that it asymptotically approaches sqrt((b - a)^2 / 12), but that this is inaccurate for small n. See these examples of the sampling distributions for those small n: http://imgur.com/a/Asgkz (note: the "population" and "observed" sd line labels were accidentally swapped) These were created using this code: totalSamples = 100000;               for n = 2:10         theSamples = rand(n,totalSamples);         setSDs = std(theSamples,[],1);         popSD = sqrt((1-0)^2/12);         nelements = hist(setSDs,50);              figure         hist(setSDs,50)         hold on         plot([popSD popSD],[0,max(nelements)],'r')         plot([mean(setSDs) mean(setSDs)],[0,max(nelements)],'g')         xlabel('Count')         ylabel('Sample Standard Deviation')         title({'Sampling Distribution of the Standard Deviation with 100000 sets';sprintf('drawn from a uniform distribution [0,1] with sample size %u',n)})         legend('Observations',sprintf('Pop SD = %f',popSD),sprintf('Observed Mean Sample SD = %f',mean(setSDs)))     end


title: Acceptable r-square value for multiple linear regression model

text: I'm currently working on my thesis, more specifically I'm analyzing some data collected from researchers about the project's they're working on. In the end, I have performed a multiple linear regression to verify which determinants (budget, team size,...) significantly influence the external impact of a research project. While the dependent variable (project impact) has been assessed from a survey (conducted among the institute's researchers), the independent variables are based on objective measures. **In such a scenario, what would be an acceptable r-square percentage?** Could a low percentage be justified by the "subjectiveness" of the dependent variable? I've heard that the scientific area matters... i.e. in physics you rather need to have a r-square > 95% while in sociology > 5% might be already interesting... Do you have any references about this?


title: statistical method for spatial correlation between images

text: I am working on analyzing a data set and I was wondering what would be the most statistically valid method of demonstrating that there is a strong spatial correlation between images. I have a data set with about 50 pairs of images of cancerous tissue samples. The first image in each pair shows the locations of gold nanoparticles, and the second image shows the locations of the blood vessels in the same tissue sample. By looking at the images it is easy to see that the locations of the nanoparticles match up with the blood vessels, but I would like to prove this statistically in the paper. This is an important point because it demonstrates that the nanoparticles bind specifically to the cancerous areas instead of the normal tissue. I have been looking at different statistics such as a simple linear correlation or something like the answer to this question: Valid method to analyze spatial correlations in images? However, I haven't found anything that would work well for correlation between images. * * * **Edit from Ladislav Nado:** I fabricated two pictures from web...the size and resolution is equal. ![Before treatment](http://i.stack.imgur.com/5ycOr.jpg) ![After treatment](http://i.stack.imgur.com/Q33W5.jpg)


title: Contrasts to ANOVA in R

text: I need to fit orthogonal contrasts to the following model. See the following example data:                  rep <- c(rep( 1, 15), rep(2, 15))         parent1 <- c(1, 2, 3, 4, 5, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4, 1, 2, 3, 4, 5, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4)         parent2 <- c(1, 2, 3, 4, 5, 2, 3, 4, 5, 3, 4, 5, 4, 5, 5,1, 2, 3, 4, 5, 2, 3, 4, 5, 3, 4, 5, 4, 5, 5)         yld <- c(10, 12, 13, 14, 15, 12, 13, 14, 15, 13, 14, 15, 14, 15, 15,11, 13, 12, 13, 15, 14, 12, 13, 15, 12, 15, 15, 13, 15, 16)         mydf <- data.frame(rep, parent1, parent2, yld)         mydf$parent1 <- as.factor(mydf$parent1)         mydf$parent2 <- as.factor(mydf$parent2)      Parent1 and Parent2 are factors.                   mode11 <- lm(yld ~ rep + parent1:parent2, data = mydf)         anova(mode11)       I need to fit contrasts, I suppose the order of the treatments for the coefficients for the interaction is:               Parent1    1, 2, 3, 4, 5, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4         Parent2    1, 2, 3, 4, 5, 2, 3, 4, 5, 3, 4, 5, 4, 5, 5           # contrasts for Parent1:Parent2 interaction      mycontrast <- c(-4, -4, -4, -4, 4,  2, 2, 2, 2, 2,  2, 2, 2, 2,2)      I want to apply above orthogonal contrasts (mycontrasts). I am used to SAS to this short of analysis however not in R.   I have two questions: (a) How can I properly apply contrast to above model? (b) I am concerned about the order of the treatments while writing coefficients. See my assumed orders: are they correct? I assume this should be either alphabetical or ordered in file. Edits: The following my own progress to answer the question               >     mode11$coefficients # to see all coefficients           >     mode11$coefficients[1:5] # just first five to save space here           (Intercept)               rep parent11:parent21 parent1**2**:parent2**1** parent1**3**:parent2**1**           1.500000e+01      1.087998e-15     -4.500000e+00                NA                NA       Should orthogonal contrast ordered in the way the comparisions are ordered here - 1-2, 2-1, 3-1 and so on. As the design is unbalanced you can see some combinations are missing, so should we provide a value like 0 for them.                 contrasts(parent1:parent2) <- cbind(c(-4, -4, -4, -4, 4,  2, 2, 2, 2, 2,      2, 2, 2, 2,2,-4, -4, -4, -4, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2,2 )) Error in `contrasts<-`(`*tmp*`, value = c(-4, -4, -4, -4, 4, 2, 2, 2, : contrasts apply only to factors In addition: Warning messages: 1: In `*tmp*`:parent2 : numerical expression has 30 elements: only the first used 2: In `*tmp*`:parent2 : numerical expression has 30 elements: only the first used


title: How is the inverse gamma distribution related to $n$ and $\sigma$?

text: Given that the posterior estimate of $\sigma'^{2}$of a normal likelihood and an inverse gamma prior on $\sigma^2$ is: $$\sigma'^{2}\sim\textrm{IG}\left(\alpha + \frac{n}{2}, \beta +\frac{\sum_{i=1}^n{(y_i-\mu)^2}}{2}\right)$$ which is equivalent to $$\sigma'^{2}\sim\textrm{IG}\left( \frac{n}{2}, \frac{n\sigma^2}{2}\right)$$ since a weak $\textrm{IG}(\alpha, \beta)$ prior on $\sigma^2$ removes $\alpha$ and $\beta$ from eqn 1: $$\sigma'^{2}\sim\textrm{IG}\left( \frac{n}{2}, \frac{\sum_{i=1}^n{(y_i-\mu)^2}}{2}\right)$$ It is apparent that the posterior estimate of $\sigma^2$ is a function of the sample size and the sum of squares of the likelihood. But what does this mean? There is a derivation on Wikipedia that I don't quite follow. I have the following questions   1. Can I get to this second equation without invoking Bayes' rule? I am curious if there is something inherent in the parameters of an IG that is related to the mean and variance independent of the normal likelihood.   2. Can I use the sample size and standard deviation from a previous study to estimate an informed prior on $\sigma^2$, and then update the prior with new data? This seems straightforward, but I can not find any examples of doing so, or rationale why this would be a legitmate approach - other than what can be seen in the posterior.   3. Is there a popular probability or statistics textbook that I can consult for further explanation?


title: Individual comparisons with logistic regression

text: I have some results which are plotted here: ![](http://dl.dropbox.com/u/22681355/graph.png) I would like to make individual comparisons e.g. does the green line in the top case significantly drop from 2 to 3 and from 3 to 4? I have created contrast coded predictors to test for these differences separately. I would like to know if it is ok to do separate logistic regressions to check for each effect.


title: Generating all possible unique combinations of dataset with two variables

text: I am fairly new to R and I have been reading all related topics but still I don't manage to solve my rather simple problem. I have a dataset with 5 cases for two variables (X,Y), for instance:                 X |  Y     --------     100 | 12     112 | 14     140 | 21         190 | 20      89 |  7      Now I want to replicate this dataset with all possible unique combinations of both variables. If I'm not mistaken there would be !5 = 120 different sets. Now, I have been trying to do this with the 'sample' command. By keeping E constant and then sampling F, I can get a new dataset with different combinations of E and F. The problem of course is that this doesn't return all possible combinations. If I would do this 120 times, I would not get unique datasets. Thanks in advance for your help. Best regards, P