1
$\begingroup$

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. enter image description hereenter image description here 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!

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ 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. $\endgroup$
    – StoryTeller0815
    Aug 17, 2019 at 12:43
  • 1
    $\begingroup$ 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. $\endgroup$
    – StoryTeller0815
    Aug 17, 2019 at 12:47
  • $\begingroup$ 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? $\endgroup$
    – martina
    Aug 17, 2019 at 13:27
  • 1
    $\begingroup$ 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. $\endgroup$
    – StoryTeller0815
    Aug 17, 2019 at 13:30
  • 1
    $\begingroup$ Does this answer your question? Accounting for heteroskedasticity in lme linear mixed model? $\endgroup$
    – user318514
    Aug 23, 2022 at 19:09

1 Answer 1

Reset to default
2
$\begingroup$

Rather than transforming data to fit a model, I suggest using a model that fits the data.

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks a lot for the suggestion! $\endgroup$
    – martina
    Aug 17, 2019 at 13:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.