• Generalized linear mixed models (GLMMs) can be useful for non-normal data with random effects
• Fitting complex GLMMs can be a challenge
• \(\texttt{glmmTMB}\) is a fast, flexible and stable package (Brooks et al. 2017)
• It has many distributions available
• Plus flexible zero-inflated models and hurdle models

• Mountaintop removal mining and valley filling (MTR/VF) is one form of land use which can be a stressor to stream ecosystems.
• In particular salamanders are sensitive to land use changes that impact streams and their catchment
• In this study, the effects of MTR/VF were evaluated on stream salamanders in eastern Kentucky (Price et al. 2016).

Southern two-lined salamander

• Counts of salamanders were recorded at 23 headwater streams (sites); 11 streams located on mined land, and 12 reference streams.

• In each stream, a single 10-m sampling transect was set up and was sampled 4 times (usually monthly) from March through June 2013
• Several variables were recorded prior to sampling, e.g. water temp, the number of cover objects (see Brenee’L et al. (2014) for more details)
• Salamanders were identified to species and in some instances life stage (i.e. adult vs. larva)
• As multiple samples were taken from each site, at different times, we have pseudo-replication. We need a random effect for site to account for this

install.packages("glmmTMB")
library(glmmTMB)
Salamanders

##  site mined       cover sample        DOP       Wtemp
##  VF-1   yes -1.44231718      1 -0.5956834 -1.22937861
##  VF-2   yes  0.29841045      1 -0.5956834  0.08476529
##  VF-3   yes  0.39788060      1 -1.1913668  1.01417627
##   R-1    no -0.44761568      1  0.0000000 -3.02335795
##   R-2    no  0.59682090      1  0.5956834 -0.14434533
##   R-3    no  1.34284703      1  0.5956834 -0.01466007
##   R-4    no -0.94496643      1 -0.5956834 -0.44694425
##   R-5    no  0.49735075      1 -0.5956834 -0.60256655
##   R-6    no  0.14920522      1  0.0000000 -0.09247122
##   R-7    no  1.88993285      1  1.7870502 -0.24809353
##  VF-4   yes -1.59152240      1 -1.1913668  2.00410704
##   R-8    no -0.14920522      1  1.7870502  1.16979857
##   R-9    no -0.04973507      1 -1.1913668 -0.79277159
##  R-10    no  1.19364180      1  1.7870502 -0.65876350
##  R-11    no  1.79046270      1 -1.1913668 -0.10111691
##         DOY spp count
##  -1.4970025  GP     0
##  -1.4970025  GP     0
##  -1.2944669  GP     0
##  -2.7122163  GP     2
##  -0.6868600  GP     2
##  -0.6868600  GP     1
##  -0.9704099  GP     1
##  -0.9704099  GP     2
##  -0.7678742  GP     4
##  -1.5780168  GP     1
##  -1.2944669  GP     0
##  -1.5780168  GP     0
##  -1.4970025  GP     1
##  -1.5780168  GP     0
##  -1.5375097  GP     0

summary(Salamanders)

##       site     mined         cover              sample    
##  R-1    : 28   yes:308   Min.   :-1.59152   Min.   :1.00  
##  R-2    : 28   no :336   1st Qu.:-0.69629   1st Qu.:1.75  
##  R-3    : 28             Median :-0.04974   Median :2.50  
##  R-4    : 28             Mean   : 0.00000   Mean   :2.50  
##  R-5    : 28             3rd Qu.: 0.59682   3rd Qu.:3.25  
##  R-6    : 28             Max.   : 1.88993   Max.   :4.00  
##  (Other):476                                              
##       DOP              Wtemp              DOY         
##  Min.   :-2.1984   Min.   :-3.0234   Min.   :-2.7122  
##  1st Qu.:-0.3018   1st Qu.:-0.6139   1st Qu.:-0.5653  
##  Median :-0.0916   Median : 0.0370   Median :-0.0590  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.0000   3rd Qu.: 0.6032   3rd Qu.: 0.9739  
##  Max.   : 3.1691   Max.   : 2.2094   Max.   : 1.4600  
##                                                       
##     spp         count       
##  GP   :92   Min.   : 0.000  
##  PR   :92   1st Qu.: 0.000  
##  DM   :92   Median : 0.000  
##  EC-A :92   Mean   : 1.323  
##  EC-L :92   3rd Qu.: 2.000  
##  DES-L:92   Max.   :36.000  
##  DF   :92

summary(poisfit)

##  Family: poisson  ( log )
## Formula:          
## count ~ spp + mined + spp:mined + (1 | site)
## Data: Salamanders
## 
##      AIC      BIC   logLik deviance df.resid 
##   1940.2   2007.3   -955.1   1910.2      629 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.3316   0.5759  
## Number of obs: 644, groups:  site, 23
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -3.3771     0.7340  -4.601 4.20e-06 ***
## sppPR              0.9163     0.8367   1.095 0.273435    
## sppDM              2.2513     0.7434   3.028 0.002458 ** 
## sppEC-A            0.6931     0.8660   0.800 0.423494    
## sppEC-L            1.7047     0.7687   2.218 0.026576 *  
## sppDES-L           2.5257     0.7348   3.437 0.000588 ***
## sppDF              2.5257     0.7348   3.437 0.000588 ***
## minedno            4.1109     0.7587   5.418 6.02e-08 ***
## sppPR:minedno     -2.4887     0.8688  -2.864 0.004178 ** 
## sppDM:minedno     -2.1526     0.7554  -2.850 0.004377 ** 
## sppEC-A:minedno   -1.5279     0.8838  -1.729 0.083853 .  
## sppEC-L:minedno   -1.1212     0.7782  -1.441 0.149670    
## sppDES-L:minedno  -1.9527     0.7448  -2.622 0.008748 ** 
## sppDF:minedno     -2.6674     0.7485  -3.563 0.000366 ***
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##  Family: nbinom2  ( log )
## Formula:          
## count ~ spp + mined + spp:mined + (1 | site)
## Data: Salamanders
## 
##      AIC      BIC   logLik deviance df.resid 
##   1663.4   1734.8   -815.7   1631.4      628 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.2842   0.5331  
## Number of obs: 644, groups:  site, 23
## 
## Overdispersion parameter for nbinom2 family ():    1 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -3.3750     0.7576  -4.455 8.40e-06 ***
## sppPR              0.9306     0.8773   1.061 0.288829    
## sppDM              2.2485     0.7878   2.854 0.004314 ** 
## sppEC-A            0.7143     0.9052   0.789 0.430029    
## sppEC-L            1.8130     0.8130   2.230 0.025741 *  
## sppDES-L           2.5111     0.7795   3.221 0.001275 ** 
## sppDF              2.5765     0.7801   3.303 0.000957 ***
## minedno            4.1619     0.7932   5.247 1.55e-07 ***
## sppPR:minedno     -2.5831     0.9328  -2.769 0.005617 ** 
## sppDM:minedno     -2.1495     0.8258  -2.603 0.009245 ** 
## sppEC-A:minedno   -1.5828     0.9461  -1.673 0.094339 .  
## sppEC-L:minedno   -1.3383     0.8493  -1.576 0.115100    
## sppDES-L:minedno  -1.9358     0.8164  -2.371 0.017729 *  
## sppDF:minedno     -2.7426     0.8217  -3.338 0.000844 ***
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note: A larger overdispersion parameter corresponds to a lower variance (which means it may be more appropriate to fit a poisson model)

simulateResiduals(fittedModel = nbfit1, plot = T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.6873376 0.7379601 0.7420719 0.5822611 0.6215828 0.2874355 0.4333506 0.6321975 0.8221086 0.4611435 0.19208 0.2704758 0.5348163 0.2272108 0.2365606 0.5751728 0.005147007 0.5171514 0.00710138 0.4862319 ...

This looks much better.

A zero-inflated model adds structural zeros to allow for the higher number of zeros than expected by a Poisson/negative binomial distribution.

zinbfit1 = glmmTMB(count~ spp + mined + spp:mined + (1|site),
                 zi= ~ mined, 
                 Salamanders, family="nbinom2")

The probability of having a structural zero is always modeled with a logit link to keep it between 0 and 1, i.e.

\[ logit(p) = \beta_0 + \beta_1 minedno \] where \(p\) is the probability of having a structural zero

summary(zinbfit1)

##  Family: nbinom2  ( log )
## Formula:          
## count ~ spp + mined + spp:mined + (1 | site)
## Zero inflation:         ~mined
## Data: Salamanders
## 
##      AIC      BIC   logLik deviance df.resid 
##   1663.1   1743.5   -813.5   1627.1      626 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.1867   0.4321  
## Number of obs: 644, groups:  site, 23
## 
## Overdispersion parameter for nbinom2 family (): 1.22 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -2.7729     0.8290  -3.345 0.000823 ***
## sppPR              0.9175     0.9027   1.016 0.309422    
## sppDM              2.2563     0.8157   2.766 0.005671 ** 
## sppEC-A            0.7027     0.9292   0.756 0.449497    
## sppEC-L            1.7392     0.8416   2.067 0.038768 *  
## sppDES-L           2.4453     0.8064   3.032 0.002426 ** 
## sppDF              2.4955     0.8075   3.091 0.001998 ** 
## minedno            3.5855     0.8470   4.233  2.3e-05 ***
## sppPR:minedno     -2.5425     0.9548  -2.663 0.007747 ** 
## sppDM:minedno     -2.1501     0.8494  -2.531 0.011359 *  
## sppEC-A:minedno   -1.5475     0.9680  -1.599 0.109917    
## sppEC-L:minedno   -1.2429     0.8764  -1.418 0.156136    
## sppDES-L:minedno  -1.8629     0.8393  -2.220 0.026441 *  
## sppDF:minedno     -2.6533     0.8453  -3.139 0.001696 ** 
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.4327     0.6480  -0.668    0.504
## minedno      -2.8269     2.1353  -1.324    0.186

The zero-inflation model estimates the probability of an extra zero, so \(\texttt{minedno} < 0\) means fewer absences in sites with no mining activity;

simulateResiduals(fittedModel = zinbfit1, plot = T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.06354439 0.2394656 0.4106619 0.5396785 0.497431 0.4459096 0.4666671 0.5751449 0.7882012 0.4401021 0.3319441 0.0320186 0.445626 0.07558496 0.06856543 0.03511994 0.1847906 0.3741568 0.888824 0.7399248 ...

Zero-inflation can be used with any of the distributions in glmmTMB, this makes it easy to compare models. For example using AIC

AIC(nbfit1)

## [1] 1663.357

AIC(zinbfit1)

## [1] 1663.074

AIC for zinbfit1 is slightly lower, however, the negative binomial model appears to be just as good at modelling the zeros in the data which is not surprising (Warton 2005).

DHARMa also has a test for zero-inflation, which compares the distribution of expected zeros from simulations against the observed zeros (Hartig 2020)

testZeroInflation(nbfit1)

## 
##  DHARMa zero-inflation test via comparison to
##  expected zeros with simulation under H0 = fitted
##  model
## 
## data:  simulationOutput
## ratioObsSim = 1.0229, p-value = 0.64
## alternative hypothesis: two.sided

The LR test is only approximate for small to moderate sample sizes when testing fixed effects in GLMMs, for a more accurate answer you can use simulation (“parametric bootstrap”).

nbfit3 <- glmmTMB(count ~ 0 + spp + DOP + Wtemp + (1|site),
                 Salamanders, family="nbinom2") #null model
nbfit4 <- update(nbfit3, . ~ . + mined + spp:mined) #alternate model
LRobs <- 2*logLik(nbfit4) - 2*logLik(nbfit3) #observed test stat

nSims <- 1000
LR <- rep(NA,nSims)
for(i in 1:nSims){  
  simCount <- simulate(nbfit3)$sim_1 #simulate data under the null
  tryCatch({
    null <- glmmTMB(simCount ~ 0 + spp + DOP + Wtemp + (1|site),
                 Salamanders, family="nbinom2") 
    alt <- update(null, . ~ . + mined + spp:mined) 
    LR[i] <- 2*logLik(alt) - 2*logLik(null) 
    },
    warning = function(w) {LR[i] = NA}, 
    error=function(e) {LR[i] = NA})  #if model doesn't converge, LR = NA 
}

p <- mean(LR > LRobs,na.rm=TRUE) #P-value

## [1] 0.01

You can see what other covariance structures available here

If we want to take into account that samples taken closer together are more correlated than samples taken further apart (temporal correlations) we could use the autoregressive covariance matrix

Salamanders$sample <- as.factor(Salamanders$sample)

nbfit2 = glmmTMB(count~spp + mined + spp:mined + ar1(sample + 0|site),
                 Salamanders, family="nbinom2")

## Warning in fitTMB(TMBStruc): Model convergence problem;
## singular convergence (7). See vignette('troubleshooting')

Looks like the model is over parameterized!

• Awesome things can be done with \(\texttt{glmmTMB}\)
• It’s fast
• Flexible: more distributions available, e.g. negative binomial (two parameterizations) and Conway-Maxwell-Poisson; zero-inflation can vary across observations and can fit hurdle models
• Can compare competing models with the same package