0
$\begingroup$

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

Example1

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

$\endgroup$
  • $\begingroup$ 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/…. $\endgroup$ – Isabella Ghement Mar 13 at 22:23
  • 1
    $\begingroup$ @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! $\endgroup$ – Morgan.Brown Mar 15 at 8:56
  • $\begingroup$ 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! 👍 $\endgroup$ – Isabella Ghement Mar 15 at 15:57
  • 1
    $\begingroup$ @IsabellaGhement My answer is up. The tutorial you shared was a lifesaver! $\endgroup$ – Morgan.Brown Mar 21 at 9:17
0
$\begingroup$

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) 
```
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.