R: How to fit a linear mixed model with a custom covariance structure for two random intercepts

Question

Suppose I have a dataset with repeated measurements on q clusters. I want to fit an LMM with two random intercepts, on the same cluster, with a non-diagonal covariance structure on the random effects (if it were diagonal I assume the model would be un-identifiable).

That is:

$$y_{ij} = x_{ij}^\top\beta + b_{1j} + b_{2j} + \varepsilon_{ij},$$

where $y_{ij}$ is the $i$-th observation in the $j$-th cluster ($j = 1, ..., q; ~i = 1, ... n_j$). Again if $b_{1j}$ and $b_{2j}$ are independent I assume the model is un-identifiable but I want them to be correlated. Let $b$ be the vector of random effects of length $2q,$ then $b$ distributes Gaussian with mean $0 $ and a $2q \times 2q$ covariance:

q <- 3
rho <- 0.5
sig2b1 <- 1
sig2b2 <- 1
D <- matrix(c(sig2b1 ,rho, rho, sig2b2),nrow = 2) %x% diag(q); D

[1,]  1.0  0.0  0.0  0.5  0.0  0.0
[2,]  0.0  1.0  0.0  0.0  0.5  0.0
[3,]  0.0  0.0  1.0  0.0  0.0  0.5
[4,]  0.5  0.0  0.0  1.0  0.0  0.0
[5,]  0.0  0.5  0.0  0.0  1.0  0.0
[6,]  0.0  0.0  0.5  0.0  0.0  1.0

(Here $q = 3$, on the diagonal $\sigma^2_{b_1} = \sigma^2_{b_2} = 1 $ and we have three variance components to estimate not including the residuals variance $\sigma^2_e$)

Now the model is a standard LMM in vectorized form:

$$y = X\beta + Zb + \varepsilon,$$

Let's get some data:

library(mvtnorm)

b <- t(rmvnorm(n = 1, mean = rep(0, 2 * q), sigma = D))
# can do: plot(b[1:q], b[(q+1):(2*q)]) to see the correlation
z <- rep(letters[1:q], times = c(3,3,2) * 10)
n <- length(z)
Z <- model.matrix(~0 + z)
Z <- cbind(Z, Z) # Notice Z is concatenated to itself, Z is of order n X 2q
x <- rnorm(n)
y <- 1 + 2 * x + Z %*% b + rnorm(n)

Using lme4 with two random intercepts, is surprisingly (sometimes) not bad:

library(lme4)

mod <- lmer(y ~ x + (1|z) + (1|z))

# can do: plot(b[1:q], ranef(mod)$z[,1])
# and: plot(b[(q + 1):(2 * q)], ranef(mod)$z[,2])
# and: VarCorr(mod)

But I want to fit the model with the D covariance structure, where there is some correlation between the random intercepts of each cluster $j.$

Answer 1

Since each $(b_{1j}, b_{2j})$ (independently) follows the same bivariate normal distribution for $j \in \{1, \ldots, q\},$ your model is observationally equivalent (in the sense of an identical marginal response distribution) to $$ y_{ij} = x_{ij}^\top\beta + \gamma_j + \varepsilon_{ij}, $$ where $\gamma_j \overset{\mathrm{i.i.d.}}{\sim} \operatorname{\mathcal N}(0, \sigma_\gamma^2),$ for $j \in \{1, \ldots, q\},$ and $\sigma_\gamma^2 \equiv \sigma_{b_1}^2 + \sigma_{b_2}^2 + 2\cdot\operatorname{Cov}(b_{11}, b_{21}).$
Therefore, your model is non-identifiable.

Answer 2

To fit your linear mixed model (LMM) with the desired covariance structure for the random effects in R, you’ll need to use a package that allows specification of the covariance structure explicitly. Unfortunately, lme4 does not directly support specifying non-diagonal covariance structures between random effects in its lmer function. However, the nlme or brms packages can be used to achieve this. Here’s how you can approach the problem:

1. Using the nlme Package The nlme package allows you to specify a custom covariance structure for random effects through the pdMat classes. For your problem: Steps:
2. Combine the two random intercept terms into a single random term with a custom covariance structure.
3. Use the pdCompSymm class to define a compound symmetry structure (if appropriate), or use pdBlocked for more general cases.

Code Example:

library(nlme)

# Create a grouping factor for each cluster
z <- rep(letters[1:q], times = c(3, 3, 2) * 10)

# Fit the model with a correlated random effects structure
mod_nlme <- lme(
  fixed = y ~ x,
  random = list(z = pdBlocked(list(
    pdSymm(diag(2))  # Symmetric covariance for the two random intercepts
  ))),
  data = data.frame(y = y, x = x, z = z)
)

# Check model summary
summary(mod_nlme)

# Extract random effects
ranef(mod_nlme)

2. Using the brms Package (Bayesian Approach)

The brms package is highly flexible and allows for modeling complex random effects structures, including specifying correlations between them. This package uses Stan under the hood for Bayesian inference.

Steps:

1. Specify the model formula using | ID | gr(ID, cov = TRUE) to include correlated random intercepts.

2. Fit the model using brm.

Code Example:

library(brms)

# Create a grouping factor for each cluster
z <- rep(letters[1:q], times = c(3, 3, 2) * 10)

# Fit the model with correlated random intercepts
mod_brms <- brm(
  y ~ x + (1 + 1 | gr(z, cov = TRUE)),
  data = data.frame(y = y, x = x, z = z),
  sample_prior = "yes"
)

# Check the results
summary(mod_brms)

# Extract random effects
ranef(mod_brms)

# Plot the correlation between the random effects
plot(mod_brms)

3. Simulating Data and Validating Once you fit the model using either of these approaches, validate it by:

• Checking how well the estimated random effects correlate with the true values.

• Verifying the covariance structure using VarCorr in nlme or posterior_summary in brms.

Checking Results:

# nlme
VarCorr(mod_nlme)

# brms
posterior_summary(mod_brms)