[R-sig-ME] p values for models with random correlation parameters

[Ben Bolker]

Babila Tachu wrote:
> Hello R programmers!
> 
> I used the package languageR for mixed effect models with the syntax
> at the end of this posting. I can use pvals.fnc to get p-values for
> model 1 and 3 (hd_lmer1 and hd_lmer3). Using this with model 2 gives
> the following error message:
> 
> p2 = pvals.fnc(hd_lmer2) Error in pvals.fnc(hd_lmer2) : MCMC sampling
> is not yet implemented in lme4_0.999375 for models with random
> correlation parameters
> 
> I'd be grateful if any one could help me out on how to get p-values
> for such models.
> 
> Models:
> 
> hd_lmer1 <- lmer(rot~ time + group + sex + gen + (1 | subject) +
> (1|rot.pre), data = data_long,REML = TRUE) hd_lmer2 <- lmer(rot~ time
> + group + sex + gen + (time | subject) + (1|rot.pre), data =
> data_long,REML = TRUE) hd_lmer3 <- lmer(rot~ time*group + sex + gen +
> (1 | subject) + (1|rot.pre), data =data_long,REML = TRUE)
> 
> rot = time spent by mice on an apparatus; gen= genotype
> 

  *If* the simulate method is working for these models (I'm not sure,
you haven't provided a reproducible example, and I'm too lazy to make
one up right now), you can get a p-value by simulating from the null
(reduced) model and fitting both the reduced and full models to the
simulated data.  For examples, to test the effect of the random effect
of time by subject:

devdiff <- numeric(1000)
newdata <- data_long
for (i in 1:1000) {
   newdata$rot <- simulate(hd_lmer3)  ## sim from model without time
   fit1 <- update(hd_lmer2,data=newdata)
   fit0 <- update(hd_lmer3,data=newdata)
   devdiff[i] <- -2*(logLik(fit1)-logLik(fit0))
}
mean(devdiff>=(logLik(hd_lmer2)-logLik(hd_lmer3))

 or something like that.

  warning, this may take a while ...

  Ben Bolker