# Difficulty of specifying a within-group correlation structure in a linear mixed-effects model

Suppose we have a dataset generated in R - the response variable val is structured with three factors sub, loc and time:

require(MASS); require(nlme)
r1 <- 0.8; r2 <- 0.3; ns <- 20; nt <- 100
S1 <- matrix(c(1,r1,
r1,1),
nrow=2, ncol=2)
S2 <- matrix(c(1,r2,
r2,1),
nrow=2, ncol=2)
tt <- mvrnorm(ns, mu=c(0, 0), Sigma=S1)
dd <- data.frame(sub =rep(paste0('s', 1:ns), each=2*nt),
time=rep(paste0('t', 1:nt), times=2*ns),
loc =rep(rep(c('l', 'r'), each=nt), times=ns),
val =0)
for(s in 1:ns)
dd[((2*s-2)*nt+1):(2*s*nt), 'val'] <- c(mvrnorm(nt, mu=tt[s,], Sigma=S2))


The dataset has the following structure:

> str(dd)
'data.frame':   4000 obs. of  4 variables:
$$sub : Factor w/ 20 levels "s1","s10","s11",..: 1 1 1 1 1 1 1 1 1 1 ...$$ time: Factor w/ 100 levels "t1","t10","t100",..: 1 13 24 35 46 57 68 79 90 2 ...
$$loc : Factor w/ 2 levels "l","r": 1 1 1 1 1 1 1 1 1 1 ...$$ val : num  0.5081 1.0235 0.3486 0.2187 -0.0908 ...


I'd like to implement the following model with the R package nlme and retrieve the parameter r2:

m <- lme(val~0+loc, random=~0+loc|sub, data=dd, correlation = corSymm(form=~0+loc|sub:time))


However, this does not seem to work as I wanted:

Error in lme.formula(val ~ 0 + loc, random = ~0 + loc | sub, data = dd, : incompatible formulas for groups in 'random' and 'correlation'

The reason for the failure seems to be that the within-group variance-covariance matrix is not structured the same as that (or a subset) of the population-level variance-covariance matrix?

Is there a way to implement the model with nlme? Or hack the function lmer in the lme4 package somehow? Thanks!

========= P.S.: an extended case========

The factor loc in the original data had 2 levels. Suppose now there are 3 levels with the factor loc:

require(MASS)
r1 <- 0.7; r2 <- 0.3; r3 <- 0.5
r4 <- 0.6; r5 <- 0.4; r6 <- 0.8
ns <- 20; nt <- 100
S1 <- matrix(c(1,r1,r2,
r1,1,r3,
r2,r3,1),
nrow=3, ncol=3)
S2 <- matrix(c(1,r4,r5,
r4,1,r6,
r5,r6,1),
nrow=3, ncol=3)
tt <- mvrnorm(ns, mu=c(0, 0, 0), Sigma=S1)
dd <- data.frame(sub =rep(paste0('s', 1:ns), each=3*nt),
time=rep(paste0('t', 1:nt), times=3*ns),
loc =rep(rep(c('l', 'm', 'r'), each=nt), times=ns),
val =0)
for(s in 1:ns)
dd[((3*s-3)*nt+1):(3*s*nt), 'val'] <- c(mvrnorm(nt, mu=tt[s,], Sigma=S2))


Can we retrieve the 3 correlations r4, r5, and r6 in this case?

The following approach sort of works, but it's not ideal:

m <- lme4::lmer(val~(0+loc|sub/time), data=dd,
control = lme4::lmerControl(check.nobs.vs.nRE = "ignore"))


Is there a better parameterization that is similar to the last solution by @Sextus Empiricus?

It looks like you can model this as a nested structure

m <- lme(val~0, random=~0+loc|sub/time, data=dd)


You have a random effect for the subjects, for which you compute 20 times a bivariate location sample.

And a random effect for each time nested within the subjects, for which you compute 20 x 100 times a bivariate location sample.

Edit:

you can also use the following formula defining the correlation structure of the locations within the groups time nested inside the subjects. This function corSymm requires that the covariate is integer-valued and in order to do this you can use as.numeric(loc) inside the formula

m <- lme(val ~ 0, random=~0+loc|sub, data=dd,
correlation = corSymm(form=~0+as.numeric(loc)|sub/time))


Expressing the random effect of the subjects as a correlation structure solves the problem with the ambiguity between the random effect and the random error that is described below.

### Modeling with lmer

You can also use lmer, but then you need to suppress the warning using the control parameter.

m <- lme4::lmer(val~(0+loc|sub/time), data=dd,
control = lme4::lmerControl(check.nobs.vs.nRE = "ignore"))


The reason for this is that your model consist of 20 samples from a bivariate distribution, one for each of the 20 subjects

$$\vec{x}_i \sim MVN\left(\mu = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \Sigma = \begin{bmatrix} 1 & \rho_1\\ \rho_1 & 1\end{bmatrix}\right)$$

and 2000 samples from a bivariate distribution, one hundred for each of the 100 times per subject.

$$\vec{y}_{ij} \sim MVN\left(\mu = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \Sigma = \begin{bmatrix} 1 & \rho_2\\ \rho_2 & 1\end{bmatrix}\right)$$

Then your dependent variable is the sum of those where the subscripts $$i$$ and $$j$$ relate to the subject id and time id.

$$\begin{bmatrix}v_1 \\ v_2 \end{bmatrix}_{ij} = \begin{bmatrix}x_1 \\ x_2 \end{bmatrix}_i + \begin{bmatrix}y_1 \\ y_2 \end{bmatrix}_{ij}$$

The linear mixed-effects models use a different model and has an additional residual term.

$$\begin{bmatrix}v_1 \\ v_2\end{bmatrix}_{ij} = \begin{bmatrix}x_1 \\ x_2\end{bmatrix}_i + \begin{bmatrix}y_1 \\ y_2\end{bmatrix}_{ij} + \begin{bmatrix}\epsilon_1 \\ \epsilon_2\end{bmatrix}_{ij}$$

This makes that the correlation structure is dependent from the $$\epsilon$$ and can not be determined uniquely (you can make the same model with different combinations of correlation matrices $$\Sigma$$ and $$\epsilon$$).

The lme function somehow ignores it and estimates an $$\epsilon$$ with a very tiny value for the variance.

### Example Demonstration

below we have a computation where we use lmer to model the part $$\vec{y}_{ij}$$

$$\vec{y}_{ij} \sim MVN\left(\mu = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \Sigma = \begin{bmatrix} 1 & \rho_2\\ \rho_2 & 1\end{bmatrix}\right)$$

as equivalent to

$$\vec{y}_{ij} = \vec{z}_{ij} + \vec{\epsilon}_{ij}$$

with $$\vec{z}_{ij}$$ fully correlated (like a random intercept term) and $$\epsilon$$ completely independent (like a random noise term).

$$\vec{z}_{ij} \sim MVN\left(\mu = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \Sigma = \begin{bmatrix} \rho_2 & \rho_2\\ \rho_2 & \rho_2\end{bmatrix}\right)$$

$$\vec{\epsilon}_{ij} \sim MVN\left(\mu = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \Sigma = \begin{bmatrix} 1-\rho_2 & 0\\ 0 & 1-\rho_2\end{bmatrix}\right)$$

m <- lme4::lmer(val~0+(0+loc|sub)+(1|sub:time), data=dd)
summary(m)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: val ~ 0 + (0 + loc | sub) + (1 | sub:time)
Data: dd

AIC      BIC   logLik deviance df.resid
11632.1  11663.6  -5811.0  11622.1     3995

Scaled residuals:
Min      1Q  Median      3Q     Max
-2.9409 -0.5866 -0.0003  0.6080  3.3442

Random effects:
Groups   Name        Variance Std.Dev. Corr
sub:time (Intercept) 0.3347   0.5785
sub      locl        0.9014   0.9494
locr        0.6643   0.8151   0.84
Residual             0.7488   0.8654
Number of obs: 4000, groups:  sub:time, 2000; sub, 20


In this case the control parameter is not necessary. The modeled covariance of $$y_{ij}$$ is computed from the lines

sub:time (Intercept) 0.3347   0.5785


...

Residual             0.7488   0.8654


$$\hat\Sigma_y = \begin{bmatrix} 0.3347+0.7488 & 0.3347 \\ 0.3347 & 0.3347+0.7488 \end{bmatrix} = \begin{bmatrix} 1.0835 & 0.3347 \\ 0.3347 & 1.0835 \end{bmatrix}$$

• Very interesting! Sorry, I should have added the model formulation in addition to the R code. Actually I first tried lmer(val~1+(0+loc|sub/time), data=dd) with lme4 and it whined with the error message the random-effects parameters and the residual variance are unidentifiable. That's why I didn't even bother this with lme(val~1, random=~0+loc|sub/time, data=dd). I'm baffled now: why does the same model formulation work with nlme, but not lme4? Feb 1, 2022 at 15:21
• You're just amazing! I really appreciate your great help! Feb 1, 2022 at 18:20
• @bluepole I have changed the final part. Using the contrast variable is not neccesary. Feb 1, 2022 at 19:52
• @bluepole I have made an edit that works with the 3 levels. It is in the end not so difficult. (I did not know yet much about these functions like pdSymm and corSymm). Feb 13, 2022 at 23:16
• @bluepole I did not know that either. The pdClass and corClass stuff sounded like magic to me before but now I get them and I could hunt down the problem more easily (which was eventually just starting to read the manual which states that the variable needs to be an integer). Feb 14, 2022 at 14:46