priorsense: Prior (and likelihood) diagnostics and sensitivity analysis

Description

The priorsense package provides functions for prior and likelihood sensitivity analysis of Bayesian models. Currently it implements methods to determine the sensitivity of the posterior to power-scaling perturbations of the prior and likelihood and is the first implementation of the method described in Kallioinen et al. (2023).

Details

The main diagnostic function provided by priorsense is powerscale_sensitivity. Given a fitted model or draws object, it computes the powerscaling sensitivity diagnostic described in Kallioinen et al. (2023). It does so by perturbing the prior and likelihood and computing the effect on the posterior, without needing to refit the model (using Pareto smoothed importance sampling and importance weighted moment matching; Vehtari et al. 2022, Paananen et al. 2021).

Author(s)

Maintainer: Noa Kallioinen noa.kallioinen@aalto.fi [copyright holder]

Authors:

• Topi Paananen

• Paul-Christian Bürkner

• Aki Vehtari

Other contributors:

• Frank Weber [contributor]

References

Kallioinen, N., Paananen, T., Bürkner, P-C., Vehtari, A. (2023). Detecting and diagnosing prior and likelihood sensitivity with power-scaling perturbations. Statistics and Computing. 34(57). doi:10.1007/s11222-023-10366-5

Vehtari, A., Simpson, D., Gelman, A., Yao, Y., and Gabry, J. (2024). Pareto smoothed importance sampling. Journal of Machine Learning Research. 25(72). https://jmlr.org/papers/v25/19-556.html

Paananen, T., Piironen, J., Bürkner, P-C., Vehtari, A. (2021). Implicitly adaptive importance sampling. Statistics and Computing. 31(16). doi:10.1007/s11222-020-09982-2

See Also

powerscale_sensitivity powerscale_sequence powerscale powerscale_plot_ecdf powerscale_plot_dens powerscale_plot_quantities

Cumulative Jensen-Shannon divergence

Description

Computes the cumulative Jensen-Shannon distance between two samples.

Usage

cjs_dist(
  x,
  y,
  x_weights = NULL,
  y_weights = NULL,
  metric = TRUE,
  unsigned = TRUE,
  ...
)

Arguments

Details

The Cumulative Jensen-Shannon distance is a symmetric metric based on the cumulative Jensen-Shannon divergence. The divergence CJS(P || Q) between two cumulative distribution functions P and Q is defined as:

CJS(P || Q) = \sum P(x) \log \frac{P(x)}{0.5 (P(x) + Q(x))} + \frac{1}{2 \ln 2} \sum (Q(x) - P(x))

The symmetric metric is defined as:

CJS_{dist}(P || Q) = \sqrt{CJS(P || Q) + CJS(Q || P)}

This has an upper bound of \sqrt{ \sum (P(x) + Q(x))}

Value

distance value based on CJS computation.

References

Nguyen H-V., Vreeken J. (2015). Non-parametric Jensen-Shannon Divergence. In: Appice A., Rodrigues P., Santos Costa V., Gama J., Jorge A., Soares C. (eds) Machine Learning and Knowledge Discovery in Databases. ECML PKDD 2015. Lecture Notes in Computer Science, vol 9285. Springer, Cham. doi:10.1007/978-3-319-23525-7_11

Examples

x <- rnorm(100)
y <- rnorm(100, 2, 2)
cjs_dist(x, y, x_weights = NULL, y_weights = NULL)

Create data structure for priorsense

Description

Create a data structure that contains all required data and functions for priorsense

Usage

create_priorsense_data(x, ...)

## Default S3 method:
create_priorsense_data(
  x,
  fit = NULL,
  log_prior_fn = log_prior_draws,
  log_lik_fn = log_lik_draws,
  log_prior = NULL,
  log_lik = NULL,
  log_ratio_fn = NULL,
  log_prior_name = "lprior",
  log_lik_name = "log_lik",
  ...
)

## S3 method for class 'stanfit'
create_priorsense_data(x, ...)

## S3 method for class 'CmdStanFit'
create_priorsense_data(x, ...)

## S3 method for class 'draws'
create_priorsense_data(x, ...)

## S3 method for class 'rjags'
create_priorsense_data(x, ...)

Arguments

Value

A priorsense_data object, which contains the data and functions to run sensitivity analyses.

Examples

x <- example_powerscale_model()
drw <- x$draws

psd <- create_priorsense_data(drw)

Example Stan model for power-scaling

Description

Provides example models (with data) that are ready for use with power-scaling.

Usage

example_powerscale_model(model = "univariate_normal")

Arguments

Value

List containing model code and corresponding data.

Examples

ex_normal <- example_powerscale_model(model = "univariate_normal")

ex_eightschools <- example_powerscale_model(model = "eight_schools")

Extract log likelihood draws

Description

Extract log likelihood from fitted model and return as a draws object.

Usage

log_lik_draws(x, ...)

## S3 method for class 'stanfit'
log_lik_draws(x, joint = FALSE, log_lik_name = "log_lik", ...)

## S3 method for class 'CmdStanFit'
log_lik_draws(x, joint = FALSE, log_lik_name = "log_lik", ...)

## S3 method for class 'draws'
log_lik_draws(x, joint = FALSE, log_lik_name = "log_lik", ...)

Arguments

Value

A draws_array object containing log_lik values.

Examples

ex <- example_powerscale_model()
drw <- ex$draws

log_lik_draws(drw)

Extract log prior draws

Description

Extract log likelihood from fitted model and return as a draws object.

Usage

log_prior_draws(x, ...)

## S3 method for class 'stanfit'
log_prior_draws(x, joint = FALSE, log_prior_name = "lprior", ...)

## S3 method for class 'CmdStanFit'
log_prior_draws(x, joint = FALSE, log_prior_name = "lprior", ...)

## S3 method for class 'draws'
log_prior_draws(x, joint = FALSE, log_prior_name = "lprior", ...)

Arguments

Value

A draws_array object containing log_prior values.

Examples

ex <- example_powerscale_model()
drw <- ex$draws

log_prior_draws(drw)

Power-scale gradients

Description

Calculate the numerical derivative of posterior quantities/divergence with respect to power-scaling the specified component (prior or likelihood). This is done using importance sampling (and optionally moment matching).

Usage

powerscale_gradients(x, ...)

## Default S3 method:
powerscale_gradients(
  x,
  log_prior_name = "lprior",
  log_lik_name = "log_lik",
  ...
)

## S3 method for class 'priorsense_data'
powerscale_gradients(
  x,
  variable = NULL,
  component = c("prior", "likelihood"),
  type = c("quantities", "divergence"),
  lower_alpha = 0.99,
  upper_alpha = 1.01,
  div_measure = "cjs_dist",
  measure_args = list(),
  moment_match = FALSE,
  k_threshold = 0.5,
  resample = FALSE,
  transform = NULL,
  prediction = NULL,
  scale = FALSE,
  prior_selection = NULL,
  likelihood_selection = NULL,
  ...
)

Arguments

Value

Maximum of the absolute derivatives above and below alpha = 1.

Examples

ex <- example_powerscale_model()
drw <- ex$draws

powerscale_gradients(drw)

Prior/likelihood power-scaling perturbation

Description

Estimate posterior draws based on power-scaling perturbations of prior or likelihood using importance sampling (and optionally moment matching).

Usage

powerscale(x, ...)

## Default S3 method:
powerscale(
  x,
  component,
  alpha,
  moment_match = FALSE,
  k_threshold = NULL,
  resample = FALSE,
  transform = NULL,
  prediction = NULL,
  variable = NULL,
  selection = NULL,
  log_prior_name = "lprior",
  log_lik_name = "log_lik",
  ...
)

## S3 method for class 'priorsense_data'
powerscale(
  x,
  component,
  alpha,
  moment_match = FALSE,
  k_threshold = NULL,
  resample = FALSE,
  transform = NULL,
  prediction = NULL,
  variable = NULL,
  selection = NULL,
  log_prior_name = "lprior",
  log_lik_name = "log_lik",
  ...
)

powerscale_sequence(x, ...)

## Default S3 method:
powerscale_sequence(
  x,
  lower_alpha = 0.8,
  upper_alpha = 1/lower_alpha,
  length = 3,
  variable = NULL,
  component = c("prior", "likelihood"),
  moment_match = FALSE,
  k_threshold = 0.5,
  resample = FALSE,
  transform = NULL,
  prediction = NULL,
  auto_alpha_range = FALSE,
  symmetric = TRUE,
  prior_selection = NULL,
  likelihood_selection = NULL,
  ...
)

## S3 method for class 'priorsense_data'
powerscale_sequence(
  x,
  lower_alpha = 0.8,
  upper_alpha = 1/lower_alpha,
  length = 3,
  variable = NULL,
  component = c("prior", "likelihood"),
  moment_match = FALSE,
  k_threshold = NULL,
  resample = FALSE,
  transform = NULL,
  prediction = NULL,
  auto_alpha_range = FALSE,
  symmetric = TRUE,
  prior_selection = NULL,
  likelihood_selection = NULL,
  ...
)

Arguments

Value

A powerscaled_draws or powerscaled_sequence object, which contains the estimated posterior draws resulting from the power-scaling perturbations and details of the perturbation and estimation methods.

References

Kallioinen, N., Paananen, T., Bürkner, P-C., Vehtari, A. (2023). Detecting and diagnosing prior and likelihood sensitivity with power-scaling perturbations. Statistics and Computing. 34(57). doi:10.1007/s11222-023-10366-5

Vehtari, A., Simpson, D., Gelman, A., Yao, Y., and Gabry, J. (2024). Pareto smoothed importance sampling. Journal of Machine Learning Research. 25(72). https://jmlr.org/papers/v25/19-556.html

Paananen, T., Piironen, J., Bürkner, P-C., Vehtari, A. (2021). Implicitly adaptive importance sampling. Statistics and Computing. 31(16). doi:10.1007/s11222-020-09982-2

Examples

ex <- example_powerscale_model()

powerscale(ex$draws, component = "prior", alpha = 0.5)

powerscale_sequence(ex$draws)

Power-scaling sensitivity analysis

Description

Calculates the prior/likelihood sensitivity based on power-scaling perturbations. This is done using importance sampling (and optionally moment matching).

Usage

powerscale_sensitivity(x, ...)

## Default S3 method:
powerscale_sensitivity(
  x,
  variable = NULL,
  lower_alpha = 0.99,
  upper_alpha = 1.01,
  div_measure = "cjs_dist",
  measure_args = list(),
  component = c("prior", "likelihood"),
  sensitivity_threshold = 0.05,
  moment_match = FALSE,
  k_threshold = 0.5,
  resample = FALSE,
  transform = NULL,
  prediction = NULL,
  prior_selection = NULL,
  likelihood_selection = NULL,
  log_prior_name = "lprior",
  log_lik_name = "log_lik",
  num_args = NULL,
  ...
)

## S3 method for class 'priorsense_data'
powerscale_sensitivity(
  x,
  variable = NULL,
  lower_alpha = 0.99,
  upper_alpha = 1.01,
  div_measure = "cjs_dist",
  measure_args = list(),
  component = c("prior", "likelihood"),
  sensitivity_threshold = 0.05,
  moment_match = FALSE,
  k_threshold = 0.5,
  resample = FALSE,
  transform = NULL,
  prediction = NULL,
  prior_selection = NULL,
  likelihood_selection = NULL,
  num_args = NULL,
  ...
)

## S3 method for class 'CmdStanFit'
powerscale_sensitivity(x, ...)

## S3 method for class 'stanfit'
powerscale_sensitivity(x, ...)

Arguments

Value

Table of sensitivity values for each specified variable.

References

Kallioinen, N., Paananen, T., Bürkner, P-C., Vehtari, A. (2023). Detecting and diagnosing prior and likelihood sensitivity with power-scaling perturbations. Statistics and Computing. 34(57). doi:10.1007/s11222-023-10366-5

Vehtari, A., Simpson, D., Gelman, A., Yao, Y., and Gabry, J. (2024). Pareto smoothed importance sampling. Journal of Machine Learning Research. 25(72). https://jmlr.org/papers/v25/19-556.html

Paananen, T., Piironen, J., Bürkner, P-C., Vehtari, A. (2021). Implicitly adaptive importance sampling. Statistics and Computing. 31(16). doi:10.1007/s11222-020-09982-2

Examples

ex <- example_powerscale_model()
powerscale_sensitivity(ex$draws)

Derivative with respect to power-scaling

Description

Calculate the analytical derivative of a quantity with respect to power-scaling prior or likelihood.

Usage

powerscale_derivative(x, log_component, quantity = "mean", ...)

Arguments

Value

Derivative of the quantity with respect to log2 of the power-scaling factor (alpha).

Examples

example_model <- example_powerscale_model()
draws <- example_model$draws
log_prior <- log_prior_draws(draws, joint = TRUE)
posterior::summarise_draws(
    posterior::subset_draws(draws, variable = c("mu", "sigma")),
    mean,
    mean_sens = ~powerscale_derivative(.x, log_prior, quantity = "mean")
)

Diagnostic plots for power-scaling sensitivity

Description

Various diagnostic plots for power-scaling sensitivity. See Plot Descriptions below for details.

Usage

powerscale_plot_dens(x, ...)

powerscale_plot_ecdf(x, ...)

## S3 method for class 'powerscaled_sequence'
powerscale_plot_ecdf(
  x,
  variable = NULL,
  resample = FALSE,
  length = 3,
  facet_rows = "component",
  help_text = getOption("priorsense.plot_help_text", TRUE),
  colors = NULL,
  variables_per_page = getOption("priorsense.plot_variables_per_page", 6),
  ...
)

powerscale_plot_quantities(x, ...)

## S3 method for class 'powerscaled_sequence'
powerscale_plot_quantities(
  x,
  variable = NULL,
  quantity = c("mean", "sd"),
  div_measure = "cjs_dist",
  resample = FALSE,
  measure_args = NULL,
  mcse = TRUE,
  quantity_args = NULL,
  help_text = getOption("priorsense.plot_help_text", TRUE),
  colors = NULL,
  variables_per_page = getOption("priorsense.plot_variables_per_page", 6),
  ...
)

Arguments

Value

A ggplot object (or a priorsense_plot object which is a list of ggplot objects if there is more than one page) that can be further customized using the ggplot2 package.

Plot Descriptions

powerscale_plot_dens()

Kernel density plot of power-scaled posterior draws with respect to power-scaling.

powerscale_plot_ecdf()

Empirical cumulative distribution function plot of power-scaled posterior draws with respect to power-scaling.

powerscale_plot_quantities()

Plot of posterior quantities with respect to power-scaling.

Examples

ex <- example_powerscale_model()

powerscale_plot_dens(ex$draws)

brms predictions as draws

Description

Create predictions using brms functions and convert them into draws format

Usage

predictions_as_draws(
  x,
  predict_fn,
  prediction_names = NULL,
  warn_dims = getOption("priorsense.warn", TRUE),
  ...
)

Arguments

Value

draws array of predictions

Examples

## Not run: 
library(brms)

if ("log_prior_draws.brmsfit" %in% methods(log_prior_draws) &&
    ("log_lik_draws.brmsfit" %in% methods(log_lik_draws))) {
  fit <- brm(
    yield ~ N * P * K,
    data = npk,
    prior = prior(normal(0, 1), class = "b"),
    refresh = 0
  )

  powerscale_sensitivity(
      fit,
      variable = "_pred",
      prediction = function(x) predictions_as_draws(
                                 x, brms::posterior_epred
                               )
  )
}

## End(Not run)