Regional annual indices of abundance

Calculate annual indices of relative abundance by year for different regions. These indices can then be used to plot population trajectories for the species, and to estimate trends.

Usage

generate_indices(
  model_output = NULL,
  quantiles = c(0.025, 0.05, 0.25, 0.75, 0.95, 0.975),
  regions = c("continent", "stratum"),
  regions_index = NULL,
  alternate_n = "n",
  gam_smooths = FALSE,
  start_year = NULL,
  max_backcast = NULL,
  drop_exclude = FALSE,
  hpdi = FALSE,
  quiet = FALSE
)

Arguments

model_output: List. Model output generated by run_model().
quantiles: Numeric. Vector of quantiles to be sampled from the posterior distribution. Default is c(0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975). Note that these quantiles will be used to create confidence interval bands in plot_indices() and by quantiles in generate_trends(), so make sure you specify the ones you want to use later.
regions: Character. Which region(s) to summarize and calculate indices for. Default is "continent" and "stratum". Options also include "country", "prov_state", "bcr", and "bcr_by_country". Note that some regions only apply to specific stratifications. You can also supply a custom region that exists as a column in the regions_index data frame (see examples for more details).
regions_index: Data frame. Custom regions to summarize. Data frame must include all strata in the original data in one column (strata_name), and any custom regions defined as categories in other columns. See examples for more details.
alternate_n: Character. Indicating the name of the alternative annual index parameter in a model, Default is "n" which for all models represents an index of the estimated annual relative abundance, scaled as the expected mean count averaged across all BBS routes and observers in a given stratum and year. For some of the models included in bbsBayes2, alternatives exist that provide a partial decomposition of the time-series. For the "gamye" models, the parameter "n_smooth" represents the smooth-only version of the annual index of relative abundance (i.e., the component of the annual index estimated by the spline-based smooth of the GAM). This "n_smooth" is identical to the "n" values for the same model, but excludes the annual fluctuations. For the "gamye" models, this "n_smooth" parameter is likely the most natural parameter to use in estimating trends. A similar option exists for the "slope" models, where the parameter "n_slope" represents the component of the population trajectory estimated by the log-linear regression slope parameters in the model. Users should be particularly cautious about interpreting this "n_slope" values for relatively long time-series. As a continuous regression slope, it assumes interpreting it as an estimate of population trajectory and using it to generate trend estimates assumes that there is a single continuous rate of population change across the entire time-series. Biologically, this may be reasonable for 10-20 year periods, but will be less reasonable for longer time-periods.
gam_smooths: Logical. Whether to return an array of posterior samples of gam-based smooths through the full time-series of annual indices. These can by treated as smooth population trajectories from any model. NOTE: these gam smooths require fitting mgcv::gam() to every posterior draw of every possible population trajectory. As a result, generate_indices can require can require many minutes to complete if set to TRUE.
start_year: Numeric. Trim the data record before calculating annual indices.
max_backcast: Numeric. The number of years to back cast stratum-level estimates before the first year that species was observed on any route in that stratum. Default is NULL, which generates annual indices for the entire time series and ignores back-casting. CWS national estimates use a back cast of 5. Note that unless drop_exclude = TRUE, problematic years are only flagged, not omitted. See Details for more specifics.
drop_exclude: Logical. Whether or not strata that exceed the max_backcast threshold should be excluded from the calculations. Default is FALSE (regions are flagged and listed but not dropped).
hpdi: Logical. Should credible intervals and limits be calculated using highest posterior density intervals instead of simple quantiles of the posterior distribution. Default is FALSE. these intervals are often a better descriptor of skewed posterior distributions, such as the predicted mean counts that the indices represent. Note hpd intervals are not stable for small percentages of the posterior distribution, and so hpdi = TRUE is ignored for quantiles values between 0.33 and 0.67 (i.e., if the quantiles value defines a limit for a centered hpd interval that would include < 33% of the posterior distribution).
quiet: Logical. Suppress progress messages? Default FALSE.

Value

A list containing

indices - data frame of calculated regional annual indices of abundances
samples - array of posterior draws from the model
meta_data - meta data defining the analysis
meta_strata - data frame listing strata meta data
raw_data - data frame of summarized counts
gam_smooth_samples - optional-array of posterior draws of the gam smooth

indices contains the following columns:

year - Year of particular index
region - Region name
region_type - Type of region
strata_included - Strata potentially included in the annual index calculations
strata_excluded - Strata potentially excluded from the annual index calculations because they have no observations of the species in the first part of the time series, see arguments max_backcast and start_year
index - Strata-weighted count index (median)
index_q_XXX - Strata-weighted count index (by different quantiles)
obs_mean - Mean observed annual counts of birds across all routes and all years. An alternative estimate of the average relative abundance of the species in the region and year. Differences between this and the annual indices are a function of the model. For composite regions (i.e., anything other than stratum-level estimates) this average count is calculated as an area-weighted average across all strata included
n_routes - Number of BBS routes that contributed data for this species, region, and year
n_routes_total - Number of BBS routes that contributed data for this species and region for all years in the selected time-series, i.e., all years since start_year
n_non_zero - Number of BBS routes on which this species was observed (i.e., count is > 0) in this region and year
backcast_flag - Approximate annual average proportion of the covered species range that is free of extrapolated population trajectories. e.g., if 1.0, data cover full time-series; if 0.75, data cover 75 percent of time-series. Only calculated if max_backcast != NULL.

Details

max_backcast is a way to deal with the fact that the species of interest may not appear in the data until several years after the start of the time-series max_backcast specifies how many years can occur before the stratum is flagged. A max_backcast of 5 will flag any stratum without a non-zero (or non-NA) observation within the first 5 years of the time- series. Note that records are only flagged unless drop_exclude = TRUE. If you find that the early data record is sparse and results in the exclusion of many strata, consider trimming the early years by specifying a start_year.

Examples

# Using the example model for Pacific Wrens

# Generate the continental and stratum indices
i <- generate_indices(pacific_wren_model)
#> Processing region continent
#> Processing region stratum

# Generate the continental and stratum indices using hpdi
i <- generate_indices(pacific_wren_model, hpdi = TRUE)
#> Processing region continent
#> Processing region stratum

# Generate only country indices
i_nat <- generate_indices(pacific_wren_model, regions = "country")
#> Processing region country

# Use a custom region specification (dummy example)
library(dplyr)
#> 
#> Attaching package: ‘dplyr’
#> The following objects are masked from ‘package:stats’:
#> 
#>     filter, lag
#> The following objects are masked from ‘package:base’:
#> 
#>     intersect, setdiff, setequal, union
ri <- bbs_strata[["bbs_cws"]]
ri <- mutate(ri, my_region = if_else(prov_state %in% "ON",
                                     "Ontario", "Rest"))

# Generate indices with these custom regions
i_custom <- generate_indices(
  pacific_wren_model,
  regions = c("country", "prov_state", "my_region"),
  regions_index = ri)
#> Processing region country
#> Processing region prov_state
#> Processing region my_region

Usage

Arguments

Value

Details

See also

Examples