[R-sig-ME] Different versions of lme4 and covariance of randomeffects

[Doran, Harold]

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
>