Skip to content

Different versions of lme4 and covariance of random effects

5 messages · Doran, Harold, David Afshartous

#
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)

##################################################################
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
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
#####################################################################
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
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
#
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?
#
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:

            
#
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
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 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
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
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
#
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: