lme specifying correlation structure for three-level model

Question

I'm new to R and to multilevel modeling. I have a data set where I have a dependent variable y and predictor x, both of which are measured one time per day over a number of days within subjects. In addition, each subject is part of a twin pair. So in terms of my nesting structure, I have my measurements, nested within subject, nested within family (FamID). I expect my measurement values to be correlated over days within subjects, so I would like to specify an autocorrelation structure of order 1. Below is how I am specifying my model:

m1 <- lme(y ~ x, 
   random = list(~1 + Subject | FamID, ~1 + x | Subject), 
   data = dataset, method="ML", 
   correlation=corAR1(,form=~1|Subject),
   na.action="na.omit")

However, I receive the error message,

incompatible formulas for groups in 'random' and 'correlation'

Might anyone be able to help me appropriately specify this model?

Answer 1

As @Roland says, the correlation and random effects formulas must match. See below ...

Simulate data:

dd <- expand.grid(
    FamID = 1:3,
    Subject = 1:2,
    time = 1:10)
set.seed(101)
dd$x <- rnorm(nrow(dd))
dd$y <- rnorm(nrow(dd))

library(nlme)
m1 <- lme(y~x,
          random = ~1|FamID/Subject,
          data = dd,
          method = "ML",
          correlation = corAR1(form = ~1|FamID/Subject))

Notes:

• 1|FamID/Subject specifies "Subject nested within FamID", which sounds like what you described. Your current random effect specification list(~1 + Subject | FamID, ~1 + x | Subject) makes little sense to me: this would indicate
  • random effects of subject within family (i.e., separate variances for each subject, and an arbitrary correlation between subjects)
  • random slopes (effects of intercept and slope, arbitrarily correlated) within family

(the simpler 1|FamID/Subject specification does imply correlation between subjects, through the shared family effect; however, this correlation must be ≥ 0, unlike the 1 + Subject | FamID specification. The 1 + Subject | FamID` specification is also a little bit weird because it implies that the twins in a family are non-exchangeable, i.e. 'twin 1' and 'twin 2' would be specified in some way ...)

This is most likely overparameterized/unidentifiable (if you do want a random-slopes model allowing for the variation in the effects of x across subjects and/or families you can use random = ~1+x|FamID/Subject to estimate slopes at both levels — I checked, and this does still work with the correlation argument. I don't know if it's possible to specify random slopes at only one level (e.g. across subjects but not families) in lme ...

• corAR1(form = ~1|FamID/Subject) seems as though it might specify two autocorrelation parameters (at the levels of both family and subject-within-family), but according to the output (below) only one is estimated.

(Note that the random effects estimated are vanishingly small because I used made-up data with no structure.)

Linear mixed-effects model fit by maximum likelihood
  Data: dd 
  Log-likelihood: -81.77192
  Fixed: y ~ x 
(Intercept)           x 
 0.08731064 -0.09266083 

Random effects:
 Formula: ~1 | FamID
         (Intercept)
StdDev: 2.303609e-05

 Formula: ~1 | Subject %in% FamID
         (Intercept)  Residual
StdDev: 2.414397e-06 0.9598456

Correlation Structure: AR(1)
 Formula: ~1 | FamID/Subject 
 Parameter estimate(s):
      Phi 
-0.181599 
Number of Observations: 60
Number of Groups: 
             FamID Subject %in% FamID 
                 3                  6