3
$\begingroup$

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

$\endgroup$

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Thank you for your reply. If instead of just $b_{1j}+b_{2j}$ the model would have been non-linear (NLMM) like $f_1(b_{1j})+ f_2(b_{2j})$, where $f_1, f_2$ are some functions, it would have been identifiable, right? Any other way which comes to mind? $\endgroup$ Commented Dec 10 at 11:22
-1
$\begingroup$

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)
  1. 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)
  1. 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)
New contributor
ريماس العنزي is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
3
  • $\begingroup$ Sorry, Rick Hass, I was editing simultaneously that resulted in automatically rejecting your proposed edit. My bad. $\endgroup$ Commented Dec 9 at 15:18
  • 3
    $\begingroup$ This reads like a ChatGPT answer $\endgroup$
    – user20650
    Commented Dec 9 at 15:33
  • 1
    $\begingroup$ It reads like a ChatGPT answer and it results in an error. $\endgroup$ Commented Dec 10 at 11:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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