# How do I prioritise model diagnostics while considering model selection and parameter uncertainty?

I have fitted a generalized linear mixed model using glmmTMB on the data (110 observations, balanced data) collected from an observational study to understand the effect of each predictor (variable) on response (variable):

fit <- glmmTMB(y ~ x1 + x2 + x3 + x4 + (1 | group), family = ..., data = ...)


Using DHARMa residuals, I scrutinised fit: the DHARMa residuals had uniform distribution and expected dispersion. There was also no pattern when DHARMa residuals were grouped according to group and plotted against fitted values.

When DHARMa residuals were plotted against the predictor x3 and x4, testQuantiles showed evidence for deviations from expected locations with p-values of 0.005 and 0.0008, respectively. The residuals had uniform distribution when plotted against the remaining predictors x1 and x2.

I am now stuck because I have been unable to see potential downfalls of some of the following options I have tried:

1. I have added certain transformation of x4 to see if it improves the model fit (while keeping x4). testQuantiles now returns p-value of 0.204 of x3, but p-value for x4 is now 0.027. We can refer to this model as fit.update.
2. After going through Brooks et al. (2019), I also compared a total of 5 models using bbmle::AICctab that I was interested in, which were also simpler in complexity e.g. fit.simple.1 had formula y ~ x1 + x3 + x4 + (1 | group). The best model (with formula y ~ x1 + x3 + (1 | group) had uniform DHARMa residuals and showed no pattern when grouped according to group. However, it had non-uniform DHARMa residuals against x3.

I am not sure whether I should still consider the best model the best because the non-uniform residuals against x3 gives evidence for model misspecification which might be making the model to have high AICc.

At the same time, there are only a few groups which have high Cook's distance (considering fit). I am concerned that these groups might be influencing model estimates drastically as the observations they contribute to data are clearly extreme (way beyond Q3 + 1.5 * IQR) which I can say using some background knowledge.

Given this information:

1. What steps can I follow to quantify the uncertainty in parameter estimates of fit (or fit.update) given that there is an evidence for model misspecification?

2. What might I miss if I focus on fit.update entirely for statistical inference and ignore the best model obtained using AICc?

Brooks, M. E., Kristensen, K., Darrigo, M. R., Rubim, P., Uriarte, M., Bruna, E., & Bolker, B. M. (2019). Statistical modeling of patterns in annual reproductive rates. Ecology, 100(7), e02706.