[R] Simulate p-value in lme4

[Michael Kubovy]

Dear r-helpers,

Spencer Graves and Manual Morales proposed the following methods to  
simulate p-values in lme4:

************preliminary************

require(lme4)
require(MASS)
summary(glm(y ~ lbase*trt + lage + V4, family = poisson, data =  
epil), cor = FALSE)
epil2 <- epil[epil$period == 1, ]
epil2["period"] <- rep(0, 59); epil2["y"] <- epil2["base"]
epil["time"] <- 1; epil2["time"] <- 4
epil2 <- rbind(epil, epil2)
epil2$pred <- unclass(epil2$trt) * (epil2$period > 0)
epil2$subject <- factor(epil2$subject)
epil3 <- aggregate(epil2, list(epil2$subject, epil2$period > 0),  
function(x) if(is.numeric(x)) sum(x) else x[1])

************simulation (SG)************
o <- with(epil3, order(subject, y))
epil3. <- epil3[o,]
norep <- with(epil3., subject[-1]!=subject[-dim(epil3)[1]])
subj1 <- which(c(T, norep))
subj.pred <- epil3.[subj1, c("subject", "pred")]
subj. <- as.character(subj.pred$subject)
pred. <- subj.pred$pred
names(pred.) <- subj.

iter <- 10
chisq.sim <- rep(NA, iter)

set.seed(1)
for(i in 1:iter){
   s.i <- sample(subj.)
# Randomize subject assignments to 'pred' groups
   epil3.$pred <- pred.[s.i][epil3.$subject]
   fit1 <- lmer(y ~ pred+(1 | subject),
                 family = poisson, data = epil3.)
   fit0 <- lmer(y ~ 1+(1 | subject),
                 family = poisson, data = epil3.)
   chisq.sim[i] <- anova(fit0, fit1)[2, "Chisq"]
}

************simulation (MM)************
iter <- 10
chisq.sim <- rep(NA, iter)

Zt <- slot(model1,"Zt") # see ?lmer-class
n.grps <- dim(ranef(model1)[[1]])[1]
sd.ran.effs <- sqrt(VarCorr(model1)[[1]][1])
X <- slot(model1,"X") # see ?lmer-class
fix.effs <- matrix(rep(fixef(model1),dim(X)[1]), nrow=dim(X)[1],
                    byrow=T)
model.parms <- X*fix.effs # This gives parameters for each case
# Generate predicted values
pred.vals <- as.vector(apply(model.parms, 1, sum))

for(i in 1:iter) {
   rand.new <- as.vector(rnorm(grps,0, sd.ran.effs)) #grps should be  
n.grps?
   rand.vals <- as.vector(rand.new%*%Zt) # Assign random effects
   mu <- pred.vals+rand.vals # Expected mean
   resp <- rpois(length(mu), exp(mu))
   sim.data <- data.frame(slot(model2,"frame"), resp) # Make data frame
   sim.model1 <- lmer(resp~1+(1|subject), data=sim.data,
                      family="poisson")
   sim.model2 <- lmer(resp~pred+(1|subject), data=sim.data,
                      family="poisson")
   chisq.sim[i] <- anova(sim.model1,sim.model2)$Chisq[[2]]
}

************end************

Unfortunately the latter fails (even if I replace grps with n.grps):  
"Error in slot(model2, "frame") : object "model2" not found"

In any event, I would be eager to hear more discussion of the pros  
and cons of these approaches.

_____________________________
Professor Michael Kubovy
University of Virginia
Department of Psychology
USPS:     P.O.Box 400400    Charlottesville, VA 22904-4400
Parcels:    Room 102        Gilmer Hall
         McCormick Road    Charlottesville, VA 22903
Office:    B011    +1-434-982-4729
Lab:        B019    +1-434-982-4751
Fax:        +1-434-982-4766
WWW:    http://www.people.virginia.edu/~mk9y/