DHARMA Diagnostics shows KS deviation in glmm with Gaussian distribution

[BaptisteSadoughi]

Hello Florian,

It's the first time I use DHARMA as a diagnostic tool and I went through the tutorial and some previous issues (eg. #181) but I am still not sure what would be the best strategy given the results.

I am fitting a mixed model with random slopes in lme4. The response is a composite score expressing an animal individual social tendency with no theoretical upper bound (but in practice in the dataset between 0 and 1.44) and a sample size = 540. Histogram of the response . The residuals looked really wrong, so I applied a square root transformation to the response. Histogram of the response after square-root transformation

I investigate the relationship between the response and individual's age, each observation belonging to a given individual, with repeated measures in animal social groups, over time (season_year). The whole model with other fixed control factors and maximal random slope structure (all predictors dummy coded and centered in the random slope part of the model) is lmer(sq.response ~ z.age + reproductive status season + z.indiviudal_rank + (1 + z.age + repro.1.cent + season.nm.cent + z.indiviudal_rank | actor) + (1 + z.age + repro.1.cent season.nm.cent + z.indiviudal_rank || group) + (1 + z.age + repro.1.cent + z.indiviudal_rank || season_year), data = t.data.out, REML = FALSE, control = lmerControl(optimizer="bobyqa", optCtrl = list(maxfun=100000))).

The first diagnostic plots shows some heteroscedasticity.

Is the deviation large enough for concerns? How would one interpret that the qqplot suggests deviation in the centre of the distribution, whereas the residuals vs. predicted suggest low values are problematic?

Plotting residuals against each fixed effects suggested similar deviation each time but no super clear pattern to me. Vs. age

vs. individual rank

vs. reproductive state

vs. season

I tried adding squared age and squared individual rank to the model, but it did not improve the outcome. Would you kindly have some advice?

Thank you in advance for your time, Baptiste

[florianhartig]

Hello Baptiste,

What your plots tell you is that there is a misfit that is not random.

What DHARMa can't tell you is how much impact that has on the results you are interested in, because this is highly data, model and also question-specific.

I take it that you are mainly interested in the slope for age. Based on what I see, I would assume that results are more or less robust for that variable, just based on the observation that the misfit in slope is minor. Of course, you could introduce a quadratic term for age to improve fit, but this will be at the cost of interpretability.

My only real concern would be the misfit in the general plot for small predictions. I'm not sure where this is coming from. Misfit at the borders always has some potential to affect slopes, so it may be useful to investigate a bit when the model is making small predictions. As we don't see much in the plots against the predictors, it may also be the sign of a non-linear interaction.

Cheers, F

[BaptisteSadoughi]

Hello Florian,

many thanks for your input! Adding square terms for both continuous predictors tends to improve a bit the misfit for small predictions. However, the normality of residuals is harmed by the inclusion of these terms and the plot vs predictors are still not great.

Model including elo, elo^2, , age and age^2

Residuals vs elo in a model including both elo and elo^2 (and age and age^2)

Residuals vs age in a model including both age and age^2 (and elo and elo^2)

The correction on misfit for small values is weaker when only elo and elo^2 are included (in a model with only simple term age) but it improves normality of residuals and residuals vs predictors

Residuals vs elo in a model including both elo and elo^2 (and age)

Residuals vs age (in a model including elo and elo^2)

The model using age and age^2 but only elo gives similar outcome, with strong corrections of misfit for small values but pretty bad heteroscedasticity and wrong residuals vs predictors for both covariates and factors in the model...

I was surprised that the quadratic term for age did not improve things so much, and given the results I would be tempted to use the last version (main term for age + elo + elo^2). Could you maybe explain what you meant by "it may be useful to investigate a bit when the model is making small predictions"? Do you expect all those datapoints may share specific attributes?

Thanks again for your time and looking forward reading from you, Cheers, Baptiste

[florianhartig]

Hello Batiste,

regarding your last question: exactly, I just mean because in the overview plot, you see that the misfit occurred in particular for small predictions, I thought I would investigate under which parameter combinations those occur. But this problem has been resolved.

It is true that the plots don't look great, but I also wouldn't say that there is a major misfit. Two thoughts:

a) regarding misfit, you can always play around with alternative transformations of x, y, interactions, or fit splines (via mgcv::gam), so in principle it is always possible to get a perfect fit, but then, the question is if the final model is still interpretable.

b) One peculiarity of DHARMa is that per default, residuals are calculated unconditional, i.e. the residual is calculated as a aggregated residual of residual error and RE error, and implicitly, your plots check fixed and random effect structure together. If there is a misfit in the RE structure (e.g. because REs are not normal) this would also show up in the residual plots. If you want to condition on the fitted REs, you can do this by supplying re.form = NULL in simulateResiduals function. In this way, residual correspond to the residual function of lme4, which is also very useful btw. More comments on residual checks for LMMs here https://theoreticalecology.github.io/AdvancedRegressionModels/2C-RandomEffects.html#observation-level-residual-plots

Cheers, F

[florianhartig]

p.s.: just to re-iterate: DHARMa highlights that deviations are significant, which is something that is not done in standard plot functions, but see my comments in the vignette: significant just tells you the pattern is not random. It doesn't say it's a big problem. I don't think that these residuals look much worse than in many published papers. Maybe check with the standard plot functions to see how this would translate to the standard normal residual plots that you may be more used to.

[BaptisteSadoughi]

Hello Florian,

wow thank you for the explanation re.form = NULL totally improved the residuals and everything looks "fine" ! I will read the specific section on the link to be sure I fully understand what happens when excluding the RE structure.

Thanks again for your time this was super helpful ! Best

[florianhartig]

That's great. Note that a possible explanation for this is that the REs deviate from normality or correlate with predictors, which you could check as well.