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),
family=Gamma(link = "log"),
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?