EcoCountHelper: an R package and analytical pipeline for the analysis of ecological count data using GLMMs, and a case study of bats in Grand Teton National Park

GLMM workflow

We have adapted a workflow for constructing and interpreting GLMMs based on the ‘glmmTMB’ package (Bolker et al., 2012; Bolker, 2016; Brooks et al., 2017b). An outline of this adapted workflow involves the following steps for each taxonomic group of interest:

1. We first decide on predictor variables and random effects structures for a response of interest.

2. We then compare the fits of Poisson and negative binomial error distributions (as well as the potential for zero-inflation) by corroborating AIC values, mean-variance plots, and examining simulated residual test plots—the latter of which also help identify outliers and uniformity.

3. We then interpret model results through both standardized coefficient visualizations and predicted changes to response variables, given a user-specified increase in a single predictor variable.

To aid researchers, and particularly land and wildlife managers, in implementing GLMM-based analyses for count data, we developed the ‘EcoCountHelper’ R package to simplify this workflow.

Package workflow

The ‘EcoCountHelper’ package (downloadable via Cole, Gomes & Barber, 2021a, 2021b) is meant to assist researchers for the portions of a GLMM-based analysis after the point at which candidate models have been generated. It is important that individuals using ‘EcoCountHelper’ take care to carefully consider candidate model structures, and to appropriately implement GLMM components such as zero-inflation formula (Martin et al., 2005) and random effects (Harrison et al., 2018; Gomes, 2022). Despite the popularity of selecting predictor variables for both conditional model zero-inflation formulas via dredging with information criteria (e.g., AIC), there are many pitfalls with using these methods (Guthery et al., 2005; Link & Barker, 2006; Bolker, 2008). While there are multiple excellent resources for learning about GLMMs (Bolker, 2008), Harrison et al. (2018) is an excellent touchstone for these analyses. Once one has already generated appropriate candidate models with ‘glmmTMB’ (Bolker et al., 2012; Bolker, 2016; Brooks et al., 2017a) and prepared all data (see Supplemental “EcoCountHelperExample” vignette), the workflow we suggest for GLMM analyses using the ‘EcoCountHelper’ functions (Table 1) can be executed.

The first step in our pipeline is to select the best-fitting model for each response group under investigation, taxonomic (e.g., species, genus) or otherwise (e.g., indistinguishable grouped taxa, foraging guilds, operational taxonomic units). Two non-exclusive techniques for deciding whether to use zero-inflated models and selecting the best error distribution family for a particular count dataset are: (1) examine AIC values to assess zero-inflation and select an error distribution, and (2) examine the mean-variance relationship within the data and determine an error distribution family that best mirrors the mean-variance relationship. The ‘EcoCountHelper’ package has a function, “ModelCompare”, for simultaneously obtaining AIC values for each group’s candidate models and creating a vector of the models with the lowest AIC scores for each group. The “DistFit” family of functions aids in visual examinations of mean-variance relationships for each group in the analysis. These “DistFit” functions allow users to specify vectors by which data should be aggregated for examining the mean-variance relationship of the data. The “DistFit” functions then generate a scatterplot displaying the mean-variance relationship for the data and draw lines through the scatterplot showing three common error distributions used for count data: Poisson, negative binomial with a linear mean-variance relationship (a.k.a. quasi-Poisson), and negative binomial with a quadratic mean-variance relationship (variance increases by the squared mean) (Bolker et al., 2012). The functions used during the model construction process that correspond with each of these error distributions are “poisson” from the ‘stats’ package, “nbinom1” from ‘glmmTMB’, and “nbinom2” from ‘glmmTMB’, respectively. For the sake of simplicity, we will use these function names to describe the associated distributions from this point on. After generating mean-variance plots using the “DistFit” family of functions, the user can then choose the model with the most appropriate error distribution for each group by visually examining the plots and choosing the model employing the error distribution that best fits. By using both of these methods (AIC and mean-variance relationships) we generally arrive at the same conclusion about which distribution best fits the data. There is, however, some subjectivity in distribution choice when both Poisson and negative binomial models appear to fit similarly. In our experience with this situation, the differences in parameter estimates are usually negligible, and thus distributional choice will likely be inconsequential. In this case, one may prefer the simpler Poisson distribution because there are fewer parameters to estimate (Poisson is a special case of negative binomial model in which the variance equals the mean, so there is only one term to estimate instead of two). It is also worth noting that using AIC values to determine whether zero-inflation formulas result in better model fit does not test for zero-inflation of data. Rather, the use of AIC values to compare candidate models either employing or omitting zero-inflation formulas provides information about whether variables specified in zero-inflation formulas explain any zero-inflation that does exist within given data. Zero-inflation formulas may be implemented to examine the effects of one or multiple predictors on zero-inflation of data, but care must be taken to ensure that predictors specified in zero-inflated formulas could reasonably contribute to zero-inflation of data.

Following the model selection process, it’s important to test the goodness-of-fit of the chosen models. Because a model fits a dataset better than all others does not mean it is adequate for making inferences or predictions. Plotting residuals (or simulated residuals) is a quick and often useful step in model criticism, but it is important to note that there are many other ways to assess goodness-of-fit. EcoCountHelper’s “ResidPlot” family of functions provide a simple way to simulate residuals for models that employ a non-Gaussian error distribution and generate plots to visually check residual uniformity, outliers, and over-/under-dispersion (important indicators of goodness of fit). These functions borrow functions from the ‘DHARMa’ R package (Hartig, 2022). It is important to note that if precise predictions are required from a model in order to make management decisions, the “gold-standard” for assessing the predictive capacity of a model is leave-one-out (a variant of k-fold) cross validation (Hawkins, Basak & Mills, 2003; Vehtari, Gelman & Gabry, 2017). Holding out part of a dataset to assess predictive capacity may be a computationally less intense alternative (Kim, 2009). These methods involve additional computational steps and are outside the scope of this article.

Assuming selected models, for all species or groups of interest, fit the data adequately, the next step is interpretation of model results by visualizing model coefficients and confidence intervals. The “EffectsPlotter” function generates coefficient plots with up to three user-specified confidence intervals surrounding each coefficient value. While many managers and practitioners are most often familiar with p-values (and an alpha of 0.05), any confidence interval chosen is arbitrary and we suggest thinking deeply about the consequences of type I and type II errors in your system and adjust confidence intervals accordingly. See the “Model Interpretation” section of the methods for additional discussion surrounding confidence interval selection. More importantly the effect size, on the scale of the original response variable, should be taken into consideration when making management decisions (Sullivan & Feinn, 2012). Because the coefficients of fitted GLMMs (except for those using an identity link, as typically used in linear (Gaussian) models) are a product of a link transformation (via the link specified during the model fitting process; with Poisson and negative binomial this is typically a log link), the coefficients cannot be used to predict the effect without back-transformation. Additionally, continuous data used during model fitting is often standardized for computational gains, increased model convergence, and comparison purposes (Gelman, 2008) which further complicates the process of making predictions with meaningful units. To simplify this process, we created the “RealEffects” family of functions in ‘EcoCountHelper’ which allows users to rapidly assess the response to given changes in predictors using untransformed “real world” values and outputs to a readable text. A flowchart integrating both conceptual and programmatic steps in the workflow proposed here can be found in Fig. 1.

‘EcoCountHelper’ also has multiple accessory functions that do not directly pertain to the model fitting, model selection, and result interpretation process outlined above. One noteworthy accessory function is “scale2”. This function is identical to the base R scale function in that it standardizes a vector, but rather than subtracting the mean from each value and dividing by the standard deviation, “scale2” subtracts the mean and then divides each value by two standard deviations (Eq. (1)) as suggested by Gelman (2008). During the model fitting process, using the “scale2” function puts continuous values on the same scale as binary categorical variables allowing a direct comparison of standardized continuous coefficients and categorical coefficients.

(1) $${\displaystyle \frac{x-\underset{\_}{\overline{x}}}{2\sigma }}$$

Package limitations

The ‘EcoCountHelper’ package was designed to make GLMMs more accessible to researchers and managers who otherwise may not have sufficient programmatic skills to carry out analyses in a timely, reproducible, and statistically responsible manner. Because many of our package’s functions are designed to return specific results (e.g., plots, test statistics) while minimizing unguided data preparation, users must adhere to relatively strict object naming schemes to ensure that pattern recognizers within the functions can identify relevant objects in the global environment. This forces users to conform to a model naming scheme that provides information about the taxonomic group the model belongs to, and also creates intermediate objects that have little or no flexibility in naming conventions. While it makes some aspects of the package inflexible, this strict naming convention allows users unfamiliar with regular expressions to forgo the exercise of identifying objects in the global environment and adding them to a vector for further use either manually or through regular expressions. Because functions that generate plots require virtually no user input for plot construction, it can be difficult to make edits to resulting plots. Again, while this may be inconvenient for programmatically savvy individuals, this inflexibility in data visualization allows for programmatically inexperienced individuals to carry out GLMM based analyses in a simple and straightforward manner.

Case study methods

Throughout the summers (June to September) of 2016 and 2017, we monitored bat activity throughout Grand Teton National Park, Wyoming, USA (referred to as Grand Teton from here on) using passive acoustic monitors. Grand Teton is composed of a relatively flat, high elevation (average ~2,073 m (National Park Service, 2019)) valley bounded by the Teton mountain range to the west, and the Gros Ventre mountain range to the east. The majority of our work occurred within the valley of the park.

To assess the impact of buildings and artificial light on the bats of Grand Teton, we used Wildlife Acoustics SM2BAT and SM3BAT units to record bat echolocation at 36 sites throughout the park with varying levels of anthropogenic infrastructure. To obtain high quality recordings, we set sampling rates for all ultrasonic recorders to 384 kHz allowing us to capture frequencies up to 192 kHz. We programmed all recorders to allow ultrasonic recordings from thirty minutes before sunset to sunrise, and to begin recording when triggered by a 16 kHz or greater signal that was 18 dB or more above background sound levels. This triggered recording scheme generates WAV files with individual bat call sequences, allowing for rapid and automated processing using call identification software (e.g., SonoBat, Kaleidoscope).

We chose monitoring sites to capture variation in multiple landscape characteristics including elevation, distance to water, land cover, and human infrastructure that may influence habitat suitability for bats (e.g., porous buildings that could serve as day roosts (Voigt et al., 2016)), the presence of non-natural light sources (Stone, Harris & Jones, 2015). In total, we monitored bat activity at 36 sites throughout the park (see Fig. 2), deploying units at 26 sites in 2016, and 27 sites in 2017 (17 of which were sites established in 2016). We deployed each acoustic monitor for five to six nights at each site during 2016 and 13 to 14 nights in 2017 due to an increased availability of units. Our sampling effort totaled 840 site-nights for both years of data collection, with 276 site-nights in 2016 and 564 site-nights in 2017. For additional information regarding data collection and preparation, see the Supplemental Information.

Analysis

We used generalized linear mixed models (GLMMs) to assess the effect of the predictors mentioned above on activity for each species of bat that was detected during at least 50 site-nights (Eptesicus fuscus, Lasiurus cinereus, Lasionycteris noctivagans, Myotis evotis, Myotis lucifugus, Myotis volans and Myotis yumanensis). We ran all models in R using the ‘glmmTMB’ package (Brooks et al., 2017a) and standardized all continuous variables using our “scale2” function using the formula.

For each bat species under investigation, we generated six models (three different distribution families and a zero-inflated version of each; see below). All models shared a common conditional model structure that included year, lunar phase, ordinal date, elevation, roost-building density index, proportion of cool lights within 50 m, summed brightness score of lights within 50 m, distance to the nearest water source, proportion of developed land cover within 50 m, and proportion of forested area within 50 m as fixed effects, and monitoring site as a random intercept term. Additionally, all models used a log link function. For each species, two models were generated for each error distribution implemented in the mean-variance plots produced by the “DistFit” family of functions (“poisson”, “nbinom1”, and “nbinom2”): one with a zero-inflation formula, and one without. We modelled zero-inflation with ordinal date and site as predictors.

For each species, we determined the most appropriate error distribution by corroborating results from the synergistic “DistFit” family of functions and the “ModelCompare” function (see Package Workflow above). We determined a priori that in the case of conflicting results from the two processes above, we would rely on visual mean-variance plots for assessing a best-fit model for the proceeding analytical steps. We also used the AIC values generated by the “ModelCompare” function to determine whether a model with a zero-inflated component was necessary. Note that only the parameters specified in the zero-inflation formula are assessed for effects on zero-inflation, and while zero-inflation may exist in one’s data, ZIP or ZINB candidate models may not result in better model fit if predictors influencing zero-inflation are not specified in a zero-inflation formula. The best-fitting model for each species was subsequently tested for goodness-of-fit using EcoCountHelper’s “ResidPlot” functions. We also checked variance inflation factors (VIF) in R using the function “check_collinearity” in the ‘performance’ package (Lüdecke et al., 2020).

Model interpretation

We first examined the scaled coefficient estimates by plotting each estimate with confidence intervals on a common scale using EcoCountHelper’s “EffectPlotter” function. These plots were used to inform our understanding of the relative importance of landscape features to species-specific bat activity. While informative in terms of relative effect sizes and confidence levels surrounding those effects, these plots do not provide meaningful absolute effect sizes that can be interpreted in an ecological context. The purpose of examining effects plots of scaled and transformed coefficient estimates is to develop an understanding of the relative magnitude of effect each predictor of interest has on a bat species’ activity throughout the study area.

Rather than creating and examining effects plots, one might want to simply examine the p-values of coefficient estimates to determine which predictors warrant further investigation. We suggest avoiding making any decisions regarding the importance of landscape features for bat conservation using p-values alone. While there is nothing inherently wrong with this approach, relying on p-values alone is insufficient for meaningful interpretations of results (Sullivan & Feinn, 2012; Lin, Lucas & Shmueli, 2013; Halsey et al., 2015; Vidgen & Yasseri, 2016). By examining scaled coefficient estimates, one can quickly assess both the relative magnitude of, and confidence in, effect sizes. It is important for ‘EcoCountHelper’ users to determine confidence intervals and effect sizes that are suitable for their purposes. Many researchers have traditionally used 95% confidence intervals to determine whether the effect of a predictor is meaningful, but any confidence interval can be used to arbitrate the validity of an effect so long as the researcher has considered the consequences of type I errors due to a low confidence threshold. In interpreting the results of our case analysis, we examined both 85% and 95% confidence intervals. While not as conservative as 95% confidence intervals, 85% confidence intervals still provide information regarding trends in data that, given more sampling or follow-up research, could be confirmed at a higher confidence level (Cumming & Finch, 2005). Similarly, one may observe a large effect size but have insufficient data to obtain 95% confidence intervals that do not overlap zero. By tempering the interpretation of results by acknowledging the confidence surrounding those estimates, managers can act on the information available if necessary while also recognizing the need for additional research on the topic of interest. Considering the low risk posed by a type I error (i.e., persistence of historic and porous structures throughout the park) and the potential benefit of increasing our understanding of bat ecology in Grand Teton with relatively high certainty, we decided interpreting trends present at approximately the 85% confidence level was acceptable. Ultimately we suggest that researchers and managers implementing our suggested workflow do not treat it as a means of hypothesis testing, but rather as a path toward understanding a system surrounded by a level of uncertainty. Both the estimates provided by models as well as the associated uncertainty should be examined within the context of the research objectives to inform one’s conclusion about a predictor’s influence on an ecological process.

After identifying landscape features that influence activity levels, we examined the real-world effects of our predictors by calculating the factor by which bat activity changes relative to a percentage change in each predictor and associated confidence intervals. In the case of models implementing a (natural) log-link, this entails back-transforming the coefficient estimate by first unscaling and then exponentiating Euler’s number (e) by that quotient/value. Then the resulting term is exponentiated by the difference of interest in the values of the predictor. For example, if one is interested in the relative change across 50 m of elevation change, this difference ( $\delta \mathrm{x}$) would be 50 (since the rest of the term has already been back-transformed to the original units and represents 1 unit change; in this example 1 m). Because we scaled continuous variables using two standard deviations, we used Eq. (2) to calculate the factor by which the response variable changes given a change in a predictor.

(2) $$\left(e^{{\displaystyle \frac{\beta }{2\sigma }}}\right)\delta \mathrm{x}$$

In Eq. (2), β represents the model coefficient of interest, σ represents the standard deviation of the unstandardized predictor, e represents Euler’s number, and $\delta \mathrm{x}$ is the difference of interest in predictor x. To facilitate these calculations, the ‘EcoCountHelper’ package contains a “RealEffect” group of functions that back-transform estimates and return unscaled response for specified changes in predictors.