R/jags_check.R

In jdtrat/mickjaggr: Easily Roll Through Bayesian Regressions with 'JAGS'

#' Get residuals
#'
#' @param .recipe A recipe object
#' @param .model The output of \code{\link{jags_model_run}}
#'
#' @return A data frame containing two columns `.pred`, which are the estimated
#'   predictor values and `.resid`, the difference between the true data and the
#'   predicted data.
#' @export
#'
#' @examples
#' # coming soon
jags_check_residuals <- function(.recipe, .model) {

  # Create an index column for all output coefficients that have beta parameters
  # colnames should be same for all chains, so we're just using first chain.
  coef_cols_index <- grepl("beta", colnames(.model[[1]]))

  # creaate a matrix with the means of each beta coef from each chain
  beta_chains <- vapply(.model[,coef_cols_index], colMeans, FUN.VALUE = numeric(sum(coef_cols_index)))
  # average across all chains (rows) for each beta coefficient
  betas <- apply(beta_chains, 1, mean)

  predictors <- jags_data_predictors(.recipe = .recipe)
  response <- jags_data_response(.recipe = .recipe)

  y_hat <- predictors %*% betas

  residual_df <- data.frame(.preds = y_hat,
                            .resid = response - y_hat)

  return(residual_df)
}


#' Calculate Bayesian P-Values for each check
#'
#' The Bayesian P-Value, in this implementation, calculates the summary
#' statistics supplied by the user, such as "max" for the observed response
#' variable. As the JAGS code runs, a prediction is generated for every
#' observation using the estimated parameters.
#'
#' This code checks to see whether the actual summary statistic in the observed
#' response vector is greater than the summary statistic generated from the
#' predicted values. If the model was an appropriate choice for the data, the
#' actual summary statistic would be greater than approximately 50% of the
#' summary statistics generated by the predicted values. That is, the closer a
#' Bayesian P-Value is to 0.5, the better the fit. The closer the Bayesian
#' P-Value is to 0 or 1, the worse the fit.
#'
#'
#' @keywords internal
#' @param .list The output of \code{\link{jags_model_run}}
#' @param .check A character vector of summary stats to perform
#' @param .response_check The actual statistics from the response variable
#'
#' @return A numeric vector with the Bayesian P value for each summary statistic
#'   for each iteration.
#'
check_pp_individual <- function(.list, .check, .response_check) {
  vapply(.list[,paste0("pp.", .check)], function(x) mean(.response_check[.check] > x), FUN.VALUE = numeric(1))
}



#' Calculate Bayesian P-Values for each check
#'
#' @keywords internal
#' @param .list The output of \code{\link{jags_model_run}}
#' @param .check A character vector of summary stats to perform
#' @param .response_check The actual statistics from the response variable
#'
#' @return A matrix with the Bayesian P value for each summary stat for every iteration
#'
check_pp_internal <- function(.list, .check, .response_check) {
  vapply(.check, FUN = check_pp_individual, FUN.VALUE = numeric(length(.list)), .list = .list, .response_check = .response_check)
}

#' Perform Posterior Predictive Checks
#'
#' @inheritParams jags_check_residuals
#' @param checks A character string (or vector) of statistical summaries such as c("min", "max", "mean", "sd").
#' @param all_iters A logical indicating whether the statistic summaries should be displayed per iteration. FALSE by default.
#'
#' @return Returns a Bayesian P Value for each posterior predictive check.
#' @export
#'
#' @examples
#' # coming soon
jags_check_pp <- function(.recipe, .model, checks, by_iter = FALSE) {

  response <- jags_data_response(.recipe = .recipe)

  response_checks <- vapply(X = checks, FUN = do.call, FUN.VALUE = numeric(1), list(response))

  all_iters <- check_pp_internal(.list = .model,
                                 .check = checks,
                                 .response_check = response_checks)

  dimnames(all_iters)[[1]] <- paste0("iter", 1:length(.model))

  if (!by_iter) {
  output <- apply(all_iters, 2, mean)
  return(output)
  } else if (by_iter) {
    return(all_iters)
  }
}