[R-sig-ME] Moving from lme4a to lme4Eigen correlation of random effects goes to 1

[Ben Bolker]

Jennifer Lyon <jennifer.s.lyon at ...> writes:

> I am modeling the score participants in three conditions achieved at
> three times.  The times are not equally spaced, with the difference
> between t2 and t3 over an order of magnitude longer than the
> difference between t1 and t2. The design is not balanced. The variable
> LMNo is unique per individual.
> 
>          time
> condition t1 t2 t3
>         A 76 75 72
>         B 18 18 17
>         C 28 27 26

> I've attached a plot of the data.

  (didn't come through -- the mailing list discards a lot
of file types)

> The literature suggests that there are individual differences
> at the different times, so I fitted the following model:
> 
> me.c.m<-lmer(score~ 1+time + condition + (1+time|LMNo), me.c, REML = 0)
> 
 When moving from lme4a to lme4Eigen, the correlations of the random
 effects all went to 1.000. I was slightly surprised by this change of
 events. Is this result indicating that the simpler model:
 
 me.c.m0<-lmer(score~ 1+time + condition + (1|LMNo), me.c, REML = 0)
 
 is sufficient? 

   It typically would.  It does mean you're somewhere on the
edge of overfitting ...
   If time were used as a continuous predictor (i.e. linear regression),
it might make sense to try to fit the model with (1|LMNo)+(0+time|LMNo),
but I suspect it *doesn't* make sense when time is a factor.

 I ask because when I do a likelihood ratio test
 using anova, the p-value is small and AIC prefers the more complex
 model while BIC prefers the simpler model. Does the correlation
 going to one also indicate a preference for the simpler model?
 
 #run in lme4Eigen
 anova(me.c.m0, me.c.m)
 Data: me.c
 Models:
 me.c.m0: score ~ 1 + time + condition + (1 | LMNo)
 me.c.m: score ~ 1 + time + condition + (1 + time | LMNo)
         Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(Chisq)
 me.c.m0  7 2069.8 2096.9 -1027.9   2055.8
 me.c.m  12 2055.3 2101.9 -1015.7   2031.3 24.445      5  0.0001782 ***

  On the other hand, this seems to say pretty definitively that
the variation in time effects across individuals is doing something ...

 I don't know if this is related, but when I run profile() on
 the models in lme4a and lme4Eigen, I get an error (which is
 shown below). I have additional information on the participants,
 such as gender and which participants are siblings, but before
 I go for a more complex model, I'd like to better understand
 what this model is telling me.

  Very wise.


> Here are the details of fitting the models in lme4a and lme4Eigen:
> 
> In lme4a_0.9996875-1
> 
> > library(lme4a)
> > me.c<-read.table("mixed-effects-data-clean.txt", header=T)
> 
> > me.c.m<-lmer(score~ 1+time + condition + (1+time|LMNo), me.c, REML = 0)
> > summary(me.c.m)

 [snip]

>    Data: me.c
>       AIC       BIC    logLik  deviance
>  2049.704  2096.236 -1012.852  2025.704
> 
> Random effects:
>  Groups   Name        Variance Std.Dev. Corr
>  LMNo     (Intercept) 41.280   6.425
>           timet2       6.007   2.451    0.452
>           timet3       5.174   2.275    0.382 0.997
>  Residual              4.314   2.077
> Number of obs: 357, groups: LMNo, 122
> 
> Fixed effects:

 [snip]
 
> > profile(me.c.m)
> Error in rbind2(fillmat(pres, lowcut, zeta, wp1), fillmat(nres, lowcut,  :
>   error in evaluating the argument 'y' in selecting a method for
> function 'rbind2': Error in while (i < nr && abs(mat[i, ".zeta"]) <=
> cutoff && mat[i, cc] >  :
>   missing value where TRUE/FALSE needed

  Hmmm.  Are you willing to send data?

  As cross-checks on lme4a and lme4Eigen, you might try
the packages
* glmmADMB (recent versions can handle Gaussian responses,
although it's not well tested)
* regress
* lmm
* sabreR

  although all (except glmmADMB) have fairly different
interfaces, unfortunately. 

> 

 [snip]
 
 > me.c<-read.table("mixed-effects-data-clean.txt", header=T)
 
 > me.c.m<-lmer(score~ 1+time + condition + (1+time|LMNo), me.c, REML = 0)
 > summary(me.c.m)
 Linear mixed model fit by maximum likelihood ['summary.mer']
 Formula: score ~ 1 + time + condition + (1 + time | LMNo)
    Data: me.c
 
       AIC       BIC    logLik  deviance
  2055.339  2101.872 -1015.669  2031.339
 
 Random effects:
  Groups   Name        Variance Std.Dev. Corr
  LMNo     (Intercept) 39.351   6.273
           timet2       2.440   1.562    1.000
           timet3       1.730   1.315    1.000 1.000
  Residual              5.507   2.347
 Number of obs: 357, groups: LMNo, 122
 
 Fixed effects:

 [snip]
 
 > profile(me.c.m)
 Warning message:
 In sqrt(ores$fval - base) : NaNs produced
 Error in rbind2(fillmat(pres, lowcut, zeta, wp1), fillmat(nres, lowcut,  :
   error in evaluating the argument 'x' in selecting a method for
 function 'rbind2': Error in while (i < nr && abs(mat[i, ".zeta"]) <=
 cutoff && mat[i, cc] >  :
   missing value where TRUE/FALSE needed
 

[snip]