R/generate-regional-trends.R

In BrandonEdwards/bbsBayes: Hierarchical Bayesian Analysis of North American BBS Data

#' Generate regional trends continent and strata and optionally for countries, states/provinces, or BCRs from analyses run on the stratifications that support these composite regions
#'
#' \code{generate_regional_trends} calculates the geometric mean annual changes in population size for composite regions.
#'
#' @param indices regional indices generated by \code{generate_regional_indices}
#' @param Min_year Minimum year to calculate trends from (e.g., 1970). Default is NULL, in which case the trend is calculated from the first year of the time-series of the supplied annual_indices file
#' @param Max_year Maximum year to calculate trends to (e.g., 2018). Default is NULL, in which case the trend is calculated up to the last year of the time-series of the supplied annual_indices file
#' @param quantiles vector of quantiles to be sampled from the posterior distribution Defaults to c(0.025,0.05,0.25,0.5,0.75,0.95,0.975)
#' @param slope Logical, if TRUE, calculates an alternative trend metric, the slope of a log-linear regression through the annual indices. Default FALSE
#' @param prob_decrease Optional vector of percent-change values to calculate the posterior probabilities that the population has decreased by at least this much (e.g., prob_decrease = c(50) would result in a calculation of the probability that the population has decreased by more than 50 percent over the period of the trend, i.e., less than half the population remains. Default is NULL, in which case no probability of decrease is calculated.
#' @param prob_increase Optional vector of percent-change values to calculate the posterior probabilities that the population has increased by at least this much (e.g., prob_increase = c(100) would result in a calculation of the probability that the population has incrased by more than 100 percent, i.e., doubled, over the period of the trend. Default is NULL, in which case no probability of increase is calculated.
#'
#'
#'
#' @return Dataframe with one row for each region included in indices object, and columns including:
#'   \item{Start_year}{first year of the trend}
#'   \item{End_year}{last year of the trend}
#'   \item{Region}{short name of the region}
#'   \item{Region_alt}{Long name for region}
#'   \item{Region_type}{Type of region including continental, national,Province_State,BCR, bcr_by_national, or stratum}
#'   \item{Strata_included}{Strata included in the trend and annual index calculations}
#'   \item{Strata_excluded}{Strata potentially excluded from the trend and annual index calculations because they have no observations of the species in the first part of the time series}
#'   \item{Trend}{Estimated mean annual percent change over the trend time-period (i.e., Start_year - End_year), according to an endpoint comparison of annual index in Start_year and the annual index in End_year}
#'   \item{Trend_Q_quantiles}{quantiles of the posterior distribution of Trend estimates, matching levels included in the quantiles argument}
#'   \item{Percent_Change}{Estimated total percent change over the trend time-period}
#'   \item{Percent_Change_Q_quantiles}{quantiles of the posterior distribution of Percent Change estimates, matching levels included in the quantiles argument}
#'   \item{Slope_Trend}{Estimated mean annual percent change over the trend time-period, according to the slope of a linear regression through the log-transformed annual indices}
#'   \item{Slope_Trend_Q_quantiles}{quantiles of the posterior distribution of Percent Change estimates, matching levels included in the quantiles argument}
#'   \item{prob_decrease_X_percent}{proportion of the posterior distribution of Percent_Change that is below the percentage values supplied in prob_decrease}
#'   \item{prob_increase_X_percent}{proportion of the posterior distribution of Percent_Change that is above the percentage values supplied in prob_increase}
#'   \item{Relative_Abundance}{Mean of the annual index values across all years. An estimate of the average relative abundance of the species in the region. Can be interepreted as the predicted average count of the species in an average year on an average route by an average observer, for the years, routes, and observers in the existing data}
#'   \item{Observed_Relative_Abundance}{Mean of the 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. 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}
#'   \item{Number_of_Strata}{The number of strata included in the region}
#'   \item{Width_of_X_percent_Credible_Interval}{Width (in percent/year) of the credible interval on the Trend calculation. Calculated for the widest credible interval requested in quantiles argument. Default is 95 percent CI (i.e., Trend_Q0.975 - Trend_Q0.025)}
#'   \item{Width_of_X_percent_Credible_Interval_Slope}{Width (in percent/year) of the credible interval on the Trend calculation for the slope-based trend. Calculated for the widest credible interval requested in quantiles argument. Default is 95 percent CI (i.e., Slope_Trend_Q0.975 - Slope_Trend_Q0.025)}
#'   \item{Number_of_Routes}{The number of unique BBS routes included in the trend calculation for this region and species}
#'   \item{Mean_Number_of_Routes}{The average number of BBS routes across years contributing data for this region and species}
#'   \item{backcast_flag}{approximate proportion of the included species range*years that are free of extrapolated population trajectories e.g., 1.0 = data cover full time-series, 0.75 = data cover 75 percent of time-series. Only calculated if max_backcast != NULL}
#'
#' @importFrom stats quantile
#' @importFrom stringr str_split
#'
#'
#' @name bbsBayes-deprecated
#' @seealso \code{\link{generate_regional_trends}}
#' @keywords internal
NULL

#' @rdname bbsBayes-deprecated
#' @section \code{generate_regional_trends}:
#'   For \code{generate_regional_trends()}, use
#'   \code{generate_trends()}.
#'
#' @export
#'

generate_regional_trends <- function(indices = NULL,
                                Min_year = NULL,
                                Max_year = NULL,
                                quantiles = c(0.025,0.05,0.25,0.75,0.95,0.975),
                                slope = FALSE,
                                prob_decrease = NULL,
                                prob_increase = NULL)
{
  .Deprecated(new = "generate_trends",
              msg = paste("generate_regional_trends is deprecated in favour of generate_trends()"))

  if (is.null(indices))
  {
    stop("No indices supplied to generate_regional_trends()."); return(NULL)
  }
  n_all = indices$samples

  if (is.null(Min_year))
  {
    min_year = indices$y_min
    Min_year <- indices$startyear
    minyn <- 1
  }else{
    min_year <- indices$y_min + (Min_year-indices$startyear)
    minyn <- 1 + (Min_year-indices$startyear)

    if(min_year < 0){
      min_year = indices$y_min
      Min_year <- indices$startyear
      minyn <- 1
    }
  }
  if (is.null(Max_year))
  {
    max_year = indices$y_max
    Max_year = indices$startyear+(max_year - indices$y_min)
    maxyn <- 1+(indices$y_max-indices$y_min)
  }else{
    max_year <- indices$y_min + (Max_year-indices$startyear)
    maxyn <- 1+(Max_year-indices$startyear)
  }


regions = indices$regions
area_weights <- indices$area_weights

all_regions = names(n_all)
dsum = indices$data_summary


trend <- data.frame(Start_year = integer(),
                    End_year = integer(),
                    Region = character(),
                    Region_alt = character(),
                    Region_type = character(),
                    Strata_included = character(),
                    Strata_excluded = character(),
                    Trend = double(),
                    stringsAsFactors = FALSE)
for(qq in quantiles){
  trend[,paste0("Trend_Q",qq)] <- double()
}
trend[,"Percent_Change"] <- double()
for(qq in quantiles){
  trend[,paste0("Percent_Change_Q",qq)] <- double()
}



for(rr in regions){ #selecting the type of composite region

  #regsest = all_regions[grep(all_regions,pattern = rr)]
  regsest = paste0(rr,"_",unique(dsum[which(dsum$Region_type == rr),"Region"]))

  for(rrs in regsest){

    reg = gsub(rrs,pattern = paste0(rr,"_"),replacement = "",fixed = TRUE)
    w_summary_rows = which(dsum$Region == reg & dsum$Region_type == rr)


n = n_all[[rrs]]

reg_alt = unique(dsum[w_summary_rows,"Region_alt"])
st_inc = unique(dsum[w_summary_rows,"Strata_included"])
nstr = length(unlist(stringr::str_split(st_inc,pattern = " ; ")))
st_exc = unique(dsum[w_summary_rows,"Strata_excluded"])

  if(slope){

    bsl = function(i){
      n = length(wy)
      sy = sum(i)
      sx = sum(wy)
      ssx = sum(wy^2)
      sxy = sum(i*wy)
      b = (n*sxy - sx*sy)/(n*ssx - sx^2)
      return(b)
    }
    wy = c(minyn:maxyn)
    ne = log(n[,wy])
    m =  t(apply(ne,1,FUN = bsl))

    sl.t = as.vector((exp(m)-1)*100)

  }

  ch = n[,maxyn]/n[,minyn]
  tr = 100*((ch^(1/(maxyn-minyn)))-1)

  trendt <- data.frame(Start_year = (indices$startyear+minyn)-1,
                      End_year = (indices$startyear+maxyn)-1,
                      Region = reg,
                      Region_alt = reg_alt,
                      Region_type = rr,
                      Strata_included = st_inc,
                      Strata_excluded = st_exc,
                      Trend = median(tr),
                      stringsAsFactors = FALSE)
  for(qq in quantiles){
    trendt[,paste0("Trend_Q",qq)] <- stats::quantile(tr,qq,names = FALSE)
  }

  ### estimated %change
  trendt[,"Percent_Change"] <- 100*(median(ch)-1)
  for(qq in quantiles){
    trendt[,paste0("Percent_Change_Q",qq)] <- 100*(stats::quantile(ch,qq,names = FALSE)-1)
  }

  ### optional slope based trends
  if(slope){
    trendt[,"Slope_Trend"] <- median(sl.t)
    for(qq in quantiles){
      trendt[,paste0("Slope_Trend_Q",qq)] <- stats::quantile(sl.t,qq,names = FALSE)
    }
  }

  #### model conditional probabilities of population change during trend period
  if(!is.null(prob_decrease)){
    pch = 100*(ch-1)
for(pp in prob_decrease){
  trendt[,paste0("prob_decrease_",pp,"_percent")] <- length(pch[which(pch < (-1*pp))])/length(pch)
}
  }
  if(!is.null(prob_increase)){
    pch = 100*(ch-1)
    for(pp in prob_increase){
      trendt[,paste0("prob_increase_",pp,"_percent")] <- length(pch[which(pch > (pp))])/length(pch)
    }
  }


  ###### reliability criteria
  trendt[,"Relative_Abundance"] <- mean(dsum[w_summary_rows,"Index"])
  trendt[,"Observed_Relative_Abundance"] <- mean(dsum[w_summary_rows,"obs_mean"])
  trendt[,"Number_of_strata"] <- nstr
  q1 = quantiles[1]
  q2 = quantiles[length(quantiles)]
  trendt[,paste0("Width_of_",(q2-q1)*100,"_percent_Credible_Interval")] <- trendt[,paste0("Trend_Q",q2)]-trendt[,paste0("Trend_Q",q1)]
if(slope){
  trendt[,paste0("Width_of_",(q2-q1)*100,"_percent_Credible_Interval_Slope")] <- trendt[,paste0("Slope_Trend_Q",q2)]-trendt[,paste0("Slope_Trend_Q",q1)]
}
  trendt[,"Number_of_Routes"] <- mean(dsum[w_summary_rows,"nrts_total"])
  trendt[,"Mean_Number_of_Routes"] <- mean(dsum[w_summary_rows,"nrts"])
  trendt[,"backcast_flag"] <- mean(dsum[w_summary_rows,"backcast_flag"])


  trend = rbind(trend,trendt)
  }
}
  return(trend)
}