[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