3
$\begingroup$

I need help fixing the model I landed on through backwards step-wise elimination. I chose a negative binomial model because my variance seems much larger than the mean, with random intercepts from the variable "Site" because my sampling scheme has a repeated measures design (fish counts recorded at the same sites each season/year). I wasn't sure how to make a reproducible dataframe as there are 1000+ rows of data - I attached a subsample of 10 rows to give an idea.

Prior to modelling question: Should SAV (% cover = 0->100) be arcsine transformed for this model?

Main question: I'm not sure what to do about outliers and potential dispersion (p=0.064). After fitting the model, I tried adding "~Site" to the dispformula, but this resulted in model convergence problems. I'm also not sure why my model summaries don't show the random effects section, was this entered correctly?

Zero-inflation doesn't seem to be a problem.

My data (CYR = calendar year, num = count of fish):

> str(toad)
'data.frame':   1262 obs. of  16 variables:
 $ use_for_analysis: chr  "Standard" "Standard" "Standard" "Standard" ...
 $ CYR             : int  2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 ...
 $ Season          : chr  "DRY" "DRY" "DRY" "DRY" ...
 $ Month           : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Site            : Factor w/ 47 levels "1","2","3","4",..: 10 11 12 13 14 15 16 17 9 1 ...
 $ area_sampled    : int  3 3 3 3 3 3 3 3 3 3 ...
 $ common_name     : chr  "Gulf Toadfish" "Gulf Toadfish" "Gulf Toadfish" "Gulf Toadfish" ...
 $ num             : int  0 0 0 0 0 1 0 0 0 0 ...
 $ den             : num  0 0 0 0 0 ...
 $ occur           : int  0 0 0 0 0 1 0 0 0 0 ...
 $ temp            : num  23.3 23.7 22.8 22.6 22.7 ...
 $ sal             : num  30.7 30.9 30.2 30 29.8 ...
 $ DO              : num  6.18 7.56 5.62 5.85 6.49 6.29 5.21 5.36 6.66 5.9 ...
 $ sed_depth       : num  23 48 12 32 43 45 49 32 23 90 ...
 $ water_depth     : num  80 72 98 95 97 85 78 90 58 74 ...
 $ SAV             : num  85.6 53.8 81.6 71.8 91.5 ...

# Variance not equal to the mean
> var(toad$num)
[1] 1.97994
> mean(toad$num)
[1] 0.3716323

toad_example <- toad %>% dplyr::select(CYR, Season, Month, Site, area_sampled, num, temp, 
                               sal, DO, sed_depth, water_depth, SAV)

# Randomly selected rows for display purposes (N = 1262 rows)
toad_example <- toad_example[sample(nrow(toad_example), 10), ] 

toad_example:
     CYR Season Month Site area_sampled num  temp   sal   DO sed_depth water_depth  SAV
1476 2020    WET     9    8            3   0 30.30 20.55 7.83         8          65 93.5
1527 2021    WET     9   40            3   0 28.80 24.28 4.83        37         117 53.5
577  2011    DRY     1   16            3   0 23.43 27.32 9.05        44          78 98.0
668  2012    DRY     1   43            3   0 20.45 25.74 7.17        37          58 87.0
864  2014    DRY     3   36            3   0 22.40 17.31 8.26        61          80 34.0
507  2010    DRY     1   24            3   0 20.41 23.67 8.65        51          20 45.5
1180 2017    DRY     3   45            3   1 23.20 19.87 7.42        95          80 94.0
861  2014    DRY     3   33            3   1 22.40 18.50 9.98        14          75 52.0
506  2010    DRY     1   23            3   0 20.08 23.41 9.55        31          35 49.0
808  2013    WET     9   47            3   1 30.70 24.88 6.45        79          71 79.5

Original model:

full <- glmmTMB(num ~ temp + sal + DO + sed_depth + water_depth + 
               SAV + CYR*Season + (1|Site) + offset(area_sampled), 
               family=nbinom1, toad)

# Warning message: In (function (start, objective, gradient = NULL, hessian = NULL,:
# NA/NaN function evaluation

> summary(full) # AIC 1599.9
 Family: nbinom1  ( log )
Formula:          num ~ temp + sal + DO + sed_depth + water_depth + SAV + CYR *  
    Season + (1 | Site) + offset(area_sampled)
Data: toad

     AIC      BIC   logLik deviance df.resid 
  1599.9   1661.6   -788.0   1575.9     1250 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Site   (Intercept) 0.3403   0.5834  
Number of obs: 1262, groups:  Site, 47

Dispersion parameter for nbinom1 family (): 2.18 

Conditional model:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -3.324e+02  7.549e+01  -4.403 1.07e-05 ***
temp          -3.557e-02  4.003e-02  -0.889   0.3742    
sal            1.027e-02  1.041e-02   0.987   0.3238    
DO             7.336e-02  3.667e-02   2.001   0.0454 *  
sed_depth      6.846e-03  2.976e-03   2.300   0.0214 *  
water_depth    1.725e-03  4.120e-03   0.419   0.6755    
SAV            2.800e-02  4.150e-03   6.748 1.50e-11 ***
CYR            1.615e-01  3.762e-02   4.293 1.76e-05 ***
SeasonWET      4.182e+02  8.526e+01   4.904 9.38e-07 ***
CYR:SeasonWET -2.070e-01  4.224e-02  -4.899 9.62e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Step-wise elimintion:

reduced <- glmmTMB(num ~ sal + DO + sed_depth + water_depth + 
                  SAV + CYR*Season + (1|Site) + offset(area_sampled), 
                  family=nbinom1, toad)

summary(reduced) # AIC 1598.7

reduced2 <- glmmTMB(num ~ DO + sed_depth + water_depth + SAV + 
                   CYR*Season + (1|Site) + offset(area_sampled), 
                   family=nbinom1, toad)

summary(reduced2) # AIC 1597.3

reduced3 <- glmmTMB(num ~ DO + sed_depth + SAV + CYR*Season + (1|Site) + 
                    offset(area_sampled), family=nbinom1, toad)

summary(reduced3) # 1595.8

reduced4 <- glmmTMB(num ~ sed_depth + SAV + CYR*Season + (1|Site) + 
                    offset(area_sampled), family=nbinom1, toad)

summary(reduced4) # AIC 1596.9 (reduced3 is slightly better)

Final model:

final <- glmmTMB(num ~ DO + sed_depth + SAV + CYR*Season + 
                (1|Site) + offset(area_sampled), family=nbinom1, toad)

> summary(final)
 Family: nbinom1  ( log )
Formula:          num ~ DO + sed_depth + SAV + CYR * Season + (1 | Site) + offset(area_sampled)
Data: toad

     AIC      BIC   logLik deviance df.resid 
  1595.8   1642.1   -788.9   1577.8     1253 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Site   (Intercept) 0.3549   0.5958  
Number of obs: 1262, groups:  Site, 47

Dispersion parameter for nbinom1 family ():  2.2 

Conditional model:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -3.054e+02  6.641e+01  -4.599 4.24e-06 ***
DO             6.258e-02  3.528e-02   1.774   0.0761 .  
sed_depth      6.814e-03  2.992e-03   2.277   0.0228 *  
SAV            2.803e-02  4.150e-03   6.754 1.44e-11 ***
CYR            1.479e-01  3.286e-02   4.501 6.77e-06 ***
SeasonWET      3.803e+02  7.735e+01   4.916 8.84e-07 ***
CYR:SeasonWET -1.883e-01  3.837e-02  -4.907 9.23e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residuals:

simulationOutput <- simulateResiduals(final)
plotResiduals(simulationOutput, main=NULL)

enter image description here

testUniformity(simulationOutput)

enter image description here

testOutliers(simulationOutput, type = 'bootstrap')

enter image description here

> testDispersion(simulationOutput, type = "DHARMa")

    DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated

data:  simulationOutput
dispersion = 1.4687, p-value = 0.064
alternative hypothesis: two.sided

Zero-inflation:

> testZeroInflation(simulationOutput, plot = T)

    DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted
    model

data:  simulationOutput
ratioObsSim = 0.99281, p-value = 0.728
alternative hypothesis: two.sided

This model worked better. I had to include Year (categorical) as a random effect to capture the relatedness of observations within years, and use Year (continuous) as a fixed effect to capture any trends in abundance over time.

Additionally, even though sample effort doesn't differ by year, I needed an offset term, because area needs to be logged (negative binomial uses a log-link function and this keeps the units consistent with what is happening to the count).

*Not sure if I need to "center" the variables (grand mean or cluster) or include an optimizer though..

reduced4 <- glmmTMB(num ~ DO + SAV + CYR*Season + (1|Site) + 
                      (1|as.factor(CYR)) + offset(log(area_sampled)), 
                    family=nbinom1, toad)

> summary(reduced4) # AIC 1568.0 (reduced3 is slightly better)
 Family: nbinom1  ( log )
Formula:          num ~ DO + SAV + CYR * Season + (1 | Site) + (1 | as.factor(CYR)) +  
    offset(log(area_sampled))
Data: toad

     AIC      BIC   logLik deviance df.resid 
  1568.0   1614.2   -775.0   1550.0     1253 

Random effects:

Conditional model:
 Groups         Name        Variance Std.Dev.
 Site           (Intercept) 0.5014   0.7081  
 as.factor(CYR) (Intercept) 0.2309   0.4805  
Number of obs: 1262, groups:  Site, 47; as.factor(CYR), 15

Dispersion parameter for nbinom1 family (): 1.88 

Conditional model:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -2.521e+02  9.091e+01  -2.773  0.00555 ** 
DO             7.934e-02  3.800e-02   2.088  0.03683 *  
SAV            3.047e-02  4.515e-03   6.748 1.49e-11 ***
CYR            1.223e-01  4.506e-02   2.714  0.00664 ** 
SeasonWET      3.568e+02  8.365e+01   4.265 2.00e-05 ***
CYR:SeasonWET -1.766e-01  4.150e-02  -4.256 2.08e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

enter image description here enter image description here

Improve this question
$\endgroup$
1

1 Answer 1

Reset to default
4
$\begingroup$

Actually, based on what I see here alone, it doesn't look like your model is in some serious danger. The DHARMa vignette explains that the visual indicators of the residuals are a lot more informative than the p-values, especially when sample sizes are large (there is a fair amount of literature on how normality tests are actually fairly problematic for this very reason alone). Your qq plots and quantile plots don't look seriously bad, as the qq plot is pretty linear and your quantiles dont have any bizarre curves. Your test of dispersion is also showing your model is underdispersed...the same vignette explains that mixed models are biased towards underdispersion and it is usually better to set the alternative argument in the function to "greater" for overdispersion tests.

Here is the vignette if you haven't read it already. I find it to be a good explanation of this package and trouble-shooting common issues. It also has a lot of examples of way more problematic visuals of DHARMa residuals:

https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.