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)+
                           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)+
                           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)+
                           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:

enter image description here

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

enter image description here

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)+
                           random = list(Var6=pdDiag(~Var8)), data)

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

enter image description here

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:

enter image description here

I found some information in :

  • $\begingroup$ 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. $\endgroup$
    – user55546
    Commented Aug 1, 2023 at 16:35
  • 1
    $\begingroup$ @user55546, You can ask to merge accounts. Don't know if both accounts need to be registered. $\endgroup$
    – Ute
    Commented Aug 2, 2023 at 10:54
  • $\begingroup$ I did it, it worked!!! $\endgroup$
    – user55546
    Commented Aug 2, 2023 at 15:39
  • $\begingroup$ 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 $\endgroup$
    – Ute
    Commented Aug 3, 2023 at 15:31

1 Answer 1


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


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

  • $\begingroup$ 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? $\endgroup$
    – user55546
    Commented Aug 4, 2023 at 20:28

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.