[R-sig-ME] Different versions of lme4 and covariance of randomeffects
David Afshartous
dafshartous at med.miami.edu
Fri Aug 8 17:40:14 CEST 2008
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