[R] Simulate p-value in lme4

[Manuel Morales]

On Sun, 2006-10-08 at 07:34 -0400, Michael Kubovy wrote:
> 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?

Yes. Change grps to 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"

You need to specify the models fit0 and fit1 from SG's example as model1
and model2. E.g., model1 <- fit0, etc.

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

One practical problem with my approach (MM) is that the fitting
algorithms for lmer would often choke  - in particular for those
randomly generated data sets that by chance included variance components
close to zero.

In any case, the MCMC approach advocated by Douglas Bates is *much*
faster. That's the approach I've been using since he suggested a
function for estimating P-values from MCMC samples for factors (or model
comparisons) with greater than 2 levels.

See:http://tolstoy.newcastle.edu.au/R/e2/help/06/09/1003.html

and the very long thread that accompanies it.

HTH,

Manuel

-- 
Manuel A. Morales
http://mutualism.williams.edu
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: This is a digitally signed message part
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20061008/e98860ed/attachment.bin