# Re: [R] Simulate p-value in lme4

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Mon 21 Aug 2006 - 04:22:16 EST

Beyond this, you are to be commended for providing such a simple, self-contained example for such a sophisticated question. I think you simulation misses one important point: It assumes the between-subject variance is zero. To overcome this, I think I might try either the bootstrap or permutation distribution scrambling the assignment of subjects to treatment groups but otherwise keeping the pairs of observations together.

To this end, consider the following:

# Build a table to translate subject into 'pred' o <- with(epil3, order(subject, y))
epil3. <- epil3[o,]
norep <- with(epil3., subject[-1]!=subject[-dim(epil3)]) 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.

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"] }
```      Hope this helps.
Spencer Graves

```

Manuel Morales wrote:
> Dear list,
>
> This is more of a stats question than an R question per se. First, I
> realize there has been a lot of discussion about the problems with
> estimating P-values from F-ratios for mixed-effects models in lme4.
> Using mcmcsamp() seems like a great alternative for evaluating the
> significance of individual coefficients, but not for groups of
> coefficients as might occur in an experimental design with 3 treatment
> levels. I'm wondering if the simulation approach I use below to estimate
> the P-value for a 3-level factor is appropriate, or if there are any
> suggestions on how else to approach this problem. The model and data in
> the example are from section 10.4 of MASS.
>
> Thanks!
> Manuel
>
> # Load req. package (see functions to generate data at end of script)
> library(lme4)
> library(MASS)
>
> # Full and reduced models - pred is a factor with 3 levels
> result.full <- lmer(y~pred+(1|subject), data=epil3, family="poisson")
> result.base <- lmer(y~1+(1|subject), data=epil3, family="poisson")
>
> # Naive P-value from LR for significance of "pred" factor
> anova(result.base,result.full)\$"Pr(>Chisq)"[] # P-value
> (test.stat <- anova(result.base,result.full)\$Chisq[]) # Chisq-stat
>
> # P-value from simulation. Note that in the simulation, I use the
> # estimated random effects for each subject rather than generating a new
> # distribution of means. I'm not sure if this is appropriate or not ...
> intercept <- fixef(result.base)
> rand.effs <- ranef(result.base)[]
> mu <- exp(rep(intercept+rand.effs[],2))
>
> p.value <- function(iter, stat) {
> chi.stat <- vector()
> for(i in 1:iter) {
> resp <- rpois(length(mu), mu) # simulate values
> sim.data <- data.frame(y=resp,subject=epil3\$subject,pred=epil3\$pred)
> result.f <- lmer(y~pred+(1|subject), data=sim.data,
> family="poisson")
> result.b <- lmer(y~1+(1|subject), data=sim.data, family="poisson")
> chi.stat[i] <- anova(result.b,result.f)\$Chisq[]
> }
> val <- sum(unlist(lapply(chi.stat, function(x) if(x>stat) 1 else
> 0)))/iter
> hist(chi.stat)
> return(val)
> }
>
> p.value(10,test.stat) # Increase to >=1000 to get a reasonable P-value!
>
> # Script to generate data, from section 10.4 of MASS
> 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)
> epil3\$pred <- factor(epil3\$pred, labels = c("base", "placebo", "drug"))
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
>
https://stat.ethz.ch/mailman/listinfo/r-help