How to correctly specify the within-group correlation matrix for linear mixed model in R?

Question

I'm using the nlme package's lme function in R to fit a random-intercept, random-slope linear mixed model for some generated test data. Although the fixed effect coefficients are estimated as expected, the variance parameter estimation yields results I do not fully understand.

Specifically, the rho parameter of the compound symmetry correlation structure specified in the model statement and correlation of the random effects are not estimated, as I would expect from the data. When not fixing the correlation in corCompSymm(), the estimated rho is almost equal to default 0, and correlation of random effects is then obviously wrong (0.64 instead of -0.5).

I'm aware, that the variance-covariance parameters are estimated marginally, but is there a way to improve the estimates, given I correctly specify the correlation matrix type (CS, AR, etc.) but do not know it's parameters?

And secondly (somehow unrelated), does R allow for fitting a linear model with no random effects specified, but assumed structure for within-subject correlation matrix (analogous to SAS's PROC MIXED with no RANDOM statement)

---

The fitted model is:

lme(data = dat, fixed = y ~ as.factor(f) * x, random = ~ 1 + x | id, correlation = corCompSymm(fixed = FALSE))

and the outcome variable is simulated as:

y <- ifelse(f == 0, beta0[1], beta0[2]) + ifelse(f == 0, beta1[1] * x, beta1[2] * x + b1 * x) + b0 + b1 * x + e, where:

• f is a subject-specific grouping factor with 2 levels
• x simulates ordinal variable over which the repeated measures have been taken (design is balanced)
• beta0 and beta1 are vectors of group-specific intercept and slope, consecutively
• e is an error-term, following multivariate normal with means 0 and compound symmetry variance-covariance matrix, having 1 diagonal and .8 off-diagonal
• b0 and b1 are subject-specific errors for intercept, and slope consecutively, simulated to follow bivariate standard normal with covariance -0.5:

\begin{bmatrix} 1 & - 0.5\\ - 0.5 & 1 \end{bmatrix}

The estimated parameters:

---

Correlation Structure: Compound symmetry Formula: ~1 | id Parameter estimate(s): Rho 4.260007e-06

> VarCorr(f3) id = pdLogChol(1 + x) Variance StdDev Corr (Intercept) 0.5682353 0.7538139 (Intr) x 2.5251728 1.5890792 0.64 Residual 5.0764053 2.2530879

---

Reproducible code:

require(reshape2)
require(nlme)
require(MASS)

set.seed(1)

n <- 1000 # number of subjects
m <- 4 # rep. measurments per subject

beta0 <- c(2, 5) # intercepts for each subject-specific factor
beta1 <- c(2, 3) # slopes for each subject-specific factor

# correlation matrix for random effects
d <- - .50
D <- matrix(c(1, d,
              d, 1), nrow = 2, byrow = T)

# correlation matrix for error terms
r <- .8
R <- c(1, r, r, r,
       r, 1, r, r,
       r, r, 1, r,
       r, r, r, 1) * 5
R <- matrix(R, nrow = sqrt(length(R)))
R <- R[1:m, 1:m]

dat <- data.frame(id = 1:n)

dat$f <- sample(1:length(beta0) - 1, n, replace = TRUE) # randomly assign a factor for each subject

# assign subject-specific random slopes and intercepts
dat <- cbind(dat, 
             setNames(data.frame(
               mvrnorm(n, mu = rep(0, nrow(D)), Sigma = D)), 
               c('b0', 'b1'))) 

dat <- dat[rep(1:n, each = m), ]
dat$x <- rep(1:m, n)

# error term
dat$e <- as.vector(t(mvrnorm(n, mu = rep(0, nrow(R)), Sigma = R))) # 

# generate outcome
dat$y <- with(dat, 
              ifelse(f == 0, beta0[1], beta0[2]) +
                ifelse(f == 0, beta1[1] * x, beta1[2] * x + b1 * x) +
                b0 + b1 * x + e)

# lmm with assumed compound symmetry and correlation not fixed
summary(f3 <- lme(data = dat, fixed = y ~ as.factor(f) * x, random = ~ 1 + x | id, 
                  correlation = corCompSymm(fixed = FALSE)))
VarCorr(f3)

# lmm with with fixed correlation
summary(f4 <- lme(data = dat, fixed = y ~ as.factor(f) * x, random = ~ 1 + x | id, 
                  correlation = corCompSymm(r, fixed = TRUE)))
VarCorr(f4)    

Answer 1

The answer to your second question (fit correlation structure without random effects) is to use nlme::gls() ("generalized least squares") - it allows the same set of heteroscedasticity (weights argument) and correlation (correlation argument) as lme.

Answer 2

I think the problem is that you are interpreting the correlation argument wrong. This argument allows you to specify the within-subjects correlation, not the correlation of the random effects. The correlation argument is useful when there is residual serial correlation that the random intercept model did not sufficiently not account for.

The semi-variogram of longitudinal data helps with this interpretation, which is described in Analysis of Longitudinal Data and this highly detailed paper. With a sufficiently complex within-subject correlation structure, we can fit a random intercept to account for between-person variance, and then also include a model of the serial correlation (e.g. a Gaussian correlation structure to account for the decay in the within-subject correlation at greater lags).

In my experience, one can either model the structure of the semi-variogram using both the random andcorrelation` arguments, or one can use a random intercept and slope model. I think of it as two different ways to model the same phenomenon. The first case models the correlation structure directly, whereas the latter conditions out the correlation structure by fitting subject-specific curves.

So in the end, it makes sense to me that the model estimates a near-0 serial correlation for the corCompSymm structure because after including the random intercept and slopes, there is no serial correlation remaining to account for.