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

Question

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.

Answer 1

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.