# Model residual diagnostics of gamma GLMM with log-link

I am trying to model fish length data (N > 115.000) which are highly right-skewed using linear mixed effects models. Actually all data make sense and the the extreme high values are valid measurements. My data look like: .

I thought of using a GLMM of the family "Gamma" (with log-link) as length is a censored response (0,∞). Specifically, I am using the function glmmTMB() to fit following model structure:

    mod <- glmmTMB(L ~ gear + (1|species) + (1|location),
data=fish_df)


The model converged without any errors/warnings and the parameter estimates are reasonable, I think. However, I wanted to check/diagnose the model, especially the residuals. I was following some posts in other threads such as "What are the assumptions of a Gamma GLM or GLMM for hypothesis testing?" and "Coping with underdispersion after DHARMa diagnosis on GLMM with Gamma link=log", as I am observing similar issues of rather extreme values in my response that lead to skewed residuals of my models, I guess. So it seems that there is some heteroscedasticity in my residuals. I mean, given the large sample size it is not too surprising to see also supposed 'outliers'.

    ypred = predict(mod)
res = residuals(mod, type = 'pearson')
plot(ypred,res)
hist(res)


I further checked with the DHARMa package for some patterns which are difficult to interpret to me.

    sim_res <- DHARMa::simulateResiduals(modX5, 250)
plot(sim_res,asFactor = F)
DHARMa::testDispersion(mod)
DHARMa::testDispersion(mod, type="PearsonChisq",
alternative="greater")


The dispersion tests do not detect any significant overdispersion in the residuals but rather low dispersion parameters (for PearsonChisq = 0.13703, for DHARMa nonparametric dispersion test = 0.28605). In fact the gamma has already a dispersion parameter, so there should be no 'overdispersed' gamma according to this post.

So I am wondering whether the I could trust the model outcomes and estimates for fixed effects? To have more robust/reliable information on parameter estimates I would use parametric bootstrapped CI e.g. of the fixed effects via bootMer (but this does not solve the model fit per se). So are there any other recommendation I could try to improve the model fit, or is the fit sufficiently okay?