[R-sig-ME] p values through resampling?

[Greg Snow]

You may be interested in this post (and the thread that it is from):
http://finzi.psych.upenn.edu/R-sig-mixed-models/2009q1/001819.html


-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111


> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-
> models-bounces at r-project.org] On Behalf Of Malcolm Fairbrother
> Sent: Thursday, January 27, 2011 11:35 AM
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] p values through resampling?
> 
> Dear all,
> 
> Given the difficulty/controversy of calculating p values for the
> coefficients in mixed models, can someone please tell me whether simple
> resampling is a reasonable alternative means of generating them?
> 
> The example below resamples a level-2 covariate, and the procedure
> seems to generate p values with corresponding t values close to those
> returned by lmer. (A similar procedure works for the level-1
> covariate.) It seemed to me that by comparing p values generated from
> the t statistic (with degrees of freedom ranging from N to n), and a p
> value generated this way, one could generate a conservative test.
> 
> I may be overlooking either (a) good statistical reasons not to do
> this, or (b) previous discussions along these lines. Or this may be
> overkill, when a "abs(t)>2" rule-of-thumb test is enough. I'd be
> interested to hear what people think.
> 
> Thanks,
> Malcolm
> 
> 
> 
> library(lme4)
> n <- 20
> N <- 50
> dta <- data.frame(X1=rbinom(n*N, 1, 0.5), X2=rep(0:1, each=(n*N)/2),
> id=rep(1:N, each=n))
> dta <- within(dta, y <- X1 + X2/10 + rep(rnorm(N), each=n) +
> rnorm(n*N))
> M1 <- lmer(y ~ X1 + X2 + (1|id), dta)
> 
> res <- c()
> resamples <- 1000
> for (i in 1:resamples) {
> 	sim.X2 <- rep(sample(dta$X2[!duplicated(dta$id)], replace=T),
> each=n)
> 	dta2 <- data.frame(dta, sim.X2=sim.X2)
> 	M <- lmer(y ~ X1 + sim.X2 + (1|id), dta2)
> 	res[i] <- fixef(M)[3]
> 	}
> plot(density(res))
> table(fixef(M1)[3]<res)
> qt(1-mean(fixef(M1)[3]<res), 50) # generates t values close to those
> from lmer
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models