# Residuals vs. Predictor: quantile deviation DHARMa

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?