R/sensitivity_to_prior.R

Defines functions .extract_priors_rstanarm .prior_new_location .extract_parameters sensitivity_to_prior.default sensitivity_to_prior.stanreg sensitivity_to_prior

Documented in .extract_priors_rstanarm .prior_new_location sensitivity_to_prior sensitivity_to_prior.stanreg

#' Sensitivity to Prior
#'
#' Computes the sensitivity to priors specification. This represents the
#' proportion of change in some indices when the model is fitted with an
#' antagonistic prior (a prior of same shape located on the opposite of the
#' effect).
#'
#' @param model A Bayesian model (`stanreg` or `brmsfit`).
#' @param index The indices from which to compute the sensitivity. Can be one or
#'   multiple names of the columns returned by `describe_posterior`. The case is
#'   important here (e.g., write 'Median' instead of 'median').
#' @param magnitude This represent the magnitude by which to shift the
#'   antagonistic prior (to test the sensitivity). For instance, a magnitude of
#'   10 (default) means that the mode wil be updated with a prior located at 10
#'   standard deviations from its original location.
#' @param ... Arguments passed to or from other methods.
#'
#' @examplesIf require("rstanarm")
#' \donttest{
#' library(bayestestR)
#'
#' # rstanarm models
#' # -----------------------------------------------
#' model <- rstanarm::stan_glm(mpg ~ wt, data = mtcars)
#' sensitivity_to_prior(model)
#'
#' model <- rstanarm::stan_glm(mpg ~ wt + cyl, data = mtcars)
#' sensitivity_to_prior(model, index = c("Median", "MAP"))
#' }
#' @seealso DescTools
#' @export
sensitivity_to_prior <- function(model, ...) {
  UseMethod("sensitivity_to_prior")
}


#' @rdname sensitivity_to_prior
#' @export
sensitivity_to_prior.stanreg <- function(model, index = "Median", magnitude = 10, ...) {
  # Original
  params <- .extract_parameters(model, index = index, ...)

  # Priors
  priors <- .extract_priors_rstanarm(model)
  new_priors <- .prior_new_location(prior = priors$prior, sign = sign(params$Median), magnitude = magnitude)
  model_updated <- stats::update(model, data = insight::get_data(model), prior = new_priors, refresh = 0)

  # New model
  params_updated <- .extract_parameters(model_updated, index = index, ...)

  # Compute index
  sensitivity <- abs(as.matrix(params_updated[-1]) - as.matrix(params[-1])) / abs(as.matrix(params[-1]))

  # Clean up
  sensitivity <- as.data.frame(sensitivity)
  names(sensitivity) <- paste0("Sensitivity_", names(params_updated)[-1])
  sensitivity <- cbind(params_updated[1], sensitivity)
  row.names(sensitivity) <- NULL
  sensitivity
}


#' @export
sensitivity_to_prior.default <- function(model, ...) {
  insight::format_error(sprintf("Models of class '%s' are not yet supported.", class(model)[1]))
}


#' @keywords internal
.extract_parameters <- function(model, index = "Median", ...) {
  # Handle BF
  test <- c("pd", "rope", "p_map")
  if (any(c("bf", "bayesfactor", "bayes_factor") %in% index)) {
    test <- c(test, "bf")
  }

  params <- suppressMessages(describe_posterior(
    model,
    centrality = "all",
    dispersion = TRUE,
    test = test,
    ...
  ))

  params <- params[params$Parameter != "(Intercept)", ]
  params[unique(c("Parameter", "Median", index))]
}




#' Set a new location for a prior
#' @keywords internal
.prior_new_location <- function(prior, sign, magnitude = 10) {
  prior$location <- -1 * sign * magnitude * prior$scale
  prior
}





#' Extract and Returns the priors formatted for rstanarm
#' @keywords internal
.extract_priors_rstanarm <- function(model, ...) {
  priors <- rstanarm::prior_summary(model)

  # Deal with adjusted scale
  if (!is.null(priors$prior$adjusted_scale)) {
    priors$prior$scale <- priors$prior$adjusted_scale
    priors$prior$adjusted_scale <- NULL
  }
  priors$prior$autoscale <- FALSE


  priors
}
DominiqueMakowski/bayestestR documentation built on Feb. 12, 2024, 8:27 a.m.