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

[Malcolm Fairbrother]

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