# Linear mixed model with fixed effects in the residual variance

I would like to know how to structure a mixed effect model that allows for a) a fixed factor or b) a covariate to influence residual variance.

I have a dataset of repeated observations from individuals (id) with regards to a behaviour (y). Individuals vary randomly in their mean y value. Additionally, individuals have phenotypic strategies (factoral, Ex1 below) or varying degrees of a personality trait (covariate, Ex2). The phenotype/personality trait does not influence the mean y value of an individual, but does influence the variability an individual shows for trait y.

Consider the following 2 examples:

    library(tidyverse)
set.seed(123)
#Ex1:  fixed effect is a factor

# -- Data simulation
id  <- rep(letters[1:9], each=10)  #10 individuals with 10 observations
strat <- rep(c("a","b", "c"), each = 30) # 3 individuals use each strategy
id.int <- rep(rnorm(9,0,5), each=10)  #individuals vary in intercept
strat.e <- c(rnorm(30,0,1), rnorm(30,0,5), rnorm(30,0,10)) ## strategy b is more variable than strategy a
y      <- id.int + strat.e
data.f  <- data.frame(id,y, strat)

ggplot(data.f, aes(y=y, x = id, fill = strat)) + geom_boxplot()


#Ex2: fixed effect is a covariate
# -- Data simulation
id  <- rep(letters[1:10], each=10)  #10 individuals with 10 observations
strat <- rep(1:10, each = 10) # individuals are linearly ranked on a personality trait
id.int <- rep(rnorm(10,0,5), each=10)  #individuals vary in intercept
strat.e <- sapply(strat, function(x) rnorm(1,0,3*x)) ## variation increases with the personality trait
y      <- id.int + strat.e
data.cov  <- data.frame(id,y, strat)

ggplot(data.cov, aes(y=y, x = id, fill = strat)) + geom_boxplot()


If the residual variance were homoscedastic, this can simply be modeled with a linear mixed effect model, where individual is included as a random effect:

lme4::lmer(y ~ (1|id), data = data.f)


However, for the examples above, I would like to test whether there is additional structure in the residual variance as a result of different strategies, but I am unsure how to specify this structure.

I suspect I need a glmm with structured dispersion or an hglm, but am not confident in my ability to properly specify the matrices. Any suggestions would be much appreciated.

• You might need to switch over to the lme() function from the nlme package for what you are trying to do. See this link for a starting point: quantdev.ssri.psu.edu/sites/qdev/files/…. Mar 13, 2019 at 22:23
• @IsabellaGhement Thank you for the link, having skimmed through the post it looks perfect! I will work through it with my simulated data and hopefully be able to post an answer! Mar 15, 2019 at 8:56
• Sounds good, Morgan! Looking forward to your update - can you tag me in a comment when you post the update? I don't want to miss it! 👍 Mar 15, 2019 at 15:57
• @IsabellaGhement My answer is up. The tutorial you shared was a lifesaver! Mar 21, 2019 at 9:17

I have found two methods of structuring the variance. The first is by structuring the residuals in a linear mixed model, using nlme package. This was based on a tutorial (https://quantdev.ssri.psu.edu/sites/qdev/files/ILD_Ch06_2017_MLMwithHeterogeneousVariance.html) suggested by @IsabellaGhement.

The second is using hglm models with fixed effect in the residual dispersion part of the model.

Both lmm and hglm produce similar results.

library(nlme)
library(tidyverse)
library(hglm)


# Ex 1: Within-individual variation changes with a factor.

Linear model approach

lmm.f = lme(fixed = y ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | strat),  #structure in the residual variance
data = data.f,
method = 'REML')
summary(lmm.f )
VarCorr(lmm.f)[1]  #variance of the random effects (between individual variance)  True value: 5^2

##Within variance
summary(lmm.f)$$sigma ##residual StdDev coef(lmm.f$$modelStruct$$varStruct, unconstrained=FALSE) ##weights #estimated residual var for group = a: (summary(lmm.f)$$sigma*1.0000)^2 ## true value: 1
#estimated residual var for group = b and c:
(summary(lmm.f)$$sigma*coef(lmm.f$$modelStruct$varStruct, uncons=FALSE))^2 ## true values: 25 and 100 ## does structuring the model improve the fit? lmm.f0 = lme(fixed = y ~ 1, random = list(id = pdSymm(form = ~ 1)), ##standard lmm structure data = data.f, method = 'REML') summary(lmm.f0) anova(lmm.f0, lmm.f) ##structured model is better  hglm approach hglm.f <- hglm(fixed = y ~ 1, random = ~ 1|id, disp = ~ strat, data = data.f, calc.like=T) summary(hglm.f) ## mean model intercept = 1.637 (true value 0) ## random effect variance = Dispersion parameter for the random effects = 18.86; true value = 25, lmm = 18.86 ## residual effects = Model estimates for the dispersion term: #a: exp(int): exp(-0.4650) = 0.628; True value: 1; lmm: 0.6281346 #b: exp(int+b): exp(-0.4650+3.4408) = 19.605; True value: 25; lmm: 19.60616 #c: exp(int+c): exp(-0.4650+4.9547) = 89.0; True value: 100; lmm: 89.09715  # Ex 2: Within-individual variation changes with a covariate. Linear model approach lmm.cov = lme(fixed = y ~ 1, random = list(id = pdSymm(form = ~ 1)), weights = varExp(form = ~ strat), data = data.cov, method = 'REML') summary(lmm.cov) VarCorr(lmm.cov) VarCorr(lmm.cov)[1] ## random effect variance: true value: 25 summary(lmm.cov)$$sigma ##Residual Std Deviation coef(lmm.cov$$modelStruct$varStruct, uncons=FALSE) # parameter in the exponential variance function

##for given value strat, we expect the residual within person variance to be..
strat <- as.matrix(1:10)
varpred <- summary(lmm.cov)$$sigma^2 * exp(coef(2*lmm.cov$$modelStruct\$varStruct, uncons=FALSE)*strat)
cbind(strat, varpred)

True.val <- (3*strat)^2
plot(True.val, varpred)

lmm.cov0  = lme(fixed = y ~ 1,
random = list(id = pdSymm(form = ~ 1)), ##standard lmm structure
data = data.cov,
method = 'REML')
anova(lmm.cov0, lmm.cov) #structured model is better


hglm approach

hglm.cov <- hglm(fixed = y ~ 1, random = ~ 1|id, disp = ~ strat, data = data.cov,calc.like=T)
summary(hglm.cov)
## mean model intercept = 2.8718 (true value 0)
## random effect variance = Dispersion parameter for the random effects = 30.23;  true value = 25, lmm =  30.23021
## residual effects = Model estimates for the dispersion term:
#  intercept: 2.8718, slope: 0.4347
varpred.h <- exp(2.8718 + strat*0.4347)
cbind(strat, varpred, varpred.h)

plot(True.val, varpred.h)
$$$$
`