# How can I adress problems of heteroscedasticity in mixed model analysis?

I am analizing pupil size data using mixed model analysis in R. I use lme() from package nlme. However, I am encountering serious problems of heteroscedasticity and violation of normality assumption. Do you have any suggestions on how to solve this problem?

I would like to test whether active neurostimulation (vs. control; Stimulation) increases pupil size over time (before stimulation vs. end stimulation; Time). We employed a between subject design.

I adopted the following strategies that did not solve the problem: 1. transformed the outcome - log(). 2. modelled the covariance-variance matrix following the instructions in this page https://rpsychologist.com/r-guide-longitudinal-lme-lmer#heteroscedasticity-at-level-1 .

Here the models: Here the models: Time (0 vs . 1) Stimulation (0 vs 1) Tonic_PS = pupil size ID = participant ID

model.1 <- lme(Tonic_PS ~ Time * Stimulation,,
random = ~ Time|ID,
weights = varIdent(form= ~ 1 | Stimulation),
data =Tonic_pupilsize_T2,
method = "ML")

model.2 <- lme(Tonic_PS ~ Time * Stimulation,
random = ~ Time|ID,
weights = varIdent(form= ~ 1 | Stimulation * Time),
data =Tonic_pupilsize_T2,
method = "ML")

model.3 <- lme(Tonic_PS ~ Time * Stimulation,
random = ~ Time|ID,
correlation = corAR1(),
data =Tonic_pupilsize_T2,
method = "ML")

model.4 <- lme(Tonic_PS ~ Time * Stimulation, random = ~ Time|ID,
weights= varIdent(form= ~ 1 | Time),
correlation = corAR1(),
data =Tonic_pupilsize_T2,
method = "ML")


Any suggestion will be very appreciated!

• Could you please be more specific regarding why you think that your second approach did not solve the problem? Because it actually should - except if group is not the cause of the heteroscedasticity. Aug 17, 2019 at 12:43
• Add: If the basis for your judgment is the second plot: Modeling heteroscedasticity does not vanish if you model it. You just tell the model that level-1-residual variance varies as a function of Stimulation and Time - nothing more and nothing less. If this is the correct model for the variance part than you simply get consistent standard errors, nothing else should change. Aug 17, 2019 at 12:47
• Thanks a lot for your response. I indeed judged based on the graph. I compared the above described models with the simpler model that includes only Random slope and Intercept. model <- lme(Tonic_PS ~ Time * Stimulation, random = ~ Time|ID, data =Tonic_pupilsize_T2, method = "ML") The model fit of the above described models is better looking at BIC AIC. Is this sufficient? Could you be more specific about how I can establish looking at the standard errors if I modelled correctly the variance part? Aug 17, 2019 at 13:27
• You actually can‘t really - just like you don’t know if you have all relevant predictors for the mean you don’t know that for the variance either. For your case, estimating separate variances for each cell seems fine to me - but you could of course check for a robust method as suggested in the other answer. Aug 17, 2019 at 13:30
• Does this answer your question? Accounting for heteroskedasticity in lme linear mixed model?
– user318514
Aug 23, 2022 at 19:09

Quantile regression does not assume homoscedasticity and quantile regression for multilevel (or mixed) models now exists. There is an R package lqmm available from CRAN that should work.