# Dealing with heavy-tailed residuals when fitting hierarchical linear models using lme4

This is my first time posting, so please excuse any issues with respect to my description of the problem and the presentation of the data and code I have supplied.

Summary of the Design

30 listeners evaluated 2 speech samples from 19 speakers. When listening to each sample, they transcribed it, and rated it for comprehensibility and accentedness using 100-point sliding scales. Based on the proportion of words correctly transcribed, the sample received an intelligibility score out of 1.

Summary of the Problem

I'm having some trouble fitting a mixed-effects model to this data set since the majority of observations of the dependent variable (~65%) occur at the maximum value (in my case, 1, with values ranging from 0-1). The model residuals have a very heavy lower tail (see image below). I have tried transforming the data to improve the normality of the residuals, first using a log transformation and then by attempting other transformations with the bestNormalize package, but none of the transformations seem to help. I also tried to implement robust mixed-effects models using the robustlmer package. Adjusting the tuning parameters following Koller (2016) did not seem to help much, but admittedly, I'm not very familiar with the package and so may not have implemented all of the functionality properly and/or optimally.

I have two primary questions.

• What is the interpretation of the model given that the residuals are not normally distributed? Specifically, does the heavy lower tail mean the model is less reliable at lower values of the dependent variable? I would appreciate any accessible references that talk about the interpretation of mixed-effects models when assumptions are violated.
• What can I do to resolve this issue? I know that I could convert the continuous variable to a binary one and fit a generalized model, which would not impose the assumption that residuals are normally distributed, but is there any way to work with the continuous data?

Sample Data and Analysis

In creating this example, I have referenced Violated Normality of Residuals Assumption in Linear Mixed Model

Variables

• participant: Categorical/Grouping, 1-19
• listener: Categorical/Grouping, 1-30
• comp: Continuous, 1-100
• accent: Continuous, 1-100
• intelligibility: Continuous, 0-1

I'm interested in predicting intelligibility as a function of comp and accent (both scaled and centered), with by-participant and by-listener random effects. I also have a number of covariates that I have not included here for the sake of streamlining.

#Load required packages
library(lme4)
library(lattice)

#Create sample data set (first 100 observations of 1140)

listener <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3)

participant <- c(500, 500, 510, 504, 502, 518, 514, 512, 521, 512,
522, 515, 502, 509, 509, 507, 501, 501, 503, 510, 508, 514, 513,
521, 516, 516, 522, 504, 523, 523, 508, 511, 507, 518, 513, 515,
503, 511, 508, 509, 504, 504, 508, 522, 502, 503, 511, 523, 521,
510, 514, 500, 523, 501, 507, 518, 514, 510, 511, 507, 500, 512,
521, 515, 515, 516, 516, 512, 518, 509, 513, 501, 522, 502, 503,
513, 509, 503, 523, 503, 507, 502, 514, 521, 510, 513, 500, 504,
516, 516, 502, 518, 501, 509, 522, 518, 512, 501, 512, 513)

accent <- c(68, 32, 12, 7, 13, 25, 5, 8, 25, 3, 55, 41, 12, 39, 80, 20, 14,
46, 3, 6, 45, 31, 73, 58, 31, 34, 28, 21, 48, 43, 31, 46, 6, 5, 45,
88, 5, 33, 10, 10, 10, 9, 10, 9, 9, 9, 10, 10, 11, 9, 9, 9, 9, 9, 9,
9, 9, 9, 8, 9, 11, 9, 12, 9, 9, 9, 10, 9, 9, 12, 9, 9, 12, 12, 9, 9,
86, 13, 27, 1, 5, 27, 11, 18, 1, 63, 51, 2, 55, 42, 29, 59, 3, 88, 1,
15, 0, 38, 1, 55)

comp <- c(68, 71, 22, 22, 44, 18, 25, 4, 80, 2, 57, 28, 91, 77, 75, 21, 64, 57,
2, 10, 57, 72, 88, 80, 47, 53, 53, 56, 75, 62, 77, 28, 7, 6, 9, 39, 34,
18, 50, 55, 45, 50, 55, 47, 55, 41, 44, 49, 53, 45, 42, 40, 48, 50, 44,
44, 50, 53, 43, 47, 55, 40, 55, 45, 42, 48, 50, 45, 40, 57, 47, 51, 55,
55, 44, 55, 100, 53, 72, 52, 10, 100, 43, 65, 53, 97, 100, 57, 100, 84,
69, 89, 70, 100, 18, 34, 23, 97, 8, 98)

intelligibility <- c(0.692307692, 1, 0.727272727, 0.8, 1, 0.8, 0.909090909,
0.714285714, 1, 0.714285714, 0.866666667, 0.75, 1, 1, 1,
.666666667, 1, 1, 0.777777778, 0.666666667, 1, 1, 1, 1, 1,
1, 0.909090909, 0.8, 1, 0.5, 0.916666667, 0.8, 0.428571429,
0.909090909, 0.666666667, 0.785714286, 0.75, 0.8, 1, 1, 1,
1, 0.833333333, 1, 1, 0.875, 1, 1, 1, 0.636363636,
.818181818, 0.846153846, 1, 1, 0.666666667, 0.909090909, 1,
0.833333333, 1, 1, 1, 0.571428571, 1, 0.928571429, 0.9375,
0.882352941, 1, 1, 0.933333333, 1, 0.916666667, 1, 1, 1, 1,
1, 1, 1, 1, 1, 0.833333333, 1, 1, 1, 0.909090909, 0.75, 1,
1, 1, 1, 1, 0.909090909, 1, 1, 0.933333333, 0.6, 1, 1,
0.857142857, 1)

data <- data.frame(as.factor(listener), as.factor(participant), accent, comp, intelligibility)

#Run model, centering accent and comp
fm <- lmer(intelligibility ~ scale(comp, center = T) + scale(accent, center = T) +
(1 | participant) +
(1 | listener), data = data)

#Generate QQ plot to examine distribution of model residuals
qqmath(fm)



Thank you in advance for any help and/or references you can provide.

Update

I tried implementing the zoib package. Here is the model I tried to fit based on the examples provided in https://journal.r-project.org/archive/2015/RJ-2015-019/RJ-2015-019.pdf:

zoib(intelligibility ~ comp + accent | 1 , data = data, zero.inflation = F, one.inflation = T, joint = F, random = 1, EUID = data$participant, n.iter = 5000, n.thin = 50, n.burn = 200) First of all, it's returning an error: "object 'zname' not found." I also have a couple of questions about specifying multiple random effects, as well as the syntax that the model takes. 1. How do I specify two random effect groupings, such as data$participant and data\$listener? All of the examples seem to only show one.
2. I'm unclear what to put to the right of the pipe. In the examples, they show multiple pipes, such as yield ~ temp | 1 | 1. I'm used to putting the random grouping to the right of the pipe in lme4, such as 1 | participant, but that does not seem to be the case for the zoib package.

Again, thanks for your help! I'm very interested in these models since they seem to be a good solution for the type of data I'm dealing with.

• This is a great first post! Given the large number of 1s in your outcome, have you considered using a zero/one inflated beta regression model? See the zoib() package, which uses Bayesian estimation to fit the models. Apparently brms(), another Bayesian package can be used to fit these as well. See journal.r-project.org/archive/2015/RJ-2015-019/RJ-2015-019.pdf and vuorre.netlify.com/post/2019/02/18/… – Erik Ruzek Jan 2 at 23:40
• I would check GLMMadaptive, it's hurdle.beta.fam() looks highly relevant. @DimitrisRizopoulos – usεr11852 Jan 3 at 0:09
• @ErikRuzek, I added an update to my original post about trying to fit a one-inflated model using the zoib package. – Charlie Nagle Jan 3 at 5:31
• @usεr11852saysReinstateMonic This looks promising, but unfortunately nearly all of my models contain more than one random effect grouping factor. – Charlie Nagle Jan 3 at 5:33
• Perhaps not a great option, but with these packages, as opposed to rolling your own Stan or JAGS model from scratch, you could treat one of your factors (listener or participant) as fixed and include dummy variables in the model for the fixed factor. The zoib() syntax requires at least 4 pipes, each representing a different aspect of the model - covariates for the function of the mean of the beta distribution, the shape parameters, the link function Pr(y=0), and the link function Pr(y=1|y>0). This is all in the documentation. – Erik Ruzek Jan 3 at 21:30

The Bayesian options mentioned in the comment thread appropriately treat the outcome as if it were drawn from a beta distribution, but because the beta distribution restricts values to be between >0 & <1, zoib and brms allow you to model the inflation you have at 1. They do so by estimating a separate model for the range >0 & <1, and then for the 1 vs. all else. This is useful if you are substantively interested in the 1 vs. all else comparison. If you are not interested in it, then you can still use a model that treats your outcome as being drawn from a beta distribution, but you need to apply a transformation to nudge the 1 values to something slightly less than 1. See Beta regression of proportion data including 1 and 0.

Below is code that I used on your data to transform and run the model. Note that I re-scaled your predictors manually and that with listener only having 3 categories in the sample data, glmmTMB had problems with a singular fit.

First transform intelligibility using Smithson & Verkulien (2006) formulation 1 - (y * (n−1) + 0.5) / n

data$$intel_t <- ((data$$intelligibility * (length(data)-1)) + 0.5) / length(data)


Next, I used glmmTMB() with family=beta_familiy()

require(glmmTMB)
fm_b1 <- glmmTMB(intel_t ~ comp_c + accent_c + listener + (1 | participant) ,
data = data, family=beta_family())
summary(fm_b1)


This resulted in the following model results, note the "fixed effects" for listener:

Family: beta  ( logit )
Formula:          intel_t ~ comp_c + accent_c + listener + (1 | participant)
Data: data

AIC      BIC   logLik deviance df.resid
-237.7   -219.5    125.9   -251.7       93

Random effects:

Conditional model:
Groups      Name        Variance Std.Dev.
participant (Intercept) 0.005701 0.0755
Number of obs: 100, groups:  participant, 19

Overdispersion parameter for beta family (): 18.9

Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.74746    0.11807  14.801  < 2e-16 ***
comp_c       0.52606    0.09512   5.530 3.19e-08 ***
accent_c    -0.21195    0.09145  -2.318   0.0205 *
listener2    0.18908    0.17815   1.061   0.2886
listener3    0.13371    0.18561   0.720   0.4713
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Regression coefficients from models with a beta family link are somewhat tricky to interpret. To aid interpretation I used ggeffects for comp_c:

require(ggeffects)
ggpredict(fm_b1, "comp_c") %>% plot()


And then also for accent_c:

ggpredict(fm_b1, "accent_c") %>% plot()


The question of residuals is a bit more complicated. I used DHARMa() to look at residuals, but it is not fully integrated with glmmTMB. Of note is that the qqplot does not display the same problematic pattern you saw in your lmer model:

library(DHARMa)
res = simulateResiduals(fm_b1)
plot(res, rank=T)


As a last note, you can also estimate this same model with GLMMadaptive and you could play with the beta.hurdle.fam() as suggested in the comments.

• This is great, thank you for your help! This will make it much easier for me to use this type of model, and the code and explanation you have provided is very clear. – Charlie Nagle Jan 5 at 18:41
• Glad it was helpful! As mentioned, you should consider the zero-one inflated version of the beta model at the very least to feel strong about any conclusions you get from the analysis in which you nudge the intelligibility values away from 1. – Erik Ruzek Jan 6 at 21:42