Linear mixed model with fixed effects in the residual variance

Question

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)
head(data.f)

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)
head(data.cov)

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.

Answer 1

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