How to account for heteroskedasticity in residuals in a fully crossed mixed-effects model (lmer)?

Question

I’m new here so please let me know if I am missing anything in this description/explanation.

I have a 4x4 repeated measures design. My dependent variable is pupil dilation, my two IVs are light level and signal-to-noise ratio (target sentence to background noise). I have 3 random effects: participant, sentence, and trial number. See a snapshot of dataframe below

I have built the model according to the “keeping it maximal” advice in “Random effects structure for confirmatory hypothesis testing: Keep it maximal” by Barr et al (2013). So I have fit the most complex model consistent with the experimental design, removing only terms required to allow a non-singular fit (Barr et al. 2013)

This is my model:

lmer(peakdilation ~ snr * lightlevel + (1 + lightlevel|participant) + (1|trial_exp) + (1|sentence), pup_data, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e9)))

This model has obvious heteroskedasticity in the residuals e.g.

I am really struggling to figure out the best course of action…

• I have looked at log transforming but know this isn’t always the best option as it can make interpretation tricky (also, I had to add 1 to make some numbers positive and after that, using log(peakdilation) in my model did fully not solve the issue).

• I have looked at switching to nlme but have read that they are not particularly suitable for fully-crossed designs

• I have looked at using glmer but am finding it hard to know how to specify the model to account for this heteroskedasticity

• I have looked at using rlmer from robustlmm but again finding it hard to specify and documentation doesn’t seem to explicitly mention heteroskedasticity.

• I have looked at using brms but I am not very familiar with the Bayesian framework at this stage (and currently time-pressured)

Feeling a bit overwhelmed with the amount of information and all the different methods out there. Could anyone here please offer any advice (in accessible language) about how to move forward and what to do?

Thanks very much in advance

Answer 1

Without getting into too much math and simply from my experience - lmer should do the job just fine. A log transform should indeed bring your continuous DV closer to a normal distribution but it doesn't necessarily solve the issue entirely, just brings you closer to BLUE assumptions. Also, from my experience and the literature, mixed models might be relatively robust to violations of BLUE assumptions, as long as they're not to gross. Also, since your residuals exhibit a clear fanning patterns, using some weighted version of the mixed linear model might help. A brief search yielded the WeMix package, which might help. Anyway, when encountering a similar issue with the repeated measures ANOVA in my doctoral research, I had tried different transformations based on the residual skeweness. If neither of them improved or resolved the skeweness, I resorted to other models. Another issue you can explore is using gls (from the nlme package) to explore how different covariance structures assumptions might affect your model (see here). Finally, I know it might be hard, especially when you're time-pressured - but try not to overthink it, and just keep it simple. As long as you can justify your decisions, you should be okay. I hope this helps, good luck!