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

Question

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.

Answer 1

The only way to get statistical properties that work fully "as advertised" is to completely pre-specify a model. When a model is "tweaked" to better fit the data, Bayesian uncertainty intervals of frequentist compatibility (confidence) intervals will be too narrow because they don't select model uncertainty. P-values will be too small. On the other hand, sticking with an ill-fitting model can be a disaster. That's why I say in Regression Modeling Strategies that using the data to select the model is almost as bad as not doing so.

The best solution to my mind is to pre-specify flexible models that contain parameters for things that matter that you don't know. Bayesian priors are the only coherent way to handle much of this. For example you hope a relationship is linear but have a skeptical prior on the amount of nonlinearity as shown here.

Faraway in his paper The Cost of Data Analysis showed that the largest impact of model uncertainty comes from how the response variable is transformed. You are using a very parametric model. You can solve many problems by using Y-transformation-invariant semiparametric models. A detailed case study of using semiparametric ordinal models for continuous Y may be found here.

Comparing a very limited number of models using AIC or its modifications can be OK, but comparing 5 models results in too low a probability that AIC will select the "right" model.