# Dropping term for correlation between random effects in lme and interpretting summary output

I want to fit a model without a correlation term between the random effects with lme. In lmer this is fairly straighforward....

# lmer without correlation term
m1 <- lmer(distance ~ (1|Subject) + age + (0+age|Subject) + Sex, data = Orthodont)
VarCorr(m1)
# Groups    Name        Std.Dev.
# Subject   (Intercept) 1.474105
# Subject.1 age         0.099979
# Residual              1.402591


With lme I think I can drop the correlation term using the following specification...

# lme without correlation term?
m2 <- lme(distance ~ age + Sex, data = Orthodont, random = list(~ 1 | Subject, ~-1+ age | Subject))
VarCorr(m2)
#             Variance            StdDev
# Subject =   pdLogChol(1)
# (Intercept) 2.172946296         1.47409169
# Subject =   pdLogChol(-1 + age)
# age         0.009996006         0.09998003
# Residual    1.967260819         1.40259075


I am not entirely convinced that these are the same models though, partly because I can not find any resources that details how to specify this specific form and partly because the output from print is a little mystifying to me...

m2
# Linear mixed-effects model fit by REML
#   Data: Orthodont
#   Log-restricted-likelihood: -218.3227
#   Fixed: distance ~ age + Sex
# (Intercept)         age   SexFemale
#  17.5806928   0.6601852  -2.0117005
#
# Random effects:
#  Formula: ~1 | Subject
#         (Intercept)
# StdDev:    1.474092
#
#  Formula: ~-1 + age | Subject %in% Subject
#                age Residual
# StdDev: 0.09998003 1.402591
#
# Number of Observations: 108
# Number of Groups:
#                Subject Subject.1 %in% Subject
#                     27                     27


In particular, what is Subject %in% Subject referring to? and why are the residual considered as part of the second random effect term?

While in principle your approach works, this is not quite the 'standard' way of making the random intercepts and slopes uncorrelated. With lme you can use pdClasses (see help(pdClasses)) to give a particular structure to the variance-covariance matrix of the random effects. Here, you want to make that matrix diagonal. You can do that with:

m3 <- lme(distance ~ age + Sex, data = Orthodont, random = list(Subject = pdDiag(~ age))
m3

Linear mixed-effects model fit by REML
Data: Orthodont
Log-restricted-likelihood: -218.3227
Fixed: distance ~ age + Sex
(Intercept)         age   SexFemale
17.5806928   0.6601852  -2.0117005

Random effects:
Formula: ~age | Subject
Structure: Diagonal
(Intercept)        age Residual
StdDev:    1.474092 0.09998003 1.402591

Number of Observations: 108
Number of Groups: 27


You will find that the parameter estimates are actually identical to model m2, but the presentation of the results is more "logical".

• Thanks, that is great. I have been exploring a little, and I think my m2 has an extra hierarchy of Subject nested within Subject, hence the unexpected output. – gjabel Feb 19 '14 at 9:47
• Yes, that's what Subject %in% Subject means. But this obviously makes little sense. The results work out to be the same, but really using pdClasses is much more sensible here. – Wolfgang Feb 19 '14 at 14:48