5
$\begingroup$

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?

$\endgroup$
7
$\begingroup$

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".

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – gjabel Feb 19 '14 at 9:47
  • $\begingroup$ 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. $\endgroup$ – Wolfgang Feb 19 '14 at 14:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.