m=50 k=5 n=m*k group=rep(c(1:m),rep(k,m)) #chooseCRANmirror(graphics=F) library(lme4) #install.packages("lme4") nsim=1000 tau2=rep(0,nsim) sigma2=tau2 truesigma2=1 for (i in 1:nsim){ U=rnorm(m,sd=sqrt(truesigma2)) eps=rnorm(n) y=eps+rep(U,rep(k,m)) fit=lmer(y~(1|group)) summary(fit) variances=VarCorr(fit,type="varcov") variances tau2[i]=attr(variances,"sc") sigma2[i]=variances$group[1,1] } mean(tau2) mean(sigma2) par(mfrow=c(1,2)) hist(tau2) hist(sigma2)#only 5 repetitions dev.copy2pdf(file="m5n5sigma21.pdf") dev.copy2pdf(file="m50n5sigma21.pdf") dev.copy2pdf(file="m50n5sigma20.pdf") #parametric bootstrap truesigma2=0.2 U=rnorm(m,sd=sqrt(truesigma2)) eps=rnorm(n) y=eps+rep(U,rep(k,m)) fit=lmer(y~(1|group)) fitH0=lm(y~1) summary(fit) summary(fitH0) sigma2star=tau2star=rep(0,nsim) Q=rep(0,1000) for (i in 1:nsim){ y=simulate(fitH0)[[1]] fitH0star=lm(y~1) fitstar=lmer(y~(1|group)) Q[i]=-2*(logLik(fitH0star)-logLik(fitstar)) variances=VarCorr(fitstar,type="varcov") tau2star[i]=attr(variances,"sc") sigma2star[i]=variances$group[1,1] } hist(sigma2star) #pvalue observedsigma2hat=VarCorr(fit)$group[1,1] observedsigma2hat 0.13 mean(sigma2star>observedsigma2hat) 0.23