[R-sig-ME] Different versions of lme4 and covariance of randomeffects
David Afshartous
dafshartous at med.miami.edu
Fri Aug 8 18:13:49 CEST 2008
I just installed the latest version of Matrix in the 2.6.2 session and also
get the same results as well. So it appears that the version of Matrix was
the issue. However, now that I have the newer version of Matrix how does
one re-produce the results that the older version was producing?
On 8/8/08 11:54 AM, "Doran, Harold" <HDoran at air.org> wrote:
> 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