# Mathematical notation in heteroscedastic mixed models

I'm trying to understand how the equations (respecting indices) and the variance and covariance matrices of some mixed models adjusted in R are expressed, which are initially expressed by the lme4 and nlme packages.

The first model is a model of just one random effect on the intercept, this one doesn't seem difficult to me. Adjusted for example as follows in lme4 and nlme respectively:

# Var6 it's a qualitative variable.
# Var7 it's a qualitative variable.
modA_lmer <- lmer(log(Var1)~log(Var2)+log(Var3)+
(Var4)+(Var5)+(1|Var6), data)
modA_lme <- lme(log(Var1)~log(Var2)+log(Var3)+
(Var4)+(Var5),
random = ~1|Var6, data)

modB_lmer <- lmer(log(Var1)~log(Var2)+log(Var3)+
(Var4)+(Var5)+(1|Var7), data)
modB_lme <- lme(log(Var1)~log(Var2)+log(Var3)+
(Var4)+(Var5),
random = ~1|Var7, data)

modC_lmer <- lmer(log(Var1)~log(Var2)+log(Var3)+
(Var4)+(Var5)+(1|Var6)+(1|Var7), data)
modC_lme <- lme(log(Var1)~log(Var2)+log(Var3)+
(Var4)+(Var5),
random= list(Var8= pdIdent(form = ~ 0 + as.factor(Var6)),
Var7 = ~1), data)



I'm not sure (would it be possible for someone to check? Mainly indexes.), but I think that the equation of the three models, respecting the indices, would be given by:

I believe that the variance and covariance matrices D of the random effects would be expressed by

Ok, assuming that everything described above is correct, now imagine the mixed model B in a heterocedastic case, assuming different variances per level of a dummy variable Var 8, which adjusted through the lme4 and nlme packages would look like:

# Var 8 it's a dummy variable
modD_lmer <- lmer(log(Var1)~log(Var2)+log(Var3)+
(Var4)+(Var5)+(Var8-1||Var6), data)
modD_lme <- lme(log(Var1)~log(Var2)+log(Var3)+
(Var4)+(Var5),
random = list(Var6=pdDiag(~Var8)), data)


In this case, what would the equation of my model look like? I thought of the following situation:

Is this correct? The equation would be exactly the same but the difference would be the introduction of an index in the sigma^2 of the error?

The other question is, well, in this case, how would the matrix of variances and covariances look like? Would I have three sigmas in the diagonal matrix? I.e:

I found some information in :

• Hello! I was supposed to have made this post on my account here (without realizing it, I ended up doing it in a way to not be registered). How can I resolve? I already sent a message to support to transfer this post to my main account here. Commented Aug 1, 2023 at 16:35
• @user55546, You can ask to merge accounts. Don't know if both accounts need to be registered.
– Ute
Commented Aug 2, 2023 at 10:54
• I did it, it worked!!! Commented Aug 2, 2023 at 15:39
• Wow, you threw all your rep from the merge into a bounty - it's a very useful question :-) Hope you get some kind of reference answer that is useful for future cases, too
– Ute
Commented Aug 3, 2023 at 15:31

There's a lot going on here.

### Models A,B,C

Models A,B,C are fine in both R-Code and formula, however i have some notes.

1. Var6 and Var7 are crossed random effects, i.e. their groupings are different/ not nested (also see here: Crossed vs nested random effects: how do they differ and how are they specified correctly in lme4?), which means that it is not sensible to speak about their covariance matrix, because they don't "realize" on the same objects.
2. lme does not inherently support crossed random effects, but you seem to have used the approach from here: https://stats.stackexchange.com/a/552180/341520 for modC_lme, which relies on the Var8 being a second intercept column, i.e. all 1s.
3. If $$k$$ is the index of the individual observation, then $$\epsilon_{ijk}$$ is redundant, but not wrong. You might be interested in the equatiomatic package which can create a Latex-formula from most lmer-objects. See the demonstration below.
library(lme4)
library(nlme)
library(equatiomatic)
n <- 700
# only two fixed effects to reduce redundancy
x1 <- rexp(n)
x2 <- runif(n)
var6 <- sample(as.factor(1:40), replace = T, size = n)
ran_eff6 <- rnorm(40, sd = 1.5)
var7 <- sample(as.factor(1:50), replace = T, size = n)
ran_eff7 <- rnorm(50, sd = 2.5)
# dummy for model D, not the right on for modC_lme!
var8 <- sample(as.factor(1:2), replace = T, size = n)
# model C being true
yC <- exp(3*log(x1) + 2*x2 + ran_eff6[as.integer(var6)] + ran_eff7[as.integer(var7)] + rnorm(n, sd = 4))
dat <- data.frame(x1, x2, var6, var7, var8, yC)

modC_lmer <- lmer(log(yC)~log(x1)+(x2)+(1|var6)+(1|var7), data = dat)

extract_eq(modC_lmer)


\begin{aligned} \operatorname{yC}_{i} &\sim N \left(\alpha_{j[i],k[i]} + \beta_{1}(\operatorname{\log(x1)}) + \beta_{2}(\operatorname{x2}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for var7 j = 1,} \dots \text{,J} \\ \alpha_{k} &\sim N \left(\mu_{\alpha_{k}}, \sigma^2_{\alpha_{k}} \right) \text{, for var6 k = 1,} \dots \text{,K} \end{aligned}

### Model D

Here you seem to have gotten confused reading https://rpubs.com/bbolker/factorvar, because your formula describes the simple case of varying variance in $$\epsilon_{ijk}$$, while your R-Code tries to solve the difficult case of varying variance in $$u_i$$. Both are correct, they just don't fit together.

Adjusting the formula to the hard case I'd propose writing $$u_{i,l}\sim \mathcal N(0, \sigma_{1, l}),\epsilon_{i,l,k}\sim \mathcal N(0, \sigma_3)$$, while if your formula is correct we can use weights argument in lme just as Ben Bolker did in your source:

# generate y with simple heteroscedacity
dat\$yD <- exp(3*log(x1) + 2*x2 + ran_eff6[as.integer(var6)] + rnorm(n, sd = ifelse(var8 == "1", 3, 6)))
modD_lme <- lme(log(yD)~log(x1)+(x2), random = ~1|var6, weights = varIdent(form=~1|var8), data = dat)
summary(modD_lme)

Linear mixed-effects model fit by REML
Data: dat
AIC      BIC    logLik
4021.07 4048.351 -2004.535

Random effects:
Formula: ~1 | var6
(Intercept) Residual
StdDev:    1.287001 5.716438

Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | var8
Parameter estimates:
2         1
1.0000000 0.5062947
...more stuff...


In rnorm(n, sd = ifelse(var8 == "1", 3, 6)) we defined $$\sigma_{3, 1} = 3$$ and $$\sigma_{3, 2} = 6$$, which correspond to the residual variance of $$\hat\sigma_{3, 2} = 5.716438$$ and the parameters under Variance function show that this for group $$l = 2$$ while group $$l = 1$$ has roughly half that standard deviation. So it fits great! lme4 has no comparable functionality see here: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#spatial-and-temporal-correlation-models-heteroscedasticity-r-side-models Ben Bolkers GLMM-FAQ is frankly always worth reading, when dealing with mixed models.

A covariance matrix for the random effects is still not sensible.

• Hi Lukas! I didn't understand some topics. In the case of model D, how many variance components would I have? Do you mean that my model D is not possible to express mathematically by writing the matrix of variance and covariance? If possible, what would it look like? The matrix? Would it or wouldn't it be the matrix D that I wrote? Commented Aug 4, 2023 at 20:28