[R-sig-ME] Different versions of lme4 and covariance of randomeffects
Doran, Harold
HDoran at air.org
Fri Aug 8 17:54:07 CEST 2008
Well, I don't get this problem. Using your code, I get the same exact
output from 2.6.2 and 2.7.1. My sessionInfo shows that I have the same
packages as you in 2.7.1, but in 2.6.2 you have an older version of
Matrix. My output showing the same estimates from both R versions is
pasted below.
## output from 2.7.1
> summary(fm.het.3)
Linear mixed model fit by REML
Formula: dv ~ time.num * drug - 1 + (0 + Dind + Pind | Patient.cross)
Data: dat.new
AIC BIC logLik deviance REMLdev
705 729.7 -344.5 687.6 689
Random effects:
Groups Name Variance Std.Dev. Corr
Patient.cross Dind 29.3378 5.4164
Pind 4.8857 2.2104 0.169
Residual 1.9651 1.4018
Number of obs: 160, groups: Patient.cross, 20
Fixed effects:
Estimate Std. Error t value
time.num 0.8175 0.1402 5.832
drugD -0.8505 1.2705 -0.669
drugP 5.9720 0.6258 9.543
time.num:drugP -1.1262 0.1982 -5.681
Correlation of Fixed Effects:
tim.nm drugD drugP
drugD -0.276
drugP 0.000 0.127
tim.nm:drgP -0.707 0.195 -0.396
> sessionInfo()
R version 2.7.1 (2008-06-23)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_0.999375-24 Matrix_0.999375-11 lattice_0.17-13
loaded via a namespace (and not attached):
[1] grid_2.7.1 tools_2.7.1
### Output from 2.6.2
> summary(fm.het.3)
Linear mixed-effects model fit by REML
Formula: dv ~ time.num * drug - 1 + (0 + Dind + Pind | Patient.cross)
Data: dat.new
AIC BIC logLik MLdeviance REMLdeviance
703 724.6 -344.5 687.6 689
Random effects:
Groups Name Variance Std.Dev. Corr
Patient.cross Dind 29.3339 5.4161
Pind 4.8854 2.2103 0.169
Residual 1.9651 1.4018
number of obs: 160, groups: Patient.cross, 20
Fixed effects:
Estimate Std. Error t value
time.num 0.8175 0.1402 5.832
drugD -0.8505 1.2705 -0.669
drugP 5.9720 0.6258 9.543
time.num:drugP -1.1262 0.1982 -5.681
Correlation of Fixed Effects:
tim.nm drugD drugP
drugD -0.276
drugP 0.000 0.127
tim.nm:drgP -0.707 0.195 -0.396
> sessionInfo()
R version 2.6.2 (2008-02-08)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_0.99875-9 Matrix_0.999375-7 lattice_0.17-4
loaded via a namespace (and not attached):
[1] grid_2.6.2
> -----Original Message-----
> From: David Afshartous [mailto:dafshartous at med.miami.edu]
> Sent: Friday, August 08, 2008 11:40 AM
> To: Doran, Harold; r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] Different versions of lme4 and
> covariance of randomeffects
>
>
> Donald,
>
> Sorry, you're right, for the reproducible example I was
> mainly focusing on the correlation term but when I look
> closer there are other differences as well. And yes, the
> same data was used, generated via the code below with same
> seed at set.seed(500).
>
> Cheers,
> David
>
>
>
>
> On 8/8/08 11:30 AM, "Doran, Harold" <HDoran at air.org> wrote:
>
> > David:
> >
> > I'm looking at this and am a little confused. You say the
> differences
> > are not too dramatic. But, they are. I can see that the syntax for
> > your model specification is the same, but the estimates of the fit
> > statistics differ (logLik), the estimates of the fixed
> effects differ,
> > as well as their standard errors, and the variance
> components differ.
> >
> > Am I missing something? Was the same data set used in both cases?
> >
> >> -----Original Message-----
> >> From: r-sig-mixed-models-bounces at r-project.org
> >> [mailto:r-sig-mixed-models-bounces at r-project.org] On
> Behalf Of David
> >> Afshartous
> >> Sent: Friday, August 08, 2008 11:06 AM
> >> To: r-sig-mixed-models at r-project.org
> >> Subject: [R-sig-ME] Different versions of lme4 and covariance of
> >> randomeffects
> >>
> >>
> >>
> >> All,
> >>
> >> I recently re-estimated a model after upgrading to Rv2.7.1
> that had
> >> been estimated Rv2.6.2 previously and was surprised to see the
> >> estimated correlation between random effects in the model
> change from
> >> 0.3 to 1.0.
> >> It appears that the reason has more to do with the version of
> >> lme4 since when I download the latest possible lme4 to
> >> Rv2.6.2 the results agree.
> >>
> >> Below is a reproducible example where the difference in results is
> >> not as dramatic. Of course, for my data the initial
> results with the
> >> correlation of .3 seem a lot more plausible than that of
> 1.0 so I'm
> >> inclined to trust those results more than the newer ones.
> It seems
> >> strange to get the correlation of 1.
> >>
> >>
> >> Cheers,
> >> David
> >>
> >> library("lme4")
> >> set.seed(500)
> >> n.timepoints <- 4 ## change for shorter examples
> n.subj.per.tx <- 20
> >> sd.d <- 5; sd.p <- 2; sd.res <- 1.3 drug
> >> <- factor(rep(c("D", "P"), each = n.timepoints, times =
> >> n.subj.per.tx))
> >> drug.baseline <- rep( c(0,5), each=n.timepoints,
> times=n.subj.per.tx
> >> ) Patient <- rep(1:(n.subj.per.tx*2), each = n.timepoints)
> >> Patient.baseline <- rep( rnorm( n.subj.per.tx*2,
> sd=c(sd.d, sd.p) ),
> >> each=n.timepoints ) time
> >> <- factor(paste("Time-", rep(1:n.timepoints, n.subj.per.tx*2),
> >> sep=""))
> >> time.baseline <-
> >> rep(1:n.timepoints,n.subj.per.tx*2)*as.numeric(drug=="D")
> >> dv <- rnorm( n.subj.per.tx*n.timepoints*2,
> >> mean=time.baseline+Patient.baseline+drug.baseline, sd=sd.res
> >> ) dat.new <- data.frame(time, drug, dv, Patient)
> >> dat.new$Patient.cross <- rep(1:(n.subj.per.tx), each =
> >> 2*n.timepoints)
> >>
> >>
> >> dat.new$Dind <- as.numeric(dat.new$drug == "D") dat.new$Pind
> >> <- as.numeric(dat.new$drug == "P") dat.new$time.num =
> >> rep(1:n.timepoints, n.subj.per.tx*2)
> >>
> >> ##################################################################
> >>
> >>> sessionInfo()
> >> R version 2.6.2 (2008-02-08)
> >> i386-pc-mingw32
> >>
> >> locale:
> >> LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
> >> States.1252;LC_MONETARY=English_United
> >> States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
> >>
> >> attached base packages:
> >> [1] stats graphics grDevices utils datasets
> methods base
> >>
> >> other attached packages:
> >> [1] lme4_0.99875-9 Matrix_0.999375-5 lattice_0.17-4
> >>
> >> loaded via a namespace (and not attached):
> >> [1] grid_2.6.2
> >>> ( fm.het.3 <- lmer( dv ~ time.num*drug -1 + ( 0 + Dind + Pind |
> >> Patient.cross ), data=dat.new ) )
> >> Linear mixed-effects model fit by REML
> >> Formula: dv ~ time.num * drug - 1 + (0 + Dind + Pind |
> Patient.cross)
> >> Data: dat.new
> >> AIC BIC logLik MLdeviance REMLdeviance
> >> 684.7 706.3 -335.4 669.2 670.7
> >> Random effects:
> >> Groups Name Variance Std.Dev. Corr
> >> Patient.cross Dind 26.8380 5.1805
> >> Pind 7.7623 2.7861 0.060
> >> Residual 1.5906 1.2612
> >> number of obs: 160, groups: Patient.cross, 20
> >>
> >> Fixed effects:
> >> Estimate Std. Error t value
> >> time.num 0.9101 0.1261 7.216
> >> drugD -0.2637 1.2088 -0.218
> >> drugP 5.0330 0.7123 7.066
> >> time.num:drugP -0.9476 0.1784 -5.313
> >>
> >> Correlation of Fixed Effects:
> >> tim.nm drugD drugP
> >> drugD -0.261
> >> drugP 0.000 0.050
> >> tim.nm:drgP -0.707 0.184 -0.313
> >>>
> >>
> >>
> #####################################################################
> >>
> >>> sessionInfo()
> >> R version 2.7.1 (2008-06-23)
> >> i386-pc-mingw32
> >>
> >> locale:
> >> LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
> >> States.1252;LC_MONETARY=English_United
> >> States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
> >>
> >> attached base packages:
> >> [1] stats graphics grDevices utils datasets
> methods base
> >>
> >> other attached packages:
> >> [1] lme4_0.999375-24 Matrix_0.999375-11 lattice_0.17-8
> >>
> >> loaded via a namespace (and not attached):
> >> [1] grid_2.7.1
> >>> ( fm.het.3 <- lmer( dv ~ time.num*drug -1 + ( 0 + Dind + Pind |
> >> Patient.cross ), data=dat.new ) )
> >> Linear mixed model fit by REML
> >> Formula: dv ~ time.num * drug - 1 + (0 + Dind + Pind |
> Patient.cross)
> >> Data: dat.new
> >> AIC BIC logLik deviance REMLdev
> >> 705 729.7 -344.5 687.6 689
> >> Random effects:
> >> Groups Name Variance Std.Dev. Corr
> >> Patient.cross Dind 29.3378 5.4164
> >> Pind 4.8857 2.2104 0.169
> >> Residual 1.9651 1.4018
> >> Number of obs: 160, groups: Patient.cross, 20
> >>
> >> Fixed effects:
> >> Estimate Std. Error t value
> >> time.num 0.8175 0.1402 5.832
> >> drugD -0.8505 1.2705 -0.669
> >> drugP 5.9720 0.6258 9.543
> >> time.num:drugP -1.1262 0.1982 -5.681
> >>
> >> Correlation of Fixed Effects:
> >> tim.nm drugD drugP
> >> drugD -0.276
> >> drugP 0.000 0.127
> >> tim.nm:drgP -0.707 0.195 -0.396
> >>>
> >>
> >> _______________________________________________
> >> R-sig-mixed-models at r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >>
>
>
More information about the R-sig-mixed-models
mailing list