title: When does it make sense to reject/accept an hypothesis?

text: In To P or not to P: on the evidential nature of P-values and their place in scientic inference, Michael Lew has shown that, at least for the t-test, the one-sided p-value and sample size can be interpreted as an "address" (my term) for a given likelihood function. I have repeated some of his figures below with slight modification. The left column shows the distribution of p-values expected due to theory for different effect sizes (difference between means/pooled sd) and sample sizes. The horizontal lines mark the "slices" from which we get the likelihood functions shown by the right panels for p=0.50 and p=0.025. ![enter image description here](http://i.stack.imgur.com/c1KYL.png) These results are consistent with monte carlo simulations. For this figure I compared two groups with n=10 via t-test at a number of different effect sizes and binned 10,000 p-values into .01 intervals for each effect size. Specifically there was one group with mean=0, sd=1 and a second with a mean that I varied from -4 to 4, also with sd=1. ![enter image description here](http://i.stack.imgur.com/wtNim.png) (The above figures can be directly compared to figures 7/8 from the paper linked above and are very similar, I found the heatmaps more informative than the "clouds" used in that paper and also wished to independently replicate his result.) If we examine the likelihood functions "indexed" by the p-values, the behaviour of rejecting/accepting hypotheses or ignoring results giving p-values greater than 0.05 based on a cut-off value (either the arbitrary 0.05 used everywhere or determined by cost-benefit) appears to be absurd. Why should I not conclude from the n=100, p=0.5 case that "the current evidence shows that any effect, if present, is small"? Current practice would be to either "accept" there is no effect (hypothesis testing) or say "more data needed" (significance testing). I fail to see why I should do either of those things. Perhaps when a theory **predicted** a precise point value, then rejecting a hypothesis could make sense. But when the hypotheses are of the form either "mean1=mean2 or mean1!=mean2" I see no value. Under the conditions these tests are often being used randomization does not guarantee all confounds are balanced across groups and there should always be the worry of lurking variables, so rejecting the hypothesis that mean1 exactly equals mean2 has no scientific value as far as I can tell. Are there cases beyond the t-test where this argument would not apply? Am I missing something of value that rejecting a hypothesis with low a priori probability provides to researchers? Ignoring results above an arbitrary cutoff seems to have lead to widespread publication bias. What useful role does ignoring results play for scientists? **Michael Lew's R code to calculate the p-value distributions**               LikeFromStudentsTP<-function(n,x,Pobs,test.type, alt='one.sided'){     # test.type can be 'one.sample', 'two.sample' or 'paired'     # n is the sample size (per group for test.type = 'two.sample')     # Pobs is the observed P-value     # h is a small number used in the trivial differentiation     h<-10^-7     PowerDn<-power.t.test('n'=n, 'delta'=x, 'sd'=1,     'sig.level' = Pobs-h, 'type'= test.type, 'alternative'=alt)     PowerUp<-power.t.test('n'=n, 'delta'=deltaOnSigma, 'sd'=1,     'sig.level' = Pobs+h, 'type'= test.type, 'alternative'=alt)     PowerSlope<-(PowerUp$power-PowerDn$power)/(h*2)     L<-PowerSlope     }      **R code for figure 1**               deltaOnSigma <- 0.01*c(-400:400)     type<-'two.sample'     alt='one.sided'     p.vals<-seq(0.001,.999,by=.001)               #dev.new()     par(mfrow=c(4,2))     for(n in c(3,5,10,100)){            m<-matrix(nrow=length(deltaOnSigma), ncol=length(p.vals))       cnt<-1       for(P in p.vals){         m[,cnt]<-LikeFromStudentsTP(n,deltaOnSigma,P,type, alt)         cnt<-cnt+1       }            #remove very small values       m[which(m/max(m,na.rm=T)<10^-5)]<-NA                      m2<-log(m)            par(mar=c(4.1,5.1,2.1,2.1))       image.plot(m2, axes=F,                  breaks=seq(min(m2, na.rm=T),max(m2, na.rm=T),length=1000), col=rainbow(999),                  xlab="Effect Size", ylab="P-value"       )       title(main=paste("n=",n))       axis(side=1, at=seq(0,1,by=.25), labels=seq(-4,4,by=2))       axis(side=2, at=seq(0,1,by=.05), labels=seq(0,1,by=.05))       axis(side=4, at =.5, labels="Log-Likelihood", pos=.95, tick=F)       abline(v=0.5, lwd=1)       abline(h=.5, lwd=3, lty=1)       abline(h=.025, lwd=3, lty=2)       par(mar=c(5.1,4.1,4.1,2.1))                 plot(deltaOnSigma,m[,which(p.vals==.025)], type="l", lwd=3, lty=2,            xlab="Effect Size", ylab="Likelihood", xlim=c(-4,4),            main=paste("Likelihood functions for","n=",n)       )       lines(deltaOnSigma,m[,which(p.vals==.5)], lwd=3, lty=1)       legend("topleft", legend=c("p=.5","p=.025"), lty=c(1,2),lwd=1, bty="n")     }      **R code for figure 2**               p.vals<-seq(0,1,by=.01)     deltaOnSigma <- 0.01*c(-400:400)     n<-10     n2<-10     sd2<-1     num.sims<-10000     sp<-sqrt((9*1^2 +(n2-1)*sd2^2)/(n+n2-2))          p.out=matrix(nrow=num.sims*length(deltaOnSigma) ,ncol=2)     m<-matrix(0,nrow=length(deltaOnSigma),ncol=length(p.vals))     pb<-txtProgressBar(min = 0, max = length(deltaOnSigma) ,style = 3)     cnt<-1     cnt2<-1     for(i in deltaOnSigma ){            for(j in 1:num.sims){              m2<-i         a<-rnorm(n,0,1)         b<-rnorm(n,m2,sd2)         p<-t.test(a,b, alternative="less")$p.value              r<-end(which(deltaOnSigma<=m2/sp))[1]              m[r,end(which(p.vals<p))[1]]<-m[r,end(which(p.vals<p))[1]]+1         p.out[cnt,]<-cbind(m2/sp,p)         cnt<-cnt+1            }       cnt2<-cnt2+1       setTxtProgressBar(pb, cnt2)     }     close(pb)               m[which(m==0)]<-NA          m2<-log(m)               dev.new()     par(mfrow=c(2,1))     par(mar=c(4.1,5.1,2.1,2.1))     image.plot(m2, axes=F,                breaks=seq(min(m2, na.rm=T),max(m2, na.rm=T),length=1000), col=rainbow(999),                xlab="Effect Size", ylab="P-value"     )     title(main=paste("n=",n))     axis(side=1, at=seq(0,1,by=.25), labels=seq(-4,4,by=2))     axis(side=2, at=seq(0,1,by=.05), labels=seq(0,1,by=.05))     axis(side=4, at =.5, labels="Log-Count", pos=.95, tick=F)     abline(h=.5, lwd=3, lty=1)     abline(h=.025, lwd=3, lty=2)     abline(v=.5, lwd=2, lty=1)     par(mar=c(5.1,4.1,4.1,2.1))               hist(p.out[which(p.out[,2]>.024 & p.out[,2]<.026),1],          xlim=c(-4,4), xlab="Effect Size", col=rgb(1,0,0,.5),           main=paste("Effect Sizes for","n=",n)     )     hist(p.out[which(p.out[,2]>(.499) & p.out[,2]<.501),1], add=T,          xlim=c(-4,4),col=rgb(0,0,1,.5)     )          legend("topleft", legend=c("0.499<p<0.501","0.024<p<0.026"),             col=c("Blue","Red"), lwd=3, bty="n")


title: Contingency table analysis to rank preferences of birds per feature

text: I have a dataset that contains observations of objects that female blackbirds carry to their nests. The birds have id tags and the objects are grouped categorically with respect to their `colour`, `size` and `texture`. The code that generates the representative `birds.data` in R is at the bottom. I first look at the contingency tables per feature to see what specific feature categories a bird is more selective on.               # Create Contingency Matrices     col_mat <- as.data.frame.matrix(table(birds.data[, c("id", "colour")]))     size_mat <- as.data.frame.matrix(table(birds.data[, c("id", "size")]))     text_mat <- as.data.frame.matrix(table(birds.data[, c("id", "texture")]))      Below is a graphical summary where I `cbind()` all matrices above. I can observe that `bird 1` is pretty indifferent to the colours of objects she collects but highly prefers `small` and `soft` objects. ![enter image description here](http://i.stack.imgur.com/l9bkl.png) I first ran $\chi^2$ independence tests as below to compare whether the colour (or any feature) preferences of a bird is different from the rest.               # Run chi^2 independence test for bird1 vs rest in colour dimension     bird_1 <-  col_mat[1,]      rest <- colSums(col_mat[2:4,])     data_test <- rbind(bird_1,rest ) # column - bind the two     chisq.test(data_test)      **Question:** However, I do not know how to rank how relatively selective a bird is when comparing different features. For instance, for `bird1`, I would report ranked preferences as : `texture` > `size` > `colour`. How can I achieve that? (I guess the context independent keyword here may be: **comparing different multinomial distributions** ) **Attempt:** Thinking I cannot compare $\chi^2$ test results between features, I decided to look at variances. Below I calculate the MLE's per feature for `bird1` and then look at the (co)variances (`getMultinomCov()`) of the estimated multinomial probabilities.               # Investigate multinom variances per feature for bird1     cov_col <- getMultinomCov(col_mat[1,])     cov_size <- getMultinomCov(size_mat[1,])     cov_text <- getMultinomCov(text_mat[1,])          # Calculating covariance matrix of multinom dist     getMultinomCov <- function (data) {       N <- sum(data)       p_hat <- as.numeric(data / N)       variance         <-  N*p_hat*(1-p_hat)       covariance       <- -outer(p_hat, p_hat)*N       diag(covariance) <-  variance       return <- covariance     }      Would comparing the `cov_col`, `cov_size` and `cov_text` matrices help in finding the "preference ranking" that I am after? If so, can you please demonstrate/explain how? * * *               # Bird Sightings Data Generation     set.seed(2013)     bird1.data <- data.frame( id =  rep(1,200),                               colour = sample(c("Red", "Blue", "Yellow", "Green"), 200, replace=T, c(0.25,0.25, 0.25, 0.25)),                               size = sample(c("Large", "Medium","Small"), 200, replace = T, c(0.1,0.3, 0.6)),                               texture = sample(c("Hard", "Soft"), 200, replace = T, prob = c(0.1,0.8))     )     bird2.data <- data.frame( id =  rep(2,500),                               colour = sample(c("Red", "Blue", "Yellow", "Green"), 500, replace=T, c(0.15,0.45, 0.20, 0.20)),                               size = sample(c("Large", "Medium","Small"), 500, replace = T, c(0.33,0.33, 0.33)),                               texture = sample(c("Hard", "Soft"), 500, replace = T, prob = c(0.55,0.45))     )     bird3.data <- data.frame( id =  rep(3,100),                               colour = sample(c("Red", "Blue", "Yellow", "Green"), 100, replace=T, c(0.9, 0.02, 0.03, 0.05)),                               size = sample(c("Large", "Medium","Small"), 100, replace = T, c(0.6,0.2, 0.2)),                               texture = sample(c("Hard", "Soft"), 100, replace = T, prob = c(0.5,0.5))     )     bird4.data <- data.frame( id =  rep(4,200),                               colour = sample(c("Red", "Blue", "Yellow", "Green"), 200, replace=T, c(0.4,0.3, 0.2, 0.1)),                               size = sample(c("Large", "Medium","Small"), 200, replace = T, c(0.1,0.1, 0.8)),                               texture = sample(c("Hard", "Soft"),200, replace = T, prob = c(0.3,0.7))     )     birds.data_c <- rbind(bird1.data ,bird2.data,bird3.data,bird4.data)# one frame to rule them all     birds.data <- birds.data_c[sample(nrow(birds.data_c)),] # and in the darkness shuffle them.


title: Forecasting using arima(1 1 1)

text: I want to forecast using `arima(1 1 1)`. I have done everything but when I used the R code:               modelforecast <- forecast(model,2)      the output says:               Error in .cbind.ts(list(e1, e2), c(deparse(substitute(e1))[1L],      deparse(substitute(e2))[1L]),  :    not all series have the same frequency.      Can someone help me how to fix this problem?


title: A gallery of charts, diagrams, and plot types

text: What would you recommend as a comprehensive gallery of data presentation techniques? A source that could be used to refer to while you're thinking about better ways of presenting your data? I've identified the following ones, but will be glad if you could add yours: Online galleries:   * http://www.mathworks.com/discovery/gallery.html   * http://www.idlcoyote.com/gallery/   * https://developers.google.com/chart/interactive/docs/gallery?csw=1   * http://www.walkingrandomly.com/?p=4788   * http://en.wikipedia.org/wiki/Category:Statistical_charts_and_diagrams (does not provide one-page graphic gallery)   * http://docs.ggplot2.org/current/   * http://www.itl.nist.gov/div898/handbook/graphgal.htm   * http://scikit-learn.org/stable/auto_examples/index.html Books (plots scattered across pages):   * Edward R. Tufte, _The Visual Display of Quantitative Information_   * Nathan Yau, _Data Points_


title: What are the statistical implications of assigning treatment and control in time series study under 2 different approaches?

text: I have a project wherein my group wants to measure the effect of counseling at-risk students on their academic outcomes. There are 1,000 students who will be enrolled during the current academic term. Resources (such as the counselors) are limited, so we can only use a `subset` of these. My boss says that we have enough funding to take the `top 6 deciles` (600 students with the worst grades) and randomly assign 50% to `control` and 50% to `treatment`. The 50% who receive `treatment` will be contacted by advisers, whereas the 50% in `control` will not. However, this is not just a `pre/post` study. Instead, we will have `k` instances of outreach (i.e. `treatment`) during the term, where `k` is between 2 and 5. My boss laid out 2 different approaches that we could take and asked which I recommend. I'm want to be sure that I understand the statistical implications of each approach, so that I can make a well-informed recommendation. Here are the two choices:   1. Randomly assign `treatment` and `control` to all students at the beginning of the study, select the top 6 `deciles`, at `k = 1`, and have the counselors talk to these same students at each of the `k` time periods. Within the 6 `deciles`, about 300 should be in `treatment` and about 300 in `control` and those groups will be compared at the end of the study.   2. Randomly assign `treatment` and `control` to all students at the beginning of the study, select the top 6 `deciles`, at `k = 1`, and have counselors reach out to the approximate 50% in the `treatment` group, just as in option #1. **BUT** at `k = 2` and each subsequent time period, recalculate the top 6 `deciles`. Many of the students will be the same, but some of those who were in the top 6 `deciles` at `k = 1` may be in the bottom 4 `deciles` at `k = 2, k = 3, etc`. Assignment to `treatment` or `control` won't change, but inclusion in the subset of the population will change for some students. My intuition is that option #1 will me more simplistic and thus shield us from potential bias that we may introduce. On the other hand, #2 would give us a `larger sample` with more `identifying variation`, though not all `treatment` students will receive treatment at each of the `k` time periods. In some ways, it seems like #2 could actually allow for a more in-depth analysis, because we could not only look at the effect of `treatment`, but also the `frequency` and `timing` of treatment, so that we could potentially pick up on `diminishing marginal returns` to treatment or variations in the effectiveness of `treatment` in terms of timing (i.e. outreach by counselors could be more effective earlier on or at the end of the course). Would #2 introduce bias or create other statistical complications? If so, should I rule out #2 or would it require additional statistical controls to account for the (semi-)`endogenous` changes in the composition of the 6 deciles?


title: Is there any way to model this grouping process?

text: I've been working on this problem for three days now, and it doesn't seem that it can be solved using pure thinking. Maybe there is some distribution that can model this but I couldn't find any solution till now. This is the process I'm trying to model: I have $K$ points randomly dispersed in an area of radius $R$ using 2D-poisson random variable. I start by putting the first point in a group. Then, the distance $d$ between the first and the second points is found. If this distance is larger than a certain threshold ($Dth$) then the second point will be put in the same group as the first point. Otherwise, the second point will be put in a group alone. The process continues until the last point $K$. A point will be put in a group if it has its distance from all the points in that group larger than $Dth$. I've derived $P(d>Dth)$ which is the probability that two points have a distance between them larger than $Dth$, so let's denote it by $p$, and $P(d<Dth)$ by $q$. What I want to find is either the average number of groups resulting given $K$ points (don't care about the radius $R$ because it is incorporated in the given probabilities) or the average number of points per group given $K$. Any idea will be greatly appreciated. Thank you very much


title: Nonparametric test for trend using Python

text: I am looking to perform a nonparametric test for trend on a continuous outcome across three groups, preferably in Python. For example height (pretend height is not normal) in 4th, 5th and 6th graders. I would like to implement something like the Cuzick method. Scipy has Wilcoxon rank sum and other nonparametric methods but only for two groups. Similarly, Scipy has a Kruskal-Wallis method for three groups but it does not indicate direction or trends. Does anything like this exist for exploring a directional trend across three groups? * * * To clarify I am trying to determing whether there is a significant shift in a continuous trait measured across three groups. The groups will be of very different size: group 1 has 1000's of samples and likely to be normally distributed, group 2 100's os samples, group 3 ~10 or less. Group 1 serves as the "control" group, and my hypothesis is that the mean value of group 1 will be shifted in either direction relative to group 0, and group 2 will be shifted further **in the same direction as group 1**. Because group 3 will always be very small compared to the other group, my instinct was to use nonparametric methods, but I am open to other suggestions. Can anyone suggest a method to explore this type of directional trend? Thanks


title: How to show that the mean is (weakly) consistent

text: How can I show that the mean is weakly consistent? Is weakly consistent the same as consistent?


title: Showing count data can be dichotomized without losing too much information

text: I am working with a set of data that has accident counts by year from taxi drivers. R code:   no.accident <\- c(16003,2355,4433,18823)   one.accident <\- c(172, 26, 232, 9)   two.accident <\- c(0,0,3,0)   df =data.frame(no.accident=no.accident,one.accident=one.accident,two.accident=two.accident) year 1 there were 16003 accident-free drivers, 172 divers had 1 accident, no driver had 2 accidents, year 3 4433 drivers were accident-free, 232 drivers had 1 accident and 3 drivers had 2 accidents. I wish to show that the data can be analyzed either as a count data or a "0 vs at least 1 accident" binomial data. Could anyone kindly suggest a way to statistically show it--any hypothesis test or regressions?


title: ANOVA failed to model data, what is a more appropriate test?

text: This is a cross-post (http://stackoverflow.com/questions/19432964/anova-error- in-levelsxx) about an error I received in R while trying to run an ANOVA on my data. But error aside, I need help understanding why an ANOVA can't deal with my data and what other statistical models could be applied instead. So here's my objective: I have 3 people (speaker) who recorded a bunch of words that I analyzed. The analysis yielded 3 continuous variables: skewness, kurtosis and Center of Gravity (CoG)*. I need to find out what combinations of these 3 variables best model the difference between each speaker. For example, are skewness and CoG together more significant than just CoG in finding the difference between speakers? I have a basic knowledge of stats, but erring on the side of assuming I'm an idiot might be better for any complex explanations. Thanks in advance!   * The skewness is a measure for how much the shape of the spectrum below the center of gravity is different from the shape above the mean frequency.   * The kurtosis is a measure for how much the shape of the spectrum around the center of gravity is different from a Gaussian shape.   * The center of gravity is a measure for how high the frequencies in a spectrum are on average weighted by their energy.


title: Symmetry of alternative hypotheses of pairwise t-tests with p-value correction?

text: So I'm using R's `pairwise.t.test` on 4 populations with a Hochberg p-value correction and an alternative hypothesis that the population on the row has a higher mean than the population on the column. I get the following output:               > pairwise.t.test(data, labels, alt='greater', p.adjust.method='hochberg')             A       B      C         B 0.00323  -      -          C 0.00019  0.4864 -          D 0.00135  0.9955 0.9970      Clearly we can reject the null hypothesis and accept the alternative hypothesis for pairs `(B, A)`, `(C, A)`, and `(D, A)`. We can also not reject the null hypothesis for pairs `(D, B)` and `(D, C)` Now about `(B, D)`: If we were using no p-value correction, I would be confident in rejecting the alternative hypothesis that B has the same mean as D and accepting the null hypothesis that B has a high mean than D. However, the p-value correction causes me to be suspicious of this. My question is of the "give me a fish and then teach me to fish" variety: Does the same symmetry exist with the p-value correction as without? Also, can anyone point me toward a resource that would deepen my understanding of p-value correction to the point where the answer to my first question is obvious.


title: How to use a fitted model parameters for forecasting other time series

text: I have fitted a ARIMA(1,1,2) to time series TS1 as below:               arima112<-arima(TS1, c(1,1,2))      now I want to use the coefficients of ar and ma that I got from `arima112` to forecast another time series (TS2). How can I apply the `arima112` model on TS2? Thanks


title: Controlling for a variable in a GEE using SPSS

text: I have a dataset which I think would suit GEE. I am trying to assess the relative contribution of a number of simultaneous predictors and their interactions to the intensity of hurt feelings. I asked people to report on one hurtful event in each of 3 relationships, as well as fill out a number of scales. I then coded those hurtful events into a number of categories. I have a measure of the person who hurt the participant as a 3 level ordinal RM variable, the type of hurtful event as a 7 level nominal variable, and 2 measures of personality traits (hurt-sensitivity and self-esteem) as continuous variables. I also have a variable which I would like to control for, which is the time since the hurtful event occurred (a continuous variable). (Please let me know if I have explained this dataset clearly.) I have run a GEE using the 4 main predictor variables (relationship with person who hurt the participant, hurtful event type, hurt-sensitivity and self-esteem) in SPSS, and it looks like it all makes sense and it worked (where ahurtat=intensity of hurt feelings, arellvl=relationship with person who hurt participant, atoe7_rc=type of hurtful event, aHPS_ctrd=hurt- sensitivity and aRSE_ctrd=self-esteem):               GENLIN ahurtat BY arellvl atoe7_rc (ORDER=ASCENDING) WITH aHPS_ctrd aRSE_ctrd       /MODEL arellvl atoe7_rc aHPS_ctrd aRSE_ctrd arellvl*atoe7_rc arellvl*aHPS_ctrd arellvl*aRSE_ctrd          atoe7_rc*aHPS_ctrd atoe7_rc*aRSE_ctrd aHPS_ctrd*aRSE_ctrd arellvl*atoe7_rc*aHPS_ctrd          arellvl*atoe7_rc*aRSE_ctrd arellvl*aHPS_ctrd*aRSE_ctrd INTERCEPT=YES      DISTRIBUTION=NORMAL LINK=IDENTITY       /CRITERIA SCALE=MLE PCONVERGE=1E-006(ABSOLUTE) SINGULAR=1E-012 ANALYSISTYPE=3(WALD) CILEVEL=95          LIKELIHOOD=FULL       /EMMEANS TABLES=arellvl SCALE=ORIGINAL       /EMMEANS TABLES=atoe7_rc SCALE=ORIGINAL       /EMMEANS TABLES=arellvl*atoe7_rc SCALE=ORIGINAL       /REPEATED SUBJECT=ID_number SORT=YES CORRTYPE=UNSTRUCTURED ADJUSTCORR=YES COVB=ROBUST          MAXITERATIONS=100 PCONVERGE=1e-006(ABSOLUTE) UPDATECORR=1       /MISSING CLASSMISSING=EXCLUDE       /PRINT CPS DESCRIPTIVES MODELINFO FIT SUMMARY SOLUTION COVB.      However, I am not sure how to control for the time since the hurtful event occurred (not included in the above model). I know you can create a Type I model with the GENLIN function, but I can't work out a way to specify that I want the time since the event occurred to be in the first level of the model and the other 4 predictors to be in the second level. Any help would be greatly appreciated.


title: Treatment of outliers in annual time series data

text: I have an annual time series of data of a growth-rate variable $X$ for 50 years. Most of the values for the variable $X$ are less than 10%. The exception are two values that are around 30%. How do I treat these two values in the regression? I would like to seek your suggestions in this regard.


title: ARIMA and linear regression

text: I have a linear regression model that is used to forecast the 'afluent natural energy' (ANE) of some region. The predictors for this model are:   * the previous month ANE (`ANE0`)   * the previous month rain volume (`PREC0`)   * the current month forecast for rain volume (`PREC1`) We have 7 years of historical data for all of these variables, for each month. The current model just runs a OLS linear regression. I feel there's a lot of improvements to be done, but i'm not a time series specialist. The first thing I notice is that the predictors are highly correlated (multicollinearity). I'm not certain of the impacts of multicollinearity on prediction confidence. I decided to try a time series approach, so I ran a ACF and PACF on the historic data: The ACF shows a sine wave pattern, and the PACF has a spike at 1 and 2. So I tried both `ARIMA (2, 0, 0)` and `ARIMA(2,0,1)` to predict 20 periods ahead. The ARIMA(2,0,1) shows good results, but I'm not certain as to how to compare it to the linear regression model. What's the best way to test the performance of these model? I'm using R as analysis tool (together with the `forecast` package).


title: How to cope with multicollinearity and interactions between IVs in generalized linear models?

text: I have made a generalised linear model with a single response variable (continuous/normally distributed) and 4 explanatory variables (3 of which are factors and the fourth is an integer). I have used a Gaussian error distribution with an identity link function. Do I need to check for multicollinearity and interactions amongst explanatory variables? If yes, how do I do this with categorical explanatory variables?


title: Expected effect of each single-nucleotide polymorphism (SNP) in a genome-wide association study (GWAS)

text: It was mentioned in a genetics class that in a genetic association analyses of a trait with all SNPs, it is possible to compute the expected effect of each SNP with the trait using the correlation structure between the SNPs and the size of the effect at the true causal variant. Can you explain why this is so?


title: Expected values of (X*Y^2) when X and Y are dependent normally distributed RVs

text: I have two _jointly_ normal variables X and Y with mean both zeros and variances $\sigma^2_{X}$ and $\sigma^2_{Y}$ separately, the covariance is $\sigma_{XY}$. Now I want to calculate the expected value of $Z=X*Y^{2}$, $E(Z)$. Any ideas? Thanks.


title: Adjustment for missing values of the categorical variables in a data set

text: I Have a data set containing about 40 categorical variables. I am trying to factor analyze them. But each categorical variable contains a good number of missing values. Some of them are simply because of non-response. The respondent did not filled up any answer option for that question. Some of them are due to questions of the following type: > 5) Do you create formal work teams in your institution? >   >   >      1="NO" 2="YES" >   > > (Please skip question number 6 and 7 whose answer to this question is > 1="NO") > > 6) How many members form the work team? (for example) > > 7) What is the criterion of selecting team members? (for example) Now those who answered "NO" for question number 5 will not answer 6 and 7. He will again start from 8. This is another source of missing information or gap in the data set. Because of specially this type of missing values if I omit missings listwise a lots of information is missed. So, I am looking for adjusting these missing values. I don't know how to adjust these missing values (of both non-response and the second type I mentioned) for categorical variables. Taking mean, median or even EM algorithm may not be appropriate for categorical variable I guess. So, what should be done and how? My actual number of observations is 212, but it reduces to only 42 when I use `na.omit(data).`


title: Why do we rely on the standard error?

text: > **Possible Duplicate:**   >  Difference between standard error and standard deviation I cannot come to terms with the fact that the standard error is used to describe the **accuracy** of a point estimate. For example, when making inferences from a sampling distribution about the **accuracy** with which the sample mean estimates the population mean, we rely on the standard error. Wouldn't the standard deviation of the sampling distribution be just as appropriate if not more so at indicating the **accuracy** of the sample mean? In other words, why use the n^(1/2) in the denominator of the standard error calculation?


title: What is difference between standard error for a large sample and population standard deviation?

text: Does it make a difference between the two terms for a very large sample-size from a target population?


title: Why does the standard deviation not decrease when I do more measurements?

text: I made 100 measurements of a certain quantity, calculated mean and standard deviation (with MySQL), and got mean=0.58, SD=0.34. The std seemed too high relative to the mean, so I made 1000 measurements. This time I got mean=0.572, SD=0.33. I got frustrated by the high standard deviation, so I made 10,000 measurements. I got mean=0.5711, SD=0.34. I thought maybe this was a bug in MySQL, so I tried to use the Excel functions, but got the same results. Why does the standard deviation remain high even though I do so many measurements?


title: variance, standard deviation and standard error

text: When describing the variance of a sample data set, we can use `variance` and `standard deviation`. `standard deviation` is the square root of `variance`. Sometimes, I also heard `standard error` What are the difference between `standard error` and `standard deviation.` Is there any function to compute `standard error` within R?


title: Standard Error vs. Standard Deviation of Sample Mean

text: I'm having difficulty figuring out the difference between (typically referred to as S.E): A. The Standard Error And B. The Standard Deviation of the Sample Mean (typically referred to as s) Are they the same thing? OR   1. S.E is the standard deviation of the the mean of the sampling distribution similar formula as the standard deviation except you use n-1 instead of n in the denominator   2. s is the standard deviation of a single random sample -- same formula as the standard deviation


title: Factor scores vs. construct mean scores in regression analysis

text: I have 48 items in my questionnaire that represent 8 constructs. After conducting an exploratory factor analysis (EFA) with a principal components extraction method and Varimax rotation method, 8 sets of factor scores (FAC_1 to 8) were computed and saved using the regression method. I've read in academic papers that those factor scores can be used as variables in regression analysis, but the problem is when I do so, the regression test yields no result! (I receive only 0 values!). Alternatively I computed the mean score of the items representing each construct (resulting in 8 new variables) and used them in the regression test. This method worked, but i'm not sure if it's a correct procedure or not. ### Questions   * **Why am I getting zero-values for regression coefficients when I use factor saved scores and meaningful regression coefficients when I use computed mean scores?**   * **Can the mean scores of each construct instead of its factor scores generated through EFA be used in multiple regression analysis?**


title: Probability of classification given two observations

text: If two doctors each correctly diagnose a disease 55% of the time, and both agree in a certain case there is disease, what is the probability of disease? Surely it is higher than 55%?


title: Random Forest checklist

text: I am building a random forest model in R. Based on my research I (hope to) have come up with quite some understanding about how they work, and more importantly when they work. I simply would like to see my understanding of RF-models cross-validated here. So here we go (true or false?):   * scaling is not necessary in Random Forest models.   * when dealing with imbalanced data, one could downsample/upsample/use weights. In package `randomForest` however, the option classwt seems to be unreliable? Therefore I simply downsample my data (50/50) using package Unbalanced. In my final model I will downsample and build the forest K times and take the average of the predictions. Is this wise?   * Random Forests have little problems with highly correlated variables. I have many ideas for new features and would like to include all of them at once and based on importance (MeanDecreaseGini) optionally decide to leave some out. But the correlation between predictors will not influence performance? in other words: the model with extra features will perform at least similar to the one with less variables?


title: Mathematical sound way to rediscretize a distribution

text: Case: I have a discrete distribution with a certain interval level. I want to rediscretize it so that the interval level is reduced, but the general characteristics of the distribution remain intact. Example: I have a discrete distribution with 10 elements, ie: x= 0 1 2 3 4 5 6 7 8 9 y= 0 1 4 9 16 25 36 49 64 81 Now I want to reduce this, so that I am left with, say, 6 elements with the same range of x: x= 0 .. .. .. .. 9 y= ... My Question: Is there established (and easily understood) mathematical theory/methods to accomplish this? And where can I find this (incl. simple examples)? My initial thoughts: Apply a kind of percentile approach relating to the number of new elements I wish to have in the new discrete distribution. Ie, if I want 5-elements, then maybe take the 0th, 25th, 50th, 75th and 100th percentile of the original distibution. (BUT... I havent found any proven theory supporting this approach)


title: Assigning meaningful cluster name automatically

text: The objective of my work is to cluster the text documents. Once the documents are clustered, traditionally the system will assign numeric value for the clustered group. For example if I have 5 categories, then the clustered documents are labeled by any one of these numerical values {1,2,3,4,5}. I would like to assign the cluster name (e.g., philosophy, biology,...) automatically rather than labeling it as {1,2,3,...} for further processing. My initial idea is to provide the cluster name by counting the highest frequency word in that cluster. I am confused if this approach is good or not. I am using k-means clustering. Currently I am excluding LDA (Latent Diriclet Allocation) or other methods.


title: How do you validate your machine learning models?

text: I am wondering what approaches are commonly used for validating machine learning models designed for classification or prediction tasks: Approaches that am using at the moment: **Using truth-sets:** \- ROCs, Bootstrapping, Accuracy, Sensitivity, Specificity, Cross-validation **Orthogonal validation:** \- Use a different class of algorithm that can perform prediction or classification task and compare results Any other suggestions ?


title: Test to show when diverging linear regression models are statistically different

text: What method will allow me to identify the values of the independent variable where two (or more) diverging linear regression models are statistically different? e.g. Take two simple linear regression models derived from some experimental test: $y_1 = -1x+100$ and $y_2 = -0.5x+100$. These could represent two time points on a dose vs survivorship relationship, or two different materials in a temperature vs strength relationship. ![enter image description here](http://i.stack.imgur.com/d5Z7Q.jpg) At $x = 0$ the two models give the same $y$. As you move away from $x = 0$ the values of $y_1$ and $y_2$ diverge, however for some period their confidence intervals will overlap suggesting that $y_1$ is indistinguishable from $y_2$ in this region. Taking the temperature vs strength relationship as an example, this suggests that there is no difference in strength of the two materials over this temperature range. Therefore, I would like to be able to calculate this range, or conversely identify the range of $x$ values when $y_1$ is statistically different to $y_2$. My understanding of predictions bands is not great, but I don't think this is the solution. The only solution I can come up with is to calculate the $y$'s and CI's at discrete values of $x$ and run a series of tests. I am looking for something more elegant!


title: Can one apply different transformations to the same variable within different analyses?

text: My stats questions concerns the appropriateness of using different transformations of the same variable within different analyses (within the same larger study). My research involves emotional expression as a dependent variable. I'm looking at 1) baseline differences between healthy controls (HCs) and individuals with a disease and 2) treatment effects within the disease group (Treatment Groups: HC-no-treatment, Disease-no-treatment, Disease-treatment 1, Disease-treatment 2). Each subject was assessed at baseline, post-treatment, and follow-up. **Research Aim 1. Baseline Differences Between Health Status Groups (Healthy Controls and Disease).**   When looking at the emotion expression data, at baseline, by health status group (2: healthy controls and Disease), the raw data is not normally distributed (based on Shapiro-Wilk). Because square root transformation normalizes the distribution, I was planning on using a t-test (as opposed to Mann-Whitney) to assess baseline differences between the healthy control and disease groups. **Research Aim 2. Treatment Effects.** When looking at the emotion expression data, by treatment group (HC, disease-no-treatment, disease-treatment1, disease-treatment2) at each time (baseline, post, follow-up), the raw is not normally distributed (based on Shapiro-Wilk). For these groupings, the square root transformation does NOT normalize the distribution. However, a log10 transformation does normalize the distribution. Main Question: Can I use a square root transformation for Aim 1 and a different transformation for Aim 2 (e.g.,log10 in a mixed model ANOVA, or raw data if I end up using nonparametric stats for Aim 2)?


title: Joint Confidence Interval

text: Given the following model: ${\rm Wage}_i=\beta_0+\beta_1 {\rm Married}_i+\beta_2{\rm Female}_i+\beta_3 {\rm Married}_i \times {\rm Female}_i + \varepsilon_i$, I am interested in finding $\beta_1+\beta_3$ for different quantiles of the dependent variable (i.e. on average, how much more/less does married females make compared to single females for a given quantile). I already have the results of the quantile regression model, and I calculated $\beta_1+\beta_3$ for all quantile levels. However, I am trying to find a confidence interval for $\beta_1+\beta_3$ to determine whether this quantity is statically significant for a given quantile level. One thing that came to mind is the Bonferroni joint confidence interval. Using this method I was able to find a lower bound for the joint confidence interval for each of $\beta_1$ and $\beta_3$, but I am not sure how to find a confidence interval for the sum of the two covariates. Would it be correct to sum the two confidence intervals that we got (I didn't find evidence that support this approach)? Is there a way to determine such a confidence interval? Appreciate your help. Thanks.


title: Alternative nonparametric test for chi-square test for independence

text: suppose we have a table like this: \begin{array}{ |l | c | r |} \hline & \mathrm{win} & \mathrm{lose} \\\ \hline \mathrm{male} & X_{11} & X_{12} \\\ \hline \mathrm{female} & X_{21} & X_{22} \\\ \hline \end{array} $%![enter image description here][1]$ and we want to know if winning or losing depends on sex (male or female) apart from chi-square test for Independence, which non-parametric tests we can use? Thanks


title: Plotting the decision boundary of a kernel SVM (RBF)

text: Suppose we are given a training set of 2D points that are linearly non- separable. I train a binary SVM with an RBF kernel in order to classify them. What I want to do is to draw the desicion boundary. The kernel SVM I train leads to a decision function of the form: $$ f(\mathbf{x})=\sum_{i=1}^{N_s}\alpha_i y_i k(\mathbf{x}, \mathbf{x}_i)+b, $$ where $N_s$ is the number of support vectors, $\mathbf{x}_i$, $\alpha_i$, and $y_i$ are the $i$-th support vector, the corresponding positive Lagrangian multiplier, and the associated truth label, respectively. I use the LIBSVM package to train the SVM, thus all the above are known from the created model file. If we extend the above decision function, we see that it's a sum of Gaussians centred at the support vectors, i.e., $$ f(\mathbf{x})=\sum_{i=1}^{N_s}\alpha_i y_i \operatorname{exp}\left(-\gamma\lVert\mathbf{x}-\mathbf{x}_i\rVert^2\right)+b, $$ where $\gamma=\frac{1}{2\sigma^2}>0$ is a user-defined parameter. How could I plot in the same space of my input data (the Euclidean 2D space), the contour of the decision function? Could you please provide some equations or ideas? I've tried some things, but I had no luck! Thanks a lot!


title: What are the statistical features of a network/graph without community structure?

text: I want to filter out networks/graphs which are just a single community rather than interconnected communities. Each community can either be a random graph or a star graph, but there should no community structure. The definition is loose considering that these data sets can have multiple communities with a varying number of edges connecting them together or nodes acting as bridging islands. (these networks are connected) I thought of looking for networks where the average path length for at least 10% of the nodes is greater than 2. I am using `MATLAB' and any reference to useful libraries would be appreciated.


title: How to deal with multicollinearity when performing variable selection?

text: I have a dataset with 9 continuous independent variables. I'm trying to select amongst these variables to fit a model to a single percentage (dependent) variable, `Score`. Unfortunately, I know there will be serious collinearity between several of the variables. I've tried using the `stepAIC()` function in R for variable selection, but that method, oddly, seems sensitive to the order in which the variables are listed in the equation... Here's my R code (because it's percentage data, I use a logit transformation for Score):               library(MASS)     library(car)          data.tst = read.table("data.txt",header=T)     data.lm = lm(logit(Score) ~ Var1 + Var2 + Var3 + Var4 + Var5 + Var6 + Var7 +                  Var8 + Var9, data = data.tst)          step = stepAIC(data.lm, direction="both")     summary(step)      For some reason, I found that the variables listed at the beginning of the equation end up being selected by the `stepAIC()` function, and the outcome can be manipulated by listing, e.g., `Var9` first (following the tilde). What is a more effective (and less controversial) way of fitting a model here? I'm not actually dead-set on using linear regression: the only thing I want is to be able to understand which of the 9 variables is truly driving the variation in the `Score` variable. Preferably, this would be some method that takes the strong potential for collinearity in these 9 variables into account.


title: Comparing principal components if number of variables changes

text: I have a database that contains daily stock returns of more than a thousand stocks for many decades. I would like to achieve the following goals:   1. Construct a **time-varying** measure of co-movement (or common component) of stock returns.   2. Compare the evolution of this measure over time, and later relate it to other variables.  The main idea is to see whether _the importance of the common component of stock returns systematically varies over the business cycle_. My first idea was to run a PCA in each quarter, save the the share of explained variance by the first principal component, and then this quarterly time series of explained variance would be the measure of co-movement over time. The problem, however, is that the number of variables for which I run the PCA is not constant over time. In the early part of the sample I have, say, 200 stocks and in the more recent part, say, 1000. I assume it doesn't make sense to compare the explained variance by PC1 if the number of variables is different over time. I would appreciate any ideas to tackle this problem. Is there any measure that can be compared if the number of variables changes over time?


title: How to "undo" feature scaling/normalization for output?

text: I'm normalizing (or standardizing or feature scaling) my neural network training inputs and training targets. I just doing linear scaling and the formula I'm using is: `I = Imin + (Imax-Imin)*(D-Dmin)/(Dmax-Dmin)` where `I` is the scaled input value, `Imin` and `Imax` are the desired min and max range of the scaled values, `D` is the original data value, and `Dmin` and `Dmax` are the min and max range of the original data values. In my case I'm setting `Imax` to 1 and `Imin` to -1. I'm trying to predict a continuous real-valued output. How do I scale the output of my network back into the "unscaled" range?


title: How to get rid of heterogeneity of regression slopes using multilevel modeling?

text: I have problem with homogeneity of regression slopes in ANCOVA. I have two species of snakes and I want to compare their tail lengths (dependent variable), that are dependent on body length (in ANCOVA used as covariate). Other ANCOVA assumptions are met, but there exist significant interaction between species*bodylength (intercepts are comparable, slopes are different). Andy Field´s SPSS guide tells that broken assumption of homogeneity of regression slopes could be cast aside using multilevel statistic model. I´m trying to correctly run this analysis whole day, but I am unable to do it. Problem is that I want to compare tail lengths between species, which is my top of hierarchy (at the bottom are individual specimen). Do I have to set variable species as subject (assuming some hierarchy)? Undoubtely I have to set body length as covariate and tail length as dependent variable. How can I compare my groups (species) independently of slopes (slopes are random)? What should be set as fixed factors and what random factors in SPSS? Field´s guide also tells that for mixed model analysis there is bug (version 17.0) related to factors dialog box such that categorical variables (in this case my species variable) should be put into covariate box (at least if it is bivariate). Is he right? I know what I want to do but dont know how to do it.


title: Comparison of two classifiers based on precision/recall/F1 only?

text: For two classifiers **h1** and **h2** I have the precision, recall and F1 score as a percentage (along with the original labeled data set that they were tested on). If I had access to which samples each classifier classified right/wrong, I would be able to do, for example, McNemar's test to evaluate significance, but unfortunately I don't. I would ideally like to be able to speak on the significance of the results obtained by **h2** , that is, whether **h2** is a significant improvement over **h1**. Am I unable to do that, or is there _something_ I can say using only precision/recall/F1 and the labeled data set?


title: Association between 2 categorical data (using SPSS)

text: In a sample of 250 cardiac patients who have had a myocardial infarction and just left hospital, I asked them how important it was to them to know several issues, such as "How important is it for you to know when to return to work, or when to revisit the hospital" etc. Answers were provided ranging from "not important at all" to "very important" (Likert-type scale with 5 answers, assuming all answers have "equal" distances between them). I would like to know:   1. How to treat the answers to the above questions (as categorical?).    2. How to look at whether there is any association between the above type of questions and other questions that are categorical (for example: "Who do you prefer to give you information": a-physician, b-nurse, c-other; or "Level of education": a-primary level, b-secondary level, c-bachelor, or equivalent d-postgraduate degree).  This questionnaire hasn't got any cummulative score or subscale scores that could be used for testing associations.


title: Critical appraisal of survival paper

text: I must make a critical appraisal of a cancer survival paper. I hope someone here can give some hints. The paper investigates person-level markers of socio-economic status as exposure and looks for association of these with survival and also comorbidities, disease/prognostic factors and treatment using Cox proportional hazards for survival and logistic regression for mortality. I think the paper is a good one and I can't find the criticisms. I am not sure if Cox proportional hazards is the best model for use; maybe a full paramteric model of relative survival is better, but I don't known. Thank you for a hint if you can. http://www.ncbi.nlm.nih.gov/pubmed/22315055


title: Polynomial kernel function

text: Consider SMV with the polynomial kernel $k(x_1,x_2)=(\langle x_1, x_2\rangle + 1)^d,$ where $d > 1.$ Is it true that if the dataset is separated with a hyperplane then the SVM (with the kernel $k$) always separates the dataset with a hyperplane (without any errors in the dataset)?


title: Standard deviation of multiple root mean square values

text: I am calculating the R.M.S. of a periodic signal that has a stochastic component to it. For every period I am able to calculate a value for R.M.S. using the following function: $R.M.S. = \sqrt{\frac{1}{N} \sum_{i=1}^N x_i^2 }$, which I understand I can use only if the time series has a mean of zero, in which case the above is equivalent to the standard deviation of that time series. However, I have multiple periods, each of which has a value for R.M.S. Presumably, I can calculate the mean and standard deviation over many periods. Mathematically, what would be the correct way to express this?


title: Probability of moving house (non-independence problem)

text: We are analyzing a dataset where we are looking at the probability of individuals moving house, and the factors that predict moving. The problem is that a lot of people are married, and (usually) if one moves house, one's spouse also moves house (although not everyone is married). We have two options, neither of which seems satisfactory:   1. Model at the household level: but then including predictors at the individual level is hard (not everyone is married).    2. Model at the individual level, and use something which corrects for clustering. Doesn't seem to account for the nearly complete clustering, and I don't feel like the stimates will be appropriate. In addition, the probability of moving is curious. If a couple moves, and a single person doesn't, was that a 0.5 probability of moving (because 1/2 households moved - as in (1) above, or 0.66, because 2/3 households moved, as in (2) above.


title: Quantifying Change in a Histogram Valued Timeseries

text: I'm attempting to do binary classification where my raw features are collections of histograms that are recorded in a time series. These histograms are scaled to sum to 1. To be more precise and define some notation, let $H_{t}$ be a histogram at some time $t$ and $H_{t}(i)$ it's value in the $i$th bucket (here the $i$th bucket is produced by binning non-negative integer values). I have such an $H_{t}$ for each one of my data points (so I could write $H^k_t$ to represent the histogram associated with data point $k$ at time $t$). It is somewhat hopeless to use the histograms as features as what seems to be most relevant to my classification is the change that is occurring in the histograms as time changes. My current strategy to quantify this change is to, for each $i$, fix bucket $i$ and then preform a discrete wavelet transform along $t$ using the haar wavelet ($t$ takes a power of 2 many values so the DWT is easy to apply). While this approach works somewhat well it is constrained by it's inability to pick up simultaneous and general changes in the histogram; for instance it is heuristically probable in my dataset that a decrease in the value of bucket $i=2,3$ is 'bad' unless it is accompanied by an increase in the values of some larger buckets, that is to say that the histogram is shifting strictly right. My question is then: are there are any general techniques/strategies that describe the change of histograms over time that capture these more general changes? I should note that the ranges of $i$ and $t$ are fairly small, less than 10 each.


title: Summated scales (after factor analysis)

text: 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/factor analysis is not supported by literature, for example if I have one or two factors that combine some variables in a different form than from what I have found in the literature, thus I also name them differently)


title: Sobol sequence for a set of discrete numbers

text: Is it possible to create a sobol sequence for a set of numbers? I have 5 discrete variables x1, x2, x3, x4 and x5 and say for example x1 = {1,2,3,4,5}, x2 = {10,30,40,50,60}, x3 = {1,3,5,7,9}, x4 = {30,45,60,90}, x5 = {2,4,6,8}. If it is possible to create a sobol sequence, how will be the VBA routine look like? I have tried to sample using Rnd function in Excel and it works fine.


title: Time interval histogram (tih), Fano factor, coefficient of variance

text: I'm given this: A pacemaker neuron fires 10 spikes/s in a fairly regular manner, with interspike intervals (ISIs) that are distributed according to $\mathcal N(100, 5)$. I'm asked to draw the time interval histogram (TIH), estimate the Fano factor (FF), and then calculate the coefficient of variation (CoV). **My questions:**   1. The TIH, this just looks like a regular "bell curve" centered at the mean ISI, 100, with $SD = \sqrt{5}$, right?    2. For the FF, which I think is about the spike count (i.e. the number of spikes/s?), the mean is 10 but how can I estimate the variance? How do I get from the ISI variance to the variance in spikes/s?    3. It seems to me that the CoV can be calculated since I'm given the variance and mean of ISIs (I.e., ${\rm CoV}= SD/{\rm mean} = \sqrt{5}/100$).


title: Loadings in factor analysis

text: 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)?


title: Poisson autocorrelation function

text: I want to draw an autocorrelation function, call it, AC(tau) where tau is the offset in the autocorrelatin. The vector I'm feeding into the AC is the spike train of a neuron whose action potential fire in a poissonian manner, i.e., independent of each other and equally probably at each time step. I imagine that at tau (time delay)=0, AC=1, and that AC thenceforth decays except for smaller and smaller peaks at multiples of 1/lamba (the average interarrival interval). Is that right?   Now, what if it's a non-homogenous poisson process? Then what will the AC look like? Is it possible to generalize?


title: Summated scales for factors (latent variables)

text: 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)?


title: How to compute a latent variable (as a composite index)

text: 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)? Or what formula can you recommend me for obtaining a latent variable by taking into consideration all the factor obtained from factor analysis?


title: center variable separately for each factor in linear mixed model?

text: I am working on a relatively simple mixed model where I have two continuous predictors with an interaction term and three sites. I am treating the two predictors as fixed effects and the site as random. If I were to center the predictors by subtracting the mean of each, would I do so separately for each site? Or better to calculate the global mean of each predictor and use that to center. Alternatively, I am thinking of centering based on a physical value, a threshold, that is different between each site. Is this valid? I am using nlme in R; the call is `lme(response ~ predictor1 + predictor2 + predictor1:predictor2, random = ~1 | site)`


title: How random are the results of the kmeans algorithm?

text: I have a question regarding the **kmeans algorithm**. I know kmeans is a randomized algorithm, but how random is it and what results can I expect. Suppose you have clustered a dataset into $4$ clusters, where each point has identity $1$, $2$, $3$, or $4$ (which tells you what cluster it belongs to). You then perform a second clustering on the same dataset with the same criteria.   1. Will all the points in a specific cluster say during the first clustering be in the same cluster the next time you apply the kmeans algorithm?   2. If not, will they most likely be in the same cluster and is there some measure for this likelihood? By some output I received in R, I believe that 1. does not hold as I get different cluster sizes for different runs on the same dataset. All help is greatly appreciated!


title: R function breakpoints finds more breaks than Fstats

text: I am using the `strucchange` package to detect structural changes in my data. The plot of my Fstats shows several peaks, yet if I use `breakpoints` one additional breakpoint is found. My F statistic at that point suggested by `breakpoints` is even below the threshold. Why does it find more?


title: Which statistical analysis should I use for the study

text: I am doing a research on 20 participants where 10 will be from computers background and rest will be from non computing background. Each participants will attend a lab session and the performance of the participants will be evaluated based on a rubric. Each row of the rubric is evaluated on scale of 1 to 4. The research question question is: does backgroundof the participant effect to play in learning about xxxx. Which statistical method should I use to analyze the data?


title: A problem of probabilities, calculation of 2 events

text: I have this problem: A student must choose exactly two out of three electives: art, French, and mathematics. He chooses art with probability 5/8, French with probability 5/8, and art and French together with probability 1/4. What is the probability that he chooses mathematics? What is the probability that he chooses either art or French? Now my rationale was this: p(A and F)+p(F and M)+p(A and M)=1 something has to happen. 1/4+p(F and M)+p(A and M)=1 p(F and M)+p(A and M)=3/4 p(M)=p(F and M)+p(A and M) p(M)=3/4 But how do I get p(A) or p(F). Thanks.


title: How to get real-valued continous output from Neural Network?

text: In most of the examples I've seen so far of neural networks, the network is used for classification and the nodes are transformed with a sigmoid function . However, I would like to use a neural network to output a continuous real value (realistically the output would usually be in the range of -5 to +5). My questions are:               1. Should I still scale the input features using feature scaling? What range?     2. What transformation function should I use in place of the sigmoid?      I'm looking to initially implement it PyBrain which describes these layer types. So I'm thinking that I should have 3 layers to start (an input, hidden, and output layer) that are all linear layers? Is that a reasonable way? Or alternatively could I "stretch" the sigmoid function over the range -5 to 5?


title: How to analyse pre and post test data when the outcome measure is Yes/No data?

text: I am doing research on 20 participants and the participants will work on a lab module. There will be a pretest before the lab module and post test after the lab module on concepts related to flowchart, pseudo code , etc., The possible outcomes for the question is Yes/No. **Which statistical method should i use to analyze the pre - post test yes/no data?**


title: Why is the sum of squared differences from the sample mean smaller than the sum of squared differences from the true mean?

text: So the description I read is: > It is important to note that the formulas for calculating the variance and > the standard deviation differ depending on whether you are working with a > distribution of scores taken from a sample or from a population. The reason > these two formulas are different is quite complex and requires more space > than allowed in a short book like this. I provide an overly brief > explanation here and then encourage you to find a more thorough explanation > in a traditional statistics textbook. Briefly, when we do not know the > population mean, we must use the sample mean as an estimate. But the sample > mean will probably differ from the population mean. Whenever we use a number > other than the actual mean to calculate the variance, we will end up with a > larger variance, and therefore a larger standard deviation, than if we had > used the actual mean. This will be true regardless of whether the number we > use in our formula is smaller or larger than our actual mean. Because the > sample mean usually differs from the population mean, the variance and > standard deviation that we calculate using the sample mean will probably be > smaller than it would have been had we used the population mean. Therefore, > when we use the sample mean to generate an estimate of the population > variance or standard deviation, we will actually underestimate the size of > the true variance in the population because if we had used the population > mean in place of the sample mean, we would have created a larger sum of > squared deviations, and a larger variance and standard deviation. To adjust > for this underestimation, we use n - 1 in the denominator of our sample > formulas. Smaller denominators produce larger overall variance and standard > deviation statistics, which will be more accurate estimates of the > population parameters. I understood none of this as it seems contradictory. It says that when using a number other than the actual mean (to me, e.g. the sample mean) will be different then the population mean causing a larger variance. In the middle it says the variance and sample mean will be smaller if the population mean had been used... Can you please explain which it is exactly and why so?


title: Expectation operator and logarithmic function

text: $$\frac{1}{C_t}=E_t\left[\beta \frac{1}{C_{t+1}} \right]R_{t+1}$$ How to log linearise the function? $C_{t+1}$ is the stochastic term; $\beta$ is known.


title: Residual autocorrelation in Poisson (Neg.Bin) models - Durbin Watson test

text: I am running some Poisson (or Neg.Binomial depending on overdispersion) models and i want to check for residual autocorrelation due to the nature of the data (monthly cases). I am using R and i am producing correlograms with the use of the "acf" command, where i can check if the autocorrelation for the different lags falls outside the produced confidence intervals. On the other hand based on published papers that have done similar studies to mine i see that they have used the Durbin-Watson test to check for presense of first or second order autocorrelation which is of interest here. In R i have found two commands for this test: the "dwtest" in package lmtest (based on what is called the pan algorithm) and "durbinWatsonTest" in package car which uses bootstrap. In both these tests the documentation suggests they are used for lm objects which is not the case for my models (nor the published ones). My question is: Is the D-W test valid to use for residual autocorrelation in Poisson (Neg.Bin) models? Is there some other test for these models? If the results between the "acf" and the two D-W tests matched i would feel more comfortable but they don't and to my understanding the "acf" is more reliable.


title: Poisson regression on the means of count data

text: I just finished a small research project about hummingbirds and the effect of hummingbird feeders. I am a bit unsure about how to proceed with the statistics. We placed 15 points in a distance gradient away from the feeders, where we sampled visitation rates, pollination and bird / floral abundance (all count data). We sampled all points twice, what I would like to do is to use Poisson regression on the mean of the two samples. Does that makes statistical sense? First, I would guess that using the means would push the data distribution towards normality. Second, can you use Poisson regression on non whole numbers?


title: Where did the words bias and variance come from?

text: I understand that a bias model is more relaxed while a model with a lot of variance is more flexible. but, where did these terms come from, and, why 'bias' and why 'variance'?


title: Using extreme value theory to estimate bounds

text: Suppose I have I have a random variable $X$ that I know is doubly bounded on support $[0,\theta]$ but I dont know $\theta$ (we don't know anything on the distribution of $X$, but assume it is not crazy so usual regularity conditions apply when you need them). Consider $Y:=X+Z$, where $Z$ is standard normal ($X$ and $Z$ are independent). And suppose I draw $Y_1,Y_2,...$ which are i.i.d. readings of $Y$. Using the answer posted (actually the appendum to that answer) for this (extreme value theory: show normal to gumbel), and assuming I did my computations right, I showed that $Y_{(n)}$ when suitably normalized with the $(a_n,b_n)$ given in the **appendum** of the solution, converges to the Gumbel as well. My question is: Is there a way to use this process to provide some estimates on $\theta$? I feel that before normalizing $Y_{(n)}$, it is possible that we have a shifted, scaled version of Gumbel with the shift reflecting the value of $\theta$ somehow (intuitively, as $n$ gets larger, the index $j$ for which $Y_j = Y_{(n)}$ is more and more likely to also be such that $X_j = X_{(n)}$ which is getting closer and closer to $\theta$, so for gigantic $n$, perhaps $Y_{(n)}$ is somewhat like $Z_{(n)} + \theta$?) I can't seem to mathematically justify this or find a procedure to extract some estimates of $\theta$. Any suggestions?


title: Ensembles of Ensembles?

text: I like the idea of ensemble learners, especially Bagging, but I always wondered as why they are not the most powerful learners since they have a clean motivation. I don't have the answer to that question but I had another idea. Normally in Bagging people use the same classifier for learning. So they divide the dataset into slices and for each slice they train a classifier of the same type (e.g. logistic regression) and then they use voting. But my question is why not to use ensembles of ensembles? Why not to create a bagging classifier of logistic regression, a bagging classifier of SVM, a bagging classifier of ANN, a bagging classifier of random-forest, and then use voting. So each classifier is an ensemble and then all the ensembles become an ensemble. Then use voting again. Has anyone tried this before? Papers? ... etc? There must be!


title: How to prevent false positive findings with repeated correlations?

text: I have to run repeated correlations. If one of them reach the significant value (as given in the Pearson table), how can I be sure that it isn't actually a false positive? I'm pretty new in the statistical analysis, so I apologize if my question seems naive. Thanks for any answers.


title: randomSurvivalForest in R

text: I'm using the randomSurvivalForest package for R, version 3.6.4. I have been using it for a project for a while, with no problem. However, now I have added some additional predictors to my dataset, and find that now some of them are not scored (importance). It looks like 100 is the magic number; before I had 95 and they were all scored; now I have 110 and the output reports importance values for 100. I have experimented with some of the parameters, including -ntree- and -mtry-; different predictors are omitted, but always 10 of them. Am I missing something obvious? Or maybe not so obvious? There is no warning or other indication that rsf() has done something with my predictors, and they are still there in the dataframe when rsf() has finished. Thanks for any ideas.


title: Sparse ELM vs SVM

text: What's the difference between SVM and Sparse Extreme Learning Machine with Gaussian kernel proposed in the following paper:http://www.ntu.edu.sg/home/egbhuang/pdf/Sparse-ELM-IEEE-T- Cybernetics.pdf As use of kernel in ELM will result in zero randomness and all nodes are treated as support vector, so what will be the difference between this type of sparse ELM and conventional SVM


title: Learning parameters of a mixture of Gaussian using MLE

text: It seems that MLE (via EM) is widely used in machine learning / statistics to learn the parameters of a mixture of Gaussians. I'm assuming we're given random samples from the mixture. My question is: Are there any proven _quantitative_ bounds on the error in terms of the number of samples (and perhaps the parameters of the Gaussian)? For example, what is the runtime required to estimate the parameters up to a certain error? Ideally, these bounds would not assume that we start in a local neighborhood of the optimal solution or any such thing. (If EM is not the method of choice and there is a better way of doing it, please point this out, as well.)


title: Automatically produce summary by factor variable in R

text: I've got a dataframe like the following:                  case simulation temp plank oxygen     1     1          1    8     7     11     2     2          1   16    10     15     ...     17   17          2   26    12     17     18   18          2   15     8     12     19   19          2   28    11     21     20   20          2   24     6     14      I'd like to get summaries split by the levels of the simulation variable. For example, I'd like the mean of `temp` for simulation==1 and simulation==2, and the same for standard deviation. At the moment I'm using the following code, which is pretty horrible:               df <- read.csv("data.csv")     attach(df)          # Create subset variables     temp1 = subset(temp, simulation==1)     temp2 = subset(temp, simulation==2)          plank1 = subset(plank, simulation==1)     plank2 = subset(plank, simulation==2)          oxygen1 = subset(oxygen, simulation==1)     oxygen2 = subset(oxygen, simulation==2)          print(sd(temp1))     print(sd(temp2))     print(sd(plank1))     print(sd(plank2))      I'm sure there must be an automatic way to do this in R, but I can't find it. I've tried using `summary(df ~ simulation)`, but that doesn't produce anything useful.


title: Best feature selection method for naive Bayes classification

text: i want to make classification with naive Bayes. I have got about 100 Features. Numerical ones as well as categorical ones. Since i want only the most relevant ones to be included for the classification task i want to find them with some kind of feature elimination. My question now is the following: what is the method to use that for (paper/reference?!) and is this method implemented in some sort of software package. Since i use R i would especially prefer some R package. Thanks in advance for your help!


title: What are the odds of drawing 7 cards that end up sequential - from a 52 card deck? (7 card poker straight)

text: While playing a friendly game of Texas Hold'em poker, a player drew a 7 card straight. Although in texas hold'em a player may only use 5 of the possible 7 cards, the discussion about odds immediately came up. What are the odds of getting a 7 card straight? For those unfamiliar with poker, the question can be asked this way: What are the odds of drawing 7 cards from a 52 card deck, with those cards ending up sequential? * I tried searching the vast knowledge of the internet, but was unable to come up with an answer. the closest i get is the probability of drawing a 5 card straight in a 5 card stud game (0.00392465), but i got lost trying to add the probability of the next 2 cards - due to the complexity of the straight (the next 2 cards can complete the straight - if the first 5 cards drew **4.5.7.8.9** and the next 2 cards were **6.10** ). Any help or pointers on this subject would be extremely helpful. calculating a straight *) An Ace card can be used to start a low straight or to complete a high straight - both **A.2.3.4.5.6.7** and **8.9.10.J.Q.K.A** are legal. but it cannot be used to wrap - _J.Q.K.A.2.3.4_ is not legal.


title: Coefficient of Determination: For the perimeter and area of a square: Why different?

text: When calculating the coefficient of determination for a square, why is it that if you use the data set for the side length of as X= (1,2,3,4) and the perimeter as Y=(4,8,12,16) the Coefficient of Determination equals 100%. But when doing the same calculation for using the area of the square with the X points (1,2,3,4) and the Y data point (1,4,9,16) the Coeff of Determ decreases to 96%? Any insight would be great. The simpler the language the better too. Thanks!!


title: Poisson probability of observing at least one zero out of k independent cases

text: I have a set of observations from a biological experiment. I have counted the occurrences of a particular event on each of the chromosome for $n$ cells and finally I have the final mean value for each chromosome.               2.2, 6.4, 3.4, 10.2, 4.45      So I have 5 chromosomes making on an average of above mentioned number of events. I have calculated the Poisson probability of any of the above chromosome becoming zero. I have used R function: `ppois(q=0,lambda=wt[1:5])` and got probabilities `1.108032e-01 1.661557e-03 3.337327e-02 3.717032e-05 1.167857e-02` , which shows there is total probability of 0.15 so to observe any chromosome with zero event. Now I want to know at what values, these probabilities become zero or negligible. How can I proceed to achieve this?


title: Include interaction regressors separately?

text: I was taught that if I include a regressor in an interaction, then I should include it separately, too. That is, if I regress $y$ on $x_1 \times x_2$, then I should also include $x_1$ and $x_2$ separately. Although I have never seen this spelled out in a text, this makes sense to me -- if I include the "cross partial," then I should include the first order effect, too. However, in published economics papers (in good journals) I see a lot of $$y_i = \alpha + \beta x_{1,i} + \gamma x_{1,i} \times x_{2,i} + \epsilon_i$$ when it seems that $$y_i = \alpha + \beta x_{1,i} + \gamma x_{1,i} \times x_{2,i} + \delta x_{2,i} + \epsilon_i$$ is more correct. Is it OK to omit the non-interacted $x_2$? Is there a good reference on this? Thanks!


title: Deployment process for Classification models (i.e. decision trees)

text: We inherited a turnover prediction process from a previous analyst. The tool we are using is SPSS Modeler. The model being used is CHAID decision tree. The target variable is binary: "stayer" vs "departer" (e.g. person stays with company, person leaves the company). The process is detailed below under BACKGROUND. **I am looking for this specific feedback from the CV community:**   1. What do you see inheritly wrong with the below process?   2. What are some better ways of arriving at a predication using a decision tree classification model? ( I am only looking at decision tree models because I am new to data mining and predictive modeling and these models are easier to understand) I should mention that accuracy has so not been a problem thus far. But what prompted this question is the fact that everyone is classified as "departer" (100% false positives). I posted a separate question about that here. I can post a link to a SPSS Modeler .str file (and related text files, all anonimized) if anyone wants to see it. **BACKGROUND** Task: Predict number of terminations for three of our directorates (e.g ~ “business units”) for upcoming fiscal year using two years historical terminations and headcount data. Datafile: Employee data with several demographic fields. Every record is an employee that either stayed during the entire two year period or terminated (i.e. "stayer" vs "departer") Our current process involves runing a CHAID Model interactively using SPSS Modeler for a recent known time period. For example 2009-2010 (file setup as decribed above). > Sidenote: This is the "training" data. The data is not partioned into > "training" and "testing" data in the traditional classification sense. This > is mainly because the count of "departers" is relatively small compared > (1:14) to "stayers". We do this as follows: 1) CHAID interactive model is run iteratively (trial and error), each time measuring the performance by this subprocess: > a. Apply the model to beginning of year population (e.g. 2011 in this case) > to come up with propensity scores for each individual. > > b. Sum the propensity scores by organization and divide by two. We divide by > two since terminations being used to train the model cover a two year period > but we are only interested in projecting a year out. This produces > predictions by directorate. > > c. Once the accuracy is above 85% move on to next step. 2) Build new model using model from #1 as a guide . New model will be similar to one in #1--business knowledge and some guesswork is applied to adjust splits...but this part is minimal. 3) Apply final model to beginning of year population (e.g. 2012..assuming that's the upcoming year) to come up with final predications by directorate. BTW I tried looking for similar posts on model deployment process but only found posts such as the following which didn't hit the mark for me help-with- deploying-a-model


title: Question about proof for luce choice axiom w.r.t. conditional probability

text: In Luce (1959) the choice axiom is definied, that for a finite subset $T$ of $U$ such that, for every $S\subset T$, $P_S$ is defined.   * If $P(x,y)\ne 0,1$ for all $x,y\in T$, then for $R\subset S\subset T$: $P_T(R) = P_S(R)P_T(S)$   * If $P(x,y)= 0$ for all $x,y\in T$, then for $S\subset T$: $P_T(S) = P_{T-\\{x\\}}(S - \\{x\\})$ Now it is written that (Lemma 2, page 7) > If $P(x,y)\ne 0,1$ for all $x,y\in T$, then this axiom is equivalent to > $P_S(R) = P_T(R|S)$, for $R\subset S\subset T$. It is written that this result is obvious except for the condition $P_T(S)>0$. **My (probably naive) question is:** why is it obvious? The conditional measure indcued by $P_T$ is only defined if $P_T(S)>0$ by: $P_T(R|S) = \frac{P_T(R\cap S)}{P_T(S)}$. Again, why is it enough to proof, that $P_T(S)>0$?


title: Implementing an ETS and ARIMA forecast

text: I've been using the R Forecast package which I have used to create fitted ETS and Arima models. I can easily predict ahead within R using the Forecast package but need to be able to do prediction inside a browser using Javascript. A tiny worked example would really help me in working through the maths. How do I implement the ETS and ARIMA model as a set of equations?               ets(y = object)             Smoothing parameters:         alpha = 0.376          gamma = 1e-04             Initial states:         l = 155.9654          s=26.8689 16.1568 6.5561 -3.4351 -22.9801 -33.7052                -54.1355 -3.8163 16.3216 6.0061 20.9192 25.2434            sigma:  0.0277          ARIMA(0,0,0)(0,1,0)[12]                              sigma^2 estimated as 62.5:  log likelihood=-83.68


title: Avoid local minima in kmeans

text: Many machine learning techniques suffer from the curse of local minima, one of them is K-means. I am using a matlab script for a computer vision task. One of the first steps I do is kmeans clustering on an image in the Lab color representation, I only consider the ab components. So I consider my pixels to be a set of ab data points, without using any gemetrical information (i.e. pixel location) of the pixels. However, I guess these are less important details; in essence I am just doing k-means. My _statistically non-deterministic algorithm_ (read Matlab script) can give a (c) correct result or a (w) wrong result. Fortunately, (c) seems to happen more often, but (w) also happens :(. I can only verify this by hand, the computer is more or less clueless that it computed a wrong answer. However, the final error of k-means seems to be slightly, but consistently higher in case of (w) compared to case (c). This is the result of running the algorithm multiple times on the same image with the same parameters (k=7). That is: Runs resulting in (c):               Replicate 1, 36 iterations, total sum of distances = 1.53591e+06.     Replicate 1, 130 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 56 iterations, total sum of distances = 1.53591e+06.     Replicate 1, 129 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 45 iterations, total sum of distances = 1.53591e+06.     Replicate 1, 122 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 42 iterations, total sum of distances = 1.53591e+06.     Replicate 1, 80 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 106 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 125 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 120 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 128 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 108 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 129 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 112 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 120 iterations, total sum of distances = 1.54903e+06.     Replicate 1, 121 iterations, total sum of distances = 1.54903e+06.      Runs resulting in (w):               Replicate 1, 80 iterations, total sum of distances = 1.61651e+06.     Replicate 1, 33 iterations, total sum of distances = 1.61651e+06.     Replicate 1, 62 iterations, total sum of distances = 1.62508e+06.      So my instinct tells me that the prior probability of a (c) is a lot higher than that of a (w), which is good. Since now I could just repeat k-means n times, resulting in a nxk-means algorithm, and select the result with lowest final error. This will likely result in a (c). However, this seems a rather crude solution. Are there smarter/faster/statistically _sounder_ ways to tackle local optima in k-means?


title: Ratio of sum of Normal to sum of cubes of Normal

text: Please help me to find the limiting distribution (as $n \rightarrow \infty$) of the following: $$ U_n = \frac{X_1 + X_2 + \ldots + X_n}{X_1^3 + X_2^3 + \ldots X_n^3},$$ where $X_i$ are iid $N(0,1)$.


title: How to calculate the variance of the GED distribution?

text: The density of the GED distribution is given by \begin{align} GED(l;\mu,\beta,\nu)&=\frac{\nu \exp\left[-(\frac{1}{2}) \left|\frac{l-\mu}{\beta \lambda}\right|^\nu \right]} {\lambda 2^{(1+1/\nu)}\beta\Gamma (1/\nu)} \end{align} where $\Gamma(\cdot)$ is the gamma function and \begin{align*} \lambda&=\left[ 2^{(-2/\nu)}\Gamma(1/\nu) / \Gamma(3/\nu)\right]^{\frac{1}{2}} \end{align*} how to calculate the variance?


title: Dependent variable maximum value contingent on independent variable

text: I am trying to create a model for debt collections. In the past I have used logistic regression to predict pay/no-pay. This has worked well but has a few unfortunate consequences. People are more likely to pay when the outstanding balance is low. Thus, the model tends to push low balances to the good tail. The real objective, then, is to maximize dollars collected at a particular depth of file. The default strategy is to order accounts by the debt amount which is considered a "dummy" solution. Improving on this dummy solution has been more difficult than I anticipated. The problem is that the dollars collected is contingent on the debt amount! When I try to model dollars collected, the debt amount completely overpowers everything else in the model. Is there another way to "subtract" the influence of debt amount or otherwise control for it? The desired result is to have the model rank order dollars collected better than simply using the debt amount. These are things I have tried:   * 2-stage model using Heckman correction predicting pay/no-pay and then dollars collected   * Built separate models predicting pay/no-pay on high, medium and low debt amounts using logistic regression. Then centered each model using their segment odds.   * Linear regression on just the payers predicting dollars paid and applied the model to the entire population.   * Straight logistic regression predicting pay/no-pay


title: "How-to": omit p-value and use N-P approach in practice?

text: I've been reading a variety of articles 'bout p-value controversy and since the best way to understand theory is to practise it, I devised an example. Suppose, I'm a psychologist and I conducted an experiment. My study was simple - the aim was to check the level of empathy in two groups: men and women. From some theory and earlier studies, I know, that the distribution of empathy is normal in both groups and standard deviation is equal to 1. I want to test that H_0: μ_men = μ_women, against H_1: μ_men not equal μ_women. If I want to use Neyman-Pearson approach and omit all that p-value issue, what should I literally do? Recently, so that I can remember, what was inside, I've read those articles:   * Lew, M.J., "To P or not to P: on the evidential nature of P-values   and their place in scientic inference"   * Lehman, 1993, "The Fisher, Neyman-Pearson Theories of Testing Hypotheses"   * Hagen, R.L., "In Praise of the Null Hypothesis Statistical Test"   * Cohen, J., "The Earth is round, p>0.05"   * Hubbard, R., "The Earth is Highly Significantly Round (p<0.0001) [this is a comment on the article of Cohen]   * Schmidt, Hunter, "Are There Benefits From NHST"   * Dienes, Z., chapters 2,3 of "Understanding Psychology as a Science" Well I'll try to calculate as @Glen_b suggests in a few days. I hope I can present my results here.


title: Partial F ratio from ANOVA table

text: In multiple regression, if you have just an ANOVA table, and nothing else, no specific data, how can you do a partial F test on X1, given X2 is already in the model? So, you have the ANOVA table:               source        df   SS    MS     F     -----------------------------------     regression    2    1.44  0.72   9.72     error         3    0.22  0.07     total         5    1.66     -----------------------------------      All values are filled in. With only this information, how can you do the partial F test where:   * F = MSR(X1|X2) / MSE(X1, X2)   * MSR(X1|X2) = SSR(X1, X2) - SS(X2) = 1.44 - ????   * MSE (X1, X2) = MSE = 0.07 SSR(X1, X2) can be obtained from the table (SS regression) MSE(X1, X2) can also be obtained from the table (just MSE) but I cannot get SS(X2) from the table, as far as I know...... As far as I know, you need specific X and Y values to do this. Any other way from just the table?


title: Test for variation using random effects model

text: I am think that it is possible to analyse **a model with just random effects** but I am not sure as I have never done it. I am looking for guidance on whether it is appropriate, what assumptions I need to be aware of, and how to do it properly. From my study of an insect;   * I have a response variable (age at death, "age")    * Two treatments ("Treat1" and "Treat2") both of which have two levels (Treat1 has "A" and "B", and Treat2 has "P" and "Q")    * There is also 40 genotypes (1-40)    * With four replicates (w,x,y,z) of each combination of Genotype/Treat1/Treat2    * Each replicate contains 50 individuals Put simply, my data looks like 32000 rows of this:               Treat1  Treat2  Genotype  Block  Individual   Age        A       P       1         w      1            23     A       P       1         w      2            35     A       P       1         w      3            44     .       .       .         .      .            .     .       .       .         .      .            .     .       .       .         .      .            .     B       Q       40        z      50           76           I would like to know if each combination of Treat1 and Treat2 (AP,AQ,BP,BQ) have genetic genetic variation - i.e. is there variation between my 40 genotypes within each treatment combination? I think I need a model for each of AP, AQ, BP, and BQ, along the lines of               Age ~ Genotype [ Treat1 == "A" & Treat2 == "P"] * Block [ Treat1 == "A" & Treat2 == "P"]      Where Genotype and Block are random effects. I hear Gamma distribtions are better to use in lifespan (time to death) models. **My questions are:** a. Is this an appropriate way to show whether or not my genotypes have variation? b. Can I build the four models as defined above or is that a really poor way of doing it? c. If possible, what functions should I be using in R (lm, glm, lmer... & summary, summary.lm, aov, anova...)? d. What should I expect, if gamma is more suitable than gaussian, to see when I compare `plot(model)` for gamma compared to gaussian? * * * This is currently my model...               AP= df$Treat1=="A" & df$Treat2=="P"     apmodel<- lmer(df$Age[AP]~(1|df$Genotype[AP])+(1|df$Block[AP]))     summary(apmodel)      Which I think is right but I'm not sure what to do with the output..               > summary(apmodel)     Linear mixed model fit by REML      Formula: df$Age[AP] ~ (1 | df$Genotype[AP]) + (1 | df$Block[AP])             AIC   BIC logLik deviance REMLdev          57343 57371 -28667    57336   57335         Random effects:          Groups           Name        Variance Std.Dev.          df$Genotype[AP]  (Intercept) 17.23798 4.15186           df$Block[AP]     (Intercept)  0.15416 0.39263           Residual                     93.18777 9.65338          Number of obs: 7757, groups: df$line[AP], 40; df$Block[AP], 4          Fixed effects:                 Estimate Std. Error t value     (Intercept)  49.9948     0.6939   72.05      **Is there genetic variance??**


title: Comparing clusters from related data

text: I am not sure what statistical test to use for comparing clustering on related data. I have two sets of measurements on each person. I cluster each person by each measurement. I then measure how much the two coincide. I want to know how I would test for how significance the agreement of the clustering is. For example I have 10 people and for each person I have two sets of attributes say vector A and Vector B both 10 of length 10. I cluster the 10 people into 2 clusters using vector A and get something like 1 2 1 1 2 1 1 1 2 2 That is person one is assigned to cluster one, person two to cluster two person three to cluster one etc I cluster the 10 examples into 2 clusters using vector B and get something like 1 2 1 1 2 1 1 1 1 1 These two clustering agree to 80% but how significant is this? What is the formula for calculating this? Is there a variant for more than 2 clusters? would it make sense if there were a different number of clusters for each measurement or if some samples only had one set of measurements? Thank you for any thoughts :)


title: Use random forest outliers to detect group of variables

text: I have a input data and an output binary variable . The y value is 1 if the patient get ill.               > summary(x_or)        Symscore1        Symscore2        exercise3     exerciseduration3                       groupchange        age3            Min.   :0.0000   Min.   :0.0000   Min.   :1.000   Min.   :0.000     Regular to Regular          : 340   Min.   :45.00        1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:3.000   1st Qu.:2.000     Regular to Menopausal       : 360   1st Qu.:49.00        Median :0.0000   Median :0.0000   Median :4.000   Median :3.000     Transitional to Transitional: 171   Median :54.00        Mean   :0.5504   Mean   :0.5941   Mean   :3.651   Mean   :2.545     Transitional to Menopausal  :1492   Mean   :54.07        3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:5.000   3rd Qu.:3.000     Menopausal to Menopausal    : 246   3rd Qu.:59.00        Max.   :5.0000   Max.   :5.0000   Max.   :5.000   Max.   :4.000                                         Max.   :66.00          packyears           bmi3            education3        Min.   : 0.000   Min.   :16.77   Basic     : 348        1st Qu.: 0.000   1st Qu.:22.32   Highschool:1013        Median : 0.000   Median :24.84   University:1248        Mean   : 4.397   Mean   :25.60                          3rd Qu.: 5.714   3rd Qu.:27.72                          Max.   :97.143   Max.   :57.09                         > summary(y_or)        0    1      2129  480       I have fitted a random forest model and computed the outlier measure                   rf = randomForest(x = x_or, y = y_or, proximity = T, ntree = 1000)     out = outlier(x = rf$proximity, cls = y_or)          plot(out, col=y_or)      From the outlier plot I can see that some of the samples with y=1 have a very low outlier measure. The general prediction performance of my model are very low. Can I deduce from this plot that there is a subgroup of patient that will very likely get ill? Basically I would like to say at the medical team that it is difficult to obtain a general model but that we are able to classify people with this characteristics... ![enter image description here](http://i.stack.imgur.com/oIg7P.png)


title: table with preselected columns?

text: Is there any way to make a table with preselected columns? For example, let's say v=c(1,2,3,2,3,2), table(v) would produce something like this, 1 2 3 1 3 2 but I want a table like this 1 2 3 4 1 3 2 0 which means it has to have columns with zero frequency...