How can I test for differences in variation between groups in a mixed model (lme4)?

I would like to test for differences in variation, not in means, between two sites. By looking at a boxplot of my data I see that bird song in one site look much more variable in length than in another site, and also it seems there is more variation within individuals in one site than in the other. In other words, more variation within and between individuals (or so it looks..).

I'm fitting this random model to test for differences:

a   <- lmer (song_length~factor(zone)+date+(1|male/songtype),
data, REML=TRUE)


The reason for the songtype random factor is that birds have different types of songs, and I take measurements of 4 songs per song type per bird.

The question would be: is there more inter- and/or intraindividual variation in site A than in site B?

Also, if there is indeed heteroscedasticity, the model is wrong, and I would need to account for that. But that is a second question I guess.

• (A version of) This analysis can be done in a Bayesian framework with brms which allows to specify a model for all parameters of the response distribution. So for a Normal response you can specify a model for sigma, eg. sigma ~ zone. Once you fit the model you'll have the posterior distribution for the sigma at each site, so you can compare them. Though getting a p-value out of this comparison is not in the Bayesian spirit. Commented Sep 18, 2023 at 23:17

A few points:

• sometimes log-transforming data can clear up heteroscedasticity nicely; this would be my first attempt, especially as you have a positive responses variable (song length), so we would expect it to be skewed if the coefficient of variation is large

• assuming that zone in your model is equivalent to "site" in your verbal description (i.e. "site A vs. site B", is coded in the zone variable in your data), you can use the weights argument in nlme::lme to add heteroscedasticity (within-individual variability) to the model:

  a <- lme(song_length~factor(zone)+date,
random=~1|male/songtype,
data, method="REML",
weights=varIdent(~zone))


You can do this a little bit more awkwardly with lmer by adding an observation-level variable and using a dummy variable to set its value to zero in the first site:

    data <- transform(data,
obs=factor(1:nrow(data)),
zonedummy=as.numeric(zone=="B"))
a <- lmer(song_length~factor(zone)+date+
(1|male/songtype) + (zonedummy-1|obs),
data, REML=TRUE,
control=lmerControl(check.nobs.vs.nRE  = "ignore",
check.nobs.vs.nlev = "ignore"))


You can do a similar thing for the among-bird variation by adding a term of the type (zonedummy-1|male); this adds an additional among-male variance term only in site B (you may have to use additional arguments to lmerControl to override some errors).

More recently, in glmmTMB you can:

a <- glmmTMB(song_length~factor(zone)+date + (1|male/songtype),
data, REML = TRUE,
dispformula = ~zone)


The brms package will also let you fit models like this.

• Hi Ben, thanks for your time. I actually used log-transformed data. But the thing is that I do not consider heteroscedasticity a problem, but the real thibng! I mean, the interesting phenomenon I see here is that there is selection for homegeneity in one site and more variation (perhaps more emigrant birds, who knows), so in this case it is variation what matters. But surely I will use your suggestion of weights to test for means! Commented May 24, 2016 at 21:30

A mixed model may not be appropriate for your research question (about the difference in variance between groups).

Mixed models are used for inference and prediction in data where random effects are assumed to be present, such as in hierarchical data. Typically the interest is in modelling fixed effects while controlling for random effects or in the random effect themselves (such in understanding the variability across levels through variance components).

In most software, such as lme4, it is assumed that there is homoscedasticity. In some software heteroscedasticity can be modeled using weights. But there is no way to test for a difference in variation between groups in the model. Your question is about variation between groups, so one approach is simply to perform an F test of equal variance between site A and site B.

Your model formulation suggests that you are interested in the association between song_length and zone (fixed effect) while controlling for the confounder (date) and the random effects (male and songtype). I don't immediately see how gender can be a random effect - I would expect this to be a fixed effect.

• male is presumably an identifier for a particular male bird ... (note this is not gender, but male ...) Commented May 23, 2016 at 2:07
• Thanks for your answer! Really useful, and yes "male" is the ID of the subject. Best wishes! Commented May 24, 2016 at 21:21
• I get your answer but still find it difficult to understand why one cannot get a step further from controllig for variation in a hierarchical model. The hierarchy here is that I obtained several samples of several song types of several males in two sites, so I have a perfect hierarchical system (samples, songtype, male, site) and it "feels" like these levels of variance that are imbricated within each other must be amenable to testing Commented May 24, 2016 at 21:25