I have run a Gamma model where I try to predict the bioturbation potential by some sediment predictors. We sampled multiple stations over 10 years. Therefore I included station and year as random factors and added an autocorrelation term.
model <- glmmTMB(BPc ~ mudfraction + medsand + X250.500um + shellfraction + (1|StationCode) + (1|Year) + ar1(Year + 0|StationCode),family = Gamma(link = "log"), data = GLMMdataset)
I get following summary output :
Family: Gamma ( log )
Formula: BPc ~ mudfraction + medsand + X250.500um + shellfraction + (1 | StationCode) + (1 | Year) + ar1(Year + 0 | StationCode)
Data: GLMMdataset
AIC BIC logLik deviance df.resid
6398.5 6439.7 -3189.3 6378.5 444
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
StationCode (Intercept) 0.05146 0.2269
Year (Intercept) 0.04565 0.2137
StationCode.1 Year2010 0.16480 0.4060 0.35 (ar1)
Number of obs: 454, groups: StationCode, 93; Year, 9
Dispersion estimate for Gamma family (sigma^2): 0.259
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 12.8665791 0.5721050 22.490 < 2e-16 ***
mudfraction -0.0543623 0.0129689 -4.192 2.77e-05 ***
medsand -0.0102631 0.0008794 -11.670 < 2e-16 ***
X250.500um -0.0464859 0.0047010 -9.889 < 2e-16 ***
shellfraction 0.0092637 0.0027312 3.392 0.000694 ***
And residual diagnostic plots made with the DHARMa package:
I would conclude these plots look OK, but when I look at residual vs predictor I get a deviation for shellfraction at the 0.75 quantile.
How can I address this deviation for shellfraction in the residual plot so my residuals are evenly distributed?