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

lme4-nlmemixed modelmultivariate analysis

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)

Best Answer

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.

Related Solutions

ANOVA – What is Compound Symmetry Explained in Plain English?

Compound symmetry is essentially the "exchangeable" correlation structure, except with a specific decomposition for the total variance. For example, if you have mixed model for the subject $i$ in cluster $j$ response, $Y_{ij}$, with only a random intercept by cluster

$$ Y_{ij} = \alpha + \gamma_{j} + \varepsilon_{ij} $$

where $\gamma_{j}$ is the cluster $j$ random effect with variance $\sigma^{2}_{\gamma}$ and $\varepsilon_{ij}$ is the subject $i$ in cluster $j$ "measurement error" with variance $\sigma^{2}_{\varepsilon}$ and $\gamma_{j}, \varepsilon_{ij}$ are independent. This model implicitly specifies the compound symmetry covariance matrix between observations in the same cluster:

$$ {\rm cov}(Y_{ij}, Y_{kj}) = \sigma^{2}_{\gamma} + \sigma^{2}_{\varepsilon} \cdot \mathcal{I}(k = i) $$

Note that the compound symmetry assumption implies that the correlation between distinct members of a cluster is $\sigma^{2}_{\gamma}/(\sigma^{2}_{\gamma} + \sigma^{2}_{\varepsilon})$.

In "plain english" you might say this covariance structure implies that all distinct members of a cluster are equally correlated with each other and the total variation, $\sigma^{2} = \sigma^{2}_{\gamma} + \sigma^{2}_{\varepsilon}$, can be partitioned into the "shared" (within a cluster) component, $\sigma^{2}_{\gamma}$ and the "unshared" component, $\sigma^{2}_{\varepsilon}$.

Edit: To aid understanding in the "plain english" sense, consider an example where individuals are clustered within families so that $Y_{ij}$ denotes the subject $i$ in family $j$ response. In this case the compound symmetry assumption means that the total variation in $Y_{ij}$ can be partitioned into the variation within a family, $\sigma^{2}_{\varepsilon}$, and the variation between families, $\sigma^{2}_{\gamma}$.

Paired t-test – Special Case of Linear Mixed-Effect Modeling

The equivalence of the models can be observed by calculating the correlation between two observations from the same individual, as follows:

As in your notation, let $Y_{ij} = \mu + \alpha_i + \beta_j + \epsilon_{ij}$, where $\beta_j \sim N(0, \sigma_p^2)$ and $\epsilon_{ij} \sim N(0, \sigma^2)$. Then $Cov(y_{ik}, y_{jk}) = Cov(\mu + \alpha_i + \beta_k + \epsilon_{ik}, \mu + \alpha_j + \beta_k + \epsilon_{jk}) = Cov(\beta_k, \beta_k) = \sigma_p^2$, because all other terms are independent or fixed, and $Var(y_{ik}) = Var(y_{jk}) = \sigma_p^2 + \sigma^2$, so the correlation is $\sigma_p^2/(\sigma_p^2 + \sigma^2)$.

Note that the models however are not quite equivalent as the random effect model forces the correlation to be positive. The CS model and the t-test/anova model do not.

EDIT: There are two other differences as well. First, the CS and random effect models assume normality for the random effect, but the t-test/anova model does not. Secondly, the CS and random effect models are fit using maximum likelihood, while the anova is fit using mean squares; when everything is balanced they will agree, but not necessarily in more complex situations. Finally, I'd be wary of using F/df/p values from the various fits as measures of how much the models agree; see Doug Bates's famous screed on df's for more details. (END EDIT)

The problem with your R code is that you're not specifying the correlation structure properly. You need to use gls with the corCompSymm correlation structure.

Generate data so that there is a subject effect:

set.seed(5)
x <- rnorm(10)
x1<-x+rnorm(10)
x2<-x+1 + rnorm(10)
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), 
                    rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

Then here's how you'd fit the random effects and the compound symmetry models.

library(nlme)
fm1 <- lme(y ~ x, random=~1|subj, data=myDat)
fm2 <- gls(y ~ x, correlation=corCompSymm(form=~1|subj), data=myDat)

The standard errors from the random effects model are:

m1.varp <- 0.5453527^2
m1.vare <- 1.084408^2

And the correlation and residual variance from the CS model is:

m2.rho <- 0.2018595
m2.var <- 1.213816^2

And they're equal to what is expected:

> m1.varp/(m1.varp+m1.vare)
[1] 0.2018594
> sqrt(m1.varp + m1.vare)
[1] 1.213816

Other correlation structures are usually not fit with random effects but simply by specifying the desired structure; one common exception is the AR(1) + random effect model, which has a random effect and AR(1) correlation between observations on the same random effect.

EDIT2: When I fit the three options, I get exactly the same results except that gls doesn't try to guess the df for the term of interest.

> summary(fm1)
...
Fixed effects: y ~ x 
                 Value Std.Error DF   t-value p-value
(Intercept) -0.5611156 0.3838423  9 -1.461839  0.1778
xx2          2.0772757 0.4849618  9  4.283380  0.0020

> summary(fm2)
...
                 Value Std.Error   t-value p-value
(Intercept) -0.5611156 0.3838423 -1.461839  0.1610
xx2          2.0772757 0.4849618  4.283380  0.0004

> m1 <- lm(y~ x + subj, data=myDat)
> summary(m1)
...
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.3154     0.8042  -0.392  0.70403   
xx2           2.0773     0.4850   4.283  0.00204 **

(The intercept is different here because with the default coding, it's not the mean of all subjects but instead the mean of the first subject.)

It's also of interest to note that the newer lme4 package gives the same results but doesn't even try to compute a p-value.

> mm1 <- lmer(y ~ x + (1|subj), data=myDat)
> summary(mm1)
...
            Estimate Std. Error t value
(Intercept)  -0.5611     0.3838  -1.462
xx2           2.0773     0.4850   4.283

Best Answer

Related Solutions

ANOVA – What is Compound Symmetry Explained in Plain English?

Paired t-test – Special Case of Linear Mixed-Effect Modeling

Related Question