R/generate-strata-trends.R

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

#' Generate a dataframe of the stratum trends
#'
#' \code{generate_strata_trends} creates a data frame of trends by year,
#'   factored by each stratum. These trends are calculated using the indices
#'   generated by \code{generate_strata_indices}.
#'
#' @param indices Stratum indices generated by \code{generate_strata_indices}
#' @param min_year Minimum year to calculate trends from
#' @param max_year Maximum year to calculate trends to
#' @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
#'
#' @return Data frame of 2 variables:
#'   \item{Stratum}{Name of the stratum}
#'   \item{Trend}{Percent change from the minimum year to maximum year}
#'
#'
#' @importFrom stats lm
#'
#' @name bbsBayes-deprecated
#' @seealso \code{\link{generate_strata_trends}}
#' @keywords internal
NULL

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

generate_strata_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)
{
  .Deprecated(new = "generate_trends",
              msg = paste("generate_strata_trends is deprecated in favour of generate_trends()"))


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

  strata_indices <- indices$area_weights$num
  area_weights <- indices$area_weights
  if (is.null(min_year))
  {
    min_year = indices$y_min
    }
  if (is.null(max_year))
  {
    max_year = indices$y_max
    }




  trend <- data.frame(Start_year = integer(),
                      End_year = integer(),
                      Region = 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 (i in strata_indices)
  {
    if(slope){

      wy = c(min_year:max_year)

      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)
      }

     ne = log(n[,i,wy])
      m =  t(apply(ne,1,FUN = bsl))

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

     }


    ch = n[,i,max_year]/n[,i,min_year]
    tr = 100*((ch^(1/(max_year-min_year)))-1)

    trend_strata <- data.frame(Start_year = (indices$startyear+min_year)-1,
                        End_year = (indices$startyear+max_year)-1,
                        Region = area_weights[which(area_weights$num == i), ]$region,
                        Trend = median(tr),
                        stringsAsFactors = FALSE)
    for(qq in quantiles){
      trend_strata[,paste0("Trend_Q",qq)] <- quantile(tr,qq,names = FALSE)
    }
    trend_strata[,"Percent_Change"] <- 100*(median(ch)-1)
    for(qq in quantiles){
      trend_strata[,paste0("Percent_Change_Q",qq)] <- 100*(quantile(ch,qq,names = FALSE)-1)
    }

    if(slope){
      trend_strata[,"Slope_Trend"] <- median(sl.t)
      for(qq in quantiles){
        trend_strata[,paste0("Slope_Trend_Q",qq)] <- quantile(sl.t,qq,names = FALSE)
      }
    }


    trend <- rbind(trend, trend_strata)
  }


  return(trend)
}