R/generate-cont-indices.R

Defines functions generate_cont_indices

Documented in generate_cont_indices

#' Generate a dataframe of the continental yearly indices
#'
#' \code{generate_cont_indices} creates a data frame of the strata-weighed
#'   continental indices by year. This data frame can then be used to
#'   generate a population trajectory of the species.
#'
#' @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_cont_indices}}
#' @keywords internal
NULL

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

generate_cont_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 = "generate_cont_indices is deprecated in favour of generate_indices(regions = \"continental\")")

  if (is.null(jags_mod))
  {
    stop("No model output supplied to generate_cont_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_weight <- n

  obs_df = data.frame(year = integer(),
                      strat = integer(),
                      obs_mean = double())
  for (j in 1:n_strata)
  {
  rawst <- raw[which(raw$strat == j),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),]

  o_mns <- as.numeric(by(rawst[,2],INDICES = rawst[,1],FUN = mean,na.rm = TRUE))

  obs_df_t <- data.frame(year = c(y_min:y_max),
                        strat = j,
                        obs_mean = o_mns*((area_weights$area_sq_km[which(area_weights$num == j)])/ sum(area_weights$area_sq_km))*(non_zero_weight[j]))

    obs_df <- rbind(obs_df,obs_df_t)
  }

  # Weight each sampled n
  for (i in 1:n_samples)
  {
    for (j in 1:n_strata)
    {
      n_weight[i,j,] <- (n_weight[i,j,] * area_weights$area_sq_km[which(area_weights$num == j)])/ sum(area_weights$area_sq_km)
    }
  }

  N = apply(n_weight, c(1,3),sum)


  n_median <- apply(N, 2, median)
  data_summary <- data.frame(Year = seq(y_min:y_max),
                             Region = "Continental",
                             Region_alt = "Continental",
                             Region_type = "continental",
                             Index = n_median)
  for(qq in quantiles){
  data_summary[,paste0("Index_q_",qq)] <- apply(N,2,stats::quantile,probs = qq)
  }

  data_summary$Year <- as.integer((data_summary$Year - 1) + min(r_year))
  data_summary$obs_mean <- as.numeric(by(obs_df[,3],INDICES = obs_df[,1],FUN = sum,na.rm = TRUE))


  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))
}
BrandonEdwards/bbsBayes documentation built on March 3, 2023, 9:55 a.m.