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

[espesser]

Ben Bolker a écrit :
> 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
>
>   
You can also run a variant of hd_lmer2 without the correlation between 
the random intercept for subject
and the random slope for subject. You can test if this correlation term 
is needed with:
anova(hd_lmer2, hd_lmer2_bis)

hd_lmer2_bis <- lmer(rot~ time
+ group + sex + gen + (1|subject) +(0+ time | subject) + (1|rot.pre), data =
data_long,REML = TRUE) 

On this model, pvals.fnc would run.