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