R/lppd.R

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

#' Calculate log posterior predictive density
#'
#' \code{lppd} Calculate log posterior preditive density (LPPD) for the supplied
#'   model.
#'
#'   NOTE: in order to calculated LPPD, the model MUST track the parameter
#'   "lambda". In species that are data-rich, such as Wood Thrush,
#'   this produces extremely large JAGS objects, and takes up a considerable
#'   amount of memory when simulating with \code{run_model}
#'
#' @param jags_data Data prepared by \code{prepare_jags_data}, used
#'   for input to the JAGS model
#' @param jags_mod JAGS list generated by \code{run_model}
#' @param pointwise If set to \code{TRUE}, a data frame is returned
#'   that contains the pointwise LPPD for each count. Defaults
#'   to \code{FALSE}
#'
#' @return Data frame of pointwise LPPD by count if \code{pointwise}
#'   is set to \code{TRUE}. Double precision numerical value of LPPD if
#'   \code{pointwise} is set to \code{FALSE}.
#'
#' @export
#'
#' @importFrom stats dpois
#'
#' @examples
#' # Toy example with Pacific Wren sample data
#' # First, stratify the sample data
#'
#' strat_data <- stratify(by = "bbs_cws", sample_data = TRUE)
#'
#' # Prepare the stratified data for use in a JAGS model.
#' jags_data <- prepare_jags_data(strat_data = strat_data,
#'                                species_to_run = "Pacific Wren",
#'                                model = "firstdiff",
#'                                min_year = 2014,
#'                                max_year = 2018)
#'
#' # Now run a JAGS model. Make sure to track the lambda parameter here
#'
#' jags_mod <- run_model(jags_data = jags_data,
#'                       n_adapt = 0,
#'                       n_burnin = 0,
#'                       n_iter = 5,
#'                       n_thin = 1,
#'                       parameters_to_save = c("n",
#'                                              "lambda"))
#'
#' # Output LPPD
#' lppd(jags_data = jags_data,
#'      jags_mod = jags_mod)
#'

lppd <- function(jags_data = NULL,
                 jags_mod = NULL,
                 pointwise = FALSE)
{
  bugs <- get_prepared_data(jags_data)

  lppd_point <- data.frame(index = 1:nrow(bugs))
  lppd_v <- vector(length = nrow(bugs))

  for(i in 1:nrow(bugs))
  {
    lppd_v[i] <- log(mean(stats::dpois(bugs[i,"Count"], jags_mod$sims.list$lambda[,i])))
  }

  lppd_point[,"lppd_point"] <- lppd_v

  if (isTRUE(pointwise))
  {
    return(lppd_point)
  }else{
    return(sum(lppd_point$lppd_point))
  }
}