This is the code for running the GLMM:s from the article The potential of semiochemical blends to monitor saproxylic beetle communities published in Journal of Insect Conservation

Note: All packages used in the Rmarkdown script are automatically downloaded if needed, and unpacked, in the first (hidden) code chunk.

To ensure the models were optimzed for each species, species-specific models were used. The example here is shown for the first species, Acanthocinus aedilis. The data comes from modified files of individuals counts of each species were years, harbors, or blends with zero findings have been removed to ease model fit.

Step 1: Load the data. Repeat this step for each species. In this example, the file is expected to be inside a folder named “Data” within the project directory. If your copy of the repository does not contain this folder, or if the file is stored elsewhere, you may need to adjust the path accordingly.

A.aedilis <- read.csv(here("Data","GLMM_Species_Modified", "Modified_GLMM_A.aedilis.csv"), 
                        sep=";", header=TRUE, as.is=TRUE)

A.aedilis <- A.aedilis %>% 
  mutate_at(c("Year", "Harbor", "Blend", "Date"), as.factor) %>% 
  mutate_at(c("Year", "Traps"), as.numeric)

A.aedilis_pooled <- A.aedilis %>%
  group_by(Year, Blend, Harbor) %>%
  summarize(
    Count = sum(Number, na.rm = TRUE),  # Sum beetle counts within year
    traps = first(Traps),                       # Keep the same trap count for each year
    .groups = 'drop'
  )

A.aedilis_pooled$ID <- 1:nrow(A.aedilis_pooled)

Step 2: Fit a model and check diagnostics. Remember to remove the factor of Harbor from the model for species that only appeared at one harbor.

glmm_TMB <- glmmTMB(Count ~ Blend + Harbor + (1|Year) +
                   (1|ID) +   #Add if overdispersion occurs
                   offset(log(traps)), 
                   data = A.aedilis_pooled, 
                   family = poisson(link = 'log'))

check_overdispersion(glmm_TMB)
## # Overdispersion test
## 
##        dispersion ratio =  0.328
##   Pearson's Chi-Squared = 15.417
##                 p-value =      1
anova_output <- car::Anova(glmm_TMB, type = c('3'))
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Count
##               Chisq Df Pr(>Chisq)    
## (Intercept)  8.3171  1  0.0039273 ** 
## Blend       15.5159  2  0.0004273 ***
## Harbor      27.4242  2  1.109e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

If the effect of blend is significiant, move on to make pairwise-analyses between the blends:

emm <- emmeans(glmm_TMB, specs = ~ Blend, 
               type = "response", adjust = "none")
##  Blend       rate    SE  df asymp.LCL asymp.UCL
##  Ips        2.752 0.933 Inf     1.417     5.348
##  Monochamus 0.309 0.166 Inf     0.108     0.888
##  Pissodes   0.665 0.317 Inf     0.262     1.691
## 
## Results are averaged over the levels of: Harbor 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale
##  contrast              ratio    SE  df null z.ratio p.value
##  Ips / Monochamus      8.896 5.245 Inf    1   3.707  0.0006
##  Ips / Pissodes        4.138 2.221 Inf    1   2.646  0.0222
##  Monochamus / Pissodes 0.465 0.292 Inf    1  -1.219  0.4418
## 
## Results are averaged over the levels of: Harbor 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale

If you are working with a species with zero findings for one blend, use a contrast to compare the blends with findings against a null hypothesis of no difference from zero. For visualisation, all three blends are here tested against the null hypothesis.

contrast_results <- contrast(emm, list("Canada vs zero" = c (1, 0, 0),
  "Pissodes vs zero" = c(0, 1, 0),
                                       "Spain vs zero" = c(0, 0, 1)),
                             adjust = "none", type="response")
##  contrast         estimate    SE  df z.ratio p.value
##  Canada vs zero      1.012 0.339 Inf   2.987  0.0028
##  Pissodes vs zero   -1.173 0.538 Inf  -2.181  0.0292
##  Spain vs zero      -0.408 0.476 Inf  -0.856  0.3919
## 
## Results are averaged over the levels of: Harbor 
## Note: contrasts are still on the log scale

Step 3: Generate predictions of mean and coefficient of variance (CV) using the fitted model

A.aedilis_pooled$predicted_beetles <- predict(glmm_TMB, type = "response")

# Calculate predicted counts per trap (total predicted / number of # traps)
A.aedilis_pooled$predicted_beetles_per_trap <- A.aedilis_pooled$predicted_beetles / A.aedilis_pooled$traps

# Calculate mean and SD of residuals for each Blend
A.aedilis_pooled %>%
  group_by(Blend) %>%
  summarize(
    mean_beetles = mean(predicted_beetles_per_trap, na.rm = TRUE),  # Mean of predicted beetles
    cv_beetles = sd(predicted_beetles_per_trap, na.rm = TRUE) / mean(predicted_beetles_per_trap, na.rm = TRUE),  # CV = SD / Mean
    .groups = 'drop'
  )
## # A tibble: 3 × 3
##   Blend      mean_beetles cv_beetles
##   <fct>             <dbl>      <dbl>
## 1 Ips               3.73        1.63
## 2 Monochamus        0.672       2.79
## 3 Pissodes          1.33        1.97

References:

## @Book{,
##   title = {An {R} Companion to Applied Regression},
##   edition = {Third},
##   author = {John Fox and Sanford Weisberg},
##   year = {2019},
##   publisher = {Sage},
##   address = {Thousand Oaks {CA}},
##   url = {https://socialsciences.mcmaster.ca/jfox/Books/Companion/},
## }
## @Article{,
##   title = {Welcome to the {tidyverse}},
##   author = {Hadley Wickham and Mara Averick and Jennifer Bryan and Winston Chang and Lucy D'Agostino McGowan and Romain François and Garrett Grolemund and Alex Hayes and Lionel Henry and Jim Hester and Max Kuhn and Thomas Lin Pedersen and Evan Miller and Stephan Milton Bache and Kirill Müller and Jeroen Ooms and David Robinson and Dana Paige Seidel and Vitalie Spinu and Kohske Takahashi and Davis Vaughan and Claus Wilke and Kara Woo and Hiroaki Yutani},
##   year = {2019},
##   journal = {Journal of Open Source Software},
##   volume = {4},
##   number = {43},
##   pages = {1686},
##   doi = {10.21105/joss.01686},
## }
## @Article{,
##   author = {Mollie E. Brooks and Kasper Kristensen and Koen J. {van Benthem} and Arni Magnusson and Casper W. Berg and Anders Nielsen and Hans J. Skaug and Martin Maechler and Benjamin M. Bolker},
##   title = {{glmmTMB} Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling},
##   year = {2017},
##   journal = {The R Journal},
##   doi = {10.32614/RJ-2017-066},
##   pages = {378--400},
##   volume = {9},
##   number = {2},
## }
## @Article{,
##   title = {{performance}: An {R} Package for Assessment, Comparison and Testing of Statistical Models},
##   author = {Daniel Lüdecke and Mattan S. Ben-Shachar and Indrajeet Patil and Philip Waggoner and Dominique Makowski},
##   year = {2021},
##   journal = {Journal of Open Source Software},
##   volume = {6},
##   number = {60},
##   pages = {3139},
##   doi = {10.21105/joss.03139},
## }
## @Manual{,
##   title = {emmeans: Estimated Marginal Means, aka Least-Squares Means},
##   author = {Russell V. Lenth},
##   year = {2024},
##   note = {R package version 1.10.1},
##   url = {https://CRAN.R-project.org/package=emmeans},
## }
## @Manual{,
##   title = {here: A Simpler Way to Find Your Files},
##   author = {Kirill Müller},
##   year = {2020},
##   note = {R package version 1.0.1},
##   url = {https://CRAN.R-project.org/package=here},
## }