[R-sig-ME] Moving from lme4a to lme4Eigen correlation of random effects goes to 1
Jennifer Lyon
jennifer.s.lyon at gmail.com
Sat Feb 18 18:44:48 CET 2012
Hi:
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.
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? 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 ***
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.
Thanks for any and all insights.
Jen
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)
Linear mixed model fit by maximum likelihood ['summary.mer']
Formula: score ~ 1 + time + condition + (1 + time | LMNo)
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:
Estimate Std. Error t value
(Intercept) 27.8800 0.7727 36.08
timet2 -3.6805 0.3492 -10.54
timet3 -4.2801 0.3429 -12.48
conditionB -1.9732 1.7585 -1.12
conditionC -10.0301 1.4833 -6.76
Correlation of Fixed Effects:
(Intr) timet2 timet3 cndtnB
timet2 0.086
timet3 0.039 0.687
conditionB -0.436 0.000 0.000
conditionC -0.517 -0.001 0.000 0.227
> 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
> sessionInfo()
R version 2.14.1 (2011-12-22)
Platform: x86_64-unknown-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=C LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4a_0.9996875-1 MatrixModels_0.3-1 minqa_1.1.18 Rcpp_0.9.10
[5] Matrix_1.0-3 lattice_0.20-0
loaded via a namespace (and not attached):
[1] codetools_0.2-8 colorspace_1.1-1 grid_2.14.1 int64_1.1.2
[5] nlme_3.1-103 splines_2.14.1
> q()
R --vanilla
> library(lme4Eigen)
> 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:
Estimate Std. Error t value
(Intercept) 27.7704 0.7639 36.35
timet2 -3.6838 0.3348 -11.00
timet3 -4.2865 0.3305 -12.97
conditionB -1.9870 1.7285 -1.15
conditionC -9.5437 1.4583 -6.54
Correlation of Fixed Effects:
(Intr) timet2 timet3 cndtnB
timet2 0.138
timet3 0.089 0.569
conditionB -0.433 0.001 0.000
conditionC -0.514 -0.002 -0.001 0.227
> 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
> sessionInfo()
R version 2.14.1 (2011-12-22)
Platform: x86_64-unknown-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=C LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4Eigen_0.9996875-9 lattice_0.20-0
loaded via a namespace (and not attached):
[1] colorspace_1.1-1 grid_2.14.1 Matrix_1.0-3 minqa_1.1.18
[5] nlme_3.1-103 Rcpp_0.9.10 splines_2.14.1
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mixed-model-participant-score-by-condition.pdf
Type: application/pdf
Size: 6002 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20120218/b788c8a3/attachment-0001.pdf>
More information about the R-sig-mixed-models
mailing list