Re: [R] Simulate p-value in lme4

From: Michael Kubovy <kubovy_at_virginia.edu>
Date: Sun 08 Oct 2006 - 11:34:07 GMT


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/

R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun Oct 08 21:40:56 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Sun 08 Oct 2006 - 18:30:09 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.