R/generate-strata-indices.R

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

#' Generate a dataframe of the stratum indices
#'
#' \code{generate_strata_indices} creates a data frame of indices by year,
#'   factored by each stratum. These indicies can be used to generate
#'   trajectory plots of the species for each strata.
#'
#' @param jags_mod JAGS list generated by \code{run_model}
#' @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 alternate_n text string indicating the name of the alternative approach to calculating an annual index, Default is "n"
#'
#' @return List of 6 objects:
#'   \item{data_summary}{dataframe with the following columns}
#'   \item{Year}{Year of particular index}
#'   \item{Region}{Region name}
#'   \item{Index}{Strata-weighted count index}
#'   \item{additional columns for each of the values in quantiles}{quantiles of the posterior distribution}
#'
#'   \item{samples}{array of all samples from the posterior distribution}
#'   \item{area-weights}{data frame of the strata names and area weights used to calculate the continental estimates}
#'   \item{y_min}{first year used in the model, scale 1:length of time-series}
#'   \item{y_max}{last year used in the model, scale 1:length of time-series}
#'   \item{startyear}{first year used in the model, scale 1966:2018}
#'
#' @importFrom stats quantile
#'
#' @name bbsBayes-deprecated
#' @seealso \code{\link{generate_strata_indices}}
#' @keywords internal
NULL

#' @rdname bbsBayes-deprecated
#' @section \code{generate_strata_indices}:
#'   For \code{generate_strata_indices()}, use
#'   \code{generate_indices(regions = "stratum")}.
#'
#' @export
#'

generate_strata_indices <- function(jags_mod = NULL,
                                  quantiles = c(0.025,0.05,0.25,0.75,0.95,0.975),
                                  alternate_n = "n")
{
  .Deprecated(new = "generate_indices",
              msg = paste("generate_strata_indices is deprecated in favour of generate_indices(region = \"stratum\")"))


  if (is.null(jags_mod))
  {
    stop("No model output supplied to generate_strata_indices()."); return(NULL)
  }

  data_list <- extract_index_data(jags_mod = jags_mod,alt_n = alternate_n)
  n <- data_list$n
  area_weights <- data_list$area_weights
  y_min = data_list$y_min
  y_max = data_list$y_max
  r_year = data_list$r_year
  bugs_data = data_list$bugs_data

  raw = data.frame(year = bugs_data$year,
                   count = bugs_data$count,
                   strat = bugs_data$strat)

  non_zero_weight = bugs_data$nonzeroweight

  n_samples <- dim(n)[1]
  n_strata <- dim(n)[2]
  n_years <- dim(n)[3]

  strata_indices <- area_weights$num

  n_median <- apply(n, c(2,3), median)


  # n_25 <- apply(n, c(2,3), stats::quantile, probs = 0.025)
  # n_975 <- apply(n, c(2,3), stats::quantile, probs = 0.975)

  data_summary <- data.frame(Year = integer(),
                             Index = double(),
                             Region = character(),
                             Region_alt = character(),
                             Region_type = character(),
                             stringsAsFactors = FALSE)

  for (i in strata_indices)
  {
    rawst = raw[which(raw$strat == i),c("year","count")]
    yrs = data.frame(year = c(y_min:y_max))
    rawst = merge(rawst,yrs,by = "year",all = TRUE)
    rawst = rawst[order(rawst$year),]


    strat_summary <- data.frame(Year = seq(y_min:y_max),
                               Index = as.numeric(as.vector(n_median[i,])),
                               Region = area_weights[which(area_weights$num == i), ]$region,
                               Region_alt = area_weights[which(area_weights$num == i), ]$region,
                               Region_type = "stratum")
    strat_summary$Year <- (strat_summary$Year - 1) + min(r_year)
    for(qq in quantiles){
      strat_summary[,paste0("Index_q_",qq)] <- apply(n[,i,],2,stats::quantile,probs = qq)
    }
    strat_summary$obs_mean = non_zero_weight[i]*as.numeric(by(rawst[,2],INDICES = rawst[,1],FUN = mean,na.rm = TRUE))

    data_summary <- rbind(data_summary, strat_summary)
  }

  data_summary$Year = as.integer(data_summary$Year)
  return(list(data_summary = data_summary,
              samples = n,
              area_weights = area_weights,
              y_min = y_min,
              y_max = y_max,
              startyear = min(r_year),
              raw_data = raw))
}