Package {bayesics}

Title:	Bayesian Analyses for One- and Two-Sample Inference and Regression Methods
Version:	3.0.0
Maintainer:	Daniel K. Sewell <daniel-sewell@uiowa.edu>
Description:	Perform fundamental analyses using Bayesian parametric and non-parametric inference (regression, anova, 1 and 2 sample inference, non-parametric tests, etc.). (Practically) no Markov chain Monte Carlo (MCMC) is used; all exact finite sample inference is completed via closed form solutions or else through posterior sampling automated to ensure precision in interval estimate bounds. Diagnostic plots for model assessment, and key inferential quantities (point and interval estimates, probability of direction, region of practical equivalence, and Bayes factors) and model visualizations are provided. Bayes factors are computed either by the Savage Dickey ratio given in Dickey (1971) <doi:10.1214/aoms/1177693507> or by Chib's method as given in <doi:10.1080/01621459.1995.10476635>. Interpretations are from Kass and Raftery (1995) <doi:10.1080/01621459.1995.10476572>. ROPE bounds are based on discussions in Kruschke (2018) <doi:10.1177/2515245918771304>. Methods for determining the number of posterior samples required are described in Doss et al. (2014) <doi:10.1214/14-EJS957>. Bayesian model averaging is done in part by Feldkircher and Zeugner (2015) <doi:10.18637/jss.v068.i04>. Methods for contingency table analysis is described in Gunel et al. (1974) <doi:10.1093/biomet/61.3.545>. Variational Bayes (VB) methods are described in Salimans and Knowles (2013) <doi:10.1214/13-BA858>. Mediation analysis uses the framework described in Imai et al. (2010) <doi:10.1037/a0020761>. The loss-likelihood bootstrap used in the non-parametric regression modeling is described in Lyddon et al. (2019) <doi:10.1093/biomet/asz006>. Non-parametric survival methods are described in Qing et al. (2023) <doi:10.1002/pst.2256>. Methods used for the Bayesian Wilcoxon signed-rank analysis is given in Chechile (2018) <doi:10.1080/03610926.2017.1388402> and for the Bayesian Wilcoxon rank sum analysis in Chechile (2020) <doi:10.1080/03610926.2018.1549247>. Correlation analysis methods are carried out by Barch and Chechile (2023) <doi:10.32614/CRAN.package.DFBA>, and described in Lindley and Phillips (1976) <doi:10.1080/00031305.1976.10479154> and Chechile and Barch (2021) <doi:10.1016/j.jmp.2021.102638>. See also Chechile (2020, ISBN: 9780262044585).
License:	GPL (≥ 3)
Encoding:	UTF-8
Depends:	R (≥ 4.1.0)
Suggests:	datasets, rstanarm, knitr, splines, testthat (≥ 3.0.0)
Imports:	tidyr, dplyr, rlang, janitor, extraDistr, mvtnorm, Matrix, future, future.apply, ggplot2, patchwork, BMS, cluster, DFBA, tibble, survival, stringr
Config/testthat/edition:	3
URL:	https://github.com/dksewell/bayesics
BugReports:	https://github.com/dksewell/bayesics/issues
Config/roxygen2/version:	8.0.0
NeedsCompilation:	no
Packaged:	2026-07-13 20:01:10 UTC; dksewell
Author:	Daniel K. Sewell [aut, cre, cph], Alan Arakkal [aut]
Repository:	CRAN
Date/Publication:	2026-07-13 20:30:02 UTC

bayesics: Bayesian Analyses for One- and Two-Sample Inference and Regression Methods

Description

The bayesics package meant to act as a Bayesian analog to many of the procedures implemented in the stats package. It includes methods for one- and two-sample inference including 2-way contingency table analyses, parametric and non-parametric regression, Bayesian model averaging, and mediation analysis.

Details

The design of bayesics emphasizes inference and model assessment, rather than algorithmic tuning or sampling diagnostics. Most modeling functions provide familiar generic functions such as print(), summary(), plot(), and predict(), while introducing Bayesian analogues to classical generics. In particular, bayesics defines new generics for regression-based inference, including credint() (replacing confint) and get_posterior_draws().

The central aim of the package is principled Bayesian inference and interpretation. Standard reported quantities include posterior point and interval estimates, probabilities that estimands fall within a region of practical equivalence (ROPE), and probabilities of direction. Interpretations always follow Bayes factors to ensure clarity.

Estimation methods relying on MCMC require knowledge of not just chain convergence, but also assessing whether sufficient accuracy is obtained for the point estimates and critically the credible interval bounds. bayesics avoids the use of MCMC and instead relies on either closed-form solutions or independent posterior draws. For example, Bayes factors are computed analytically using either the Savage-Dickey ratio or via Chib's method. (Additionally, fractional Bayes factors are also implemented- again, analytically- for linear models.) when posterior draws are required to perform inference, the number of posterior draws is automatically selected in order to ensure sufficiently accurate results (through the user-specified argument mc_error). This entirely eliminates any need for users to perform algorithmic assessments.

Model assessments, however, are always critical, and towards this bayesics provides bayes_pvalue(), a function designed to assess both linear and generalized linear models through Bayesian p-values. While the deviance is the default, any test statistic can be incorporated.

When model diagnostics fail, non-parametric methods may be utilized instead. bayesics provides the function np_glm_b which implements the loss-likelihood bootstrap, a general Bayes inferential method (Lyddon et al., 2019).

Besides parametric and non-parametric regression techniques, non-regression methods are also implemented, such as tests of correlation or comparing two samples of count data, giving Bayesian equivalents to well-used functions such as cor_test_b (replacing cor.test), t_test_b (replacing t.test), and many others.

Author(s)

Daniel K. Sewell

References

Barch DH, Chechile RA (2023). DFBA: Distribution-Free Bayesian Analysis. doi:10.32614/CRAN.package.DFBA

Chechile, R. A. (2018) A Bayesian analysis for the Wilcoxon signed-rank statistic. Communications in Statistics - Theory and Methods, https://doi.org/10.1080/03610926.2017.1388402

Chechile, R.A. (2020). Bayesian Statistics for Experimental Scientists: A General Introduction Using Distribution_Free Statistics. Cambridge: MIT Press.

Chechile, R.A. (2020). A Bayesian analysis for the Mann-Whitney statistic. Communications in Statistics – Theory and Methods 49(3): 670-696. https://doi.org/10.1080/03610926.2018.1549247.

Chechile, R.A., & Barch, D.H. (2021). A distribution-free, Bayesian goodness-of-fit method for assessing similar scientific prediction equations. Journal of Mathematical Psychology. https://doi.org/10.1016/j.jmp.2021.102638

Chib, S. (1995). Marginal Likelihood from the Gibbs Output. Journal of the American Statistical Association, 90(432), 1313–1321. https://doi.org/10.1080/01621459.1995.10476635

James M. Dickey. "The Weighted Likelihood Ratio, Linear Hypotheses on Normal Location Parameters." Ann. Math. Statist. 42 (1) 204 - 223, February, 1971. https://doi.org/10.1214/aoms/1177693507

Charles R. Doss, James M. Flegal, Galin L. Jones, Ronald C. Neath "Markov chain Monte Carlo estimation of quantiles," Electronic Journal of Statistics, Electron. J. Statist. 8(2), 2448-2478, (2014)

Feldkircher, M. and S. Zeugner (2015): Bayesian Model Averaging Employing Fixed and Flexible Priors: The BMS Package for R, Journal of Statistical Software 68(4).

Gunel, Erdogan & Dickey, James (1974). Bayes factors for independence in contingency tables, Biometrika, 61(3), Pages 545–557, https://doi.org/10.1093/biomet/61.3.545

Imai, Kosuke, et al. “A General Approach to Causal Mediation Analysis.” Psychological Methods, vol. 15, no. 4, 2010, pp. 309–34, https://doi.org/10.1037/a0020761.

Kass, R. E., & Raftery, A. E. (1995). Bayes Factors. Journal of the American Statistical Association, 90(430), 773–795.

Kruschke JK. Rejecting or Accepting Parameter Values in Bayesian Estimation. Advances in Methods and Practices in Psychological Science. 2018;1(2):270-280. doi:10.1177/2515245918771304

Lindley, D. V., & Phillips, L. D. (1976). Inference for a Bernoulli process (a Bayesian view). The American Statistician, 30, 112-119.

S P Lyddon, C C Holmes, S G Walker, General Bayesian updating and the loss-likelihood bootstrap, Biometrika, Volume 106, Issue 2, June 2019, Pages 465–478, https://doi.org/10.1093/biomet/asz006

O’Hagan, Anthony. “Fractional Bayes Factors for Model Comparison.” Journal of the Royal Statistical Society. Series B (Methodological), vol. 57, no. 1, 1995, pp. 99–138. https://doi.org/10.1111/j.2517-6161.1995.tb02017.x

Qing Y, Thall PF, Yuan Y. A Bayesian piecewise exponential phase II design for monitoring a time-to-event endpoint. Pharm Stat. 2023 Jan;22(1):34-44. doi: 10.1002/pst.2256. Epub 2022

Tim Salimans. David A. Knowles. "Fixed-Form Variational Posterior Approximation through Stochastic Linear Regression." Bayesian Anal. 8 (4) 837 - 882, December 2013. https://doi.org/10.1214/13-BA858

See Also

Useful links:

• https://github.com/dksewell/bayesics

• Report bugs at https://github.com/dksewell/bayesics/issues

Compute AIC, BIC, DIC, or WAIC for aov_b or lm_b Objects. (Lower is Better.)

Description

Compute AIC, BIC, DIC, or WAIC for aov_b or lm_b Objects. (Lower is Better.)

Usage

DIC(object, ...)

## S3 method for class 'lm_b'
BIC(object, ...)

## S3 method for class 'lm_b'
AIC(object, ...)

## S3 method for class 'lm_b'
DIC(object, seed = 1, mc_error = 0.5, ...)

## S3 method for class 'aov_b'
DIC(object, ...)

## S3 method for class 'lm_b'
WAIC(object, seed = 1, mc_error = 0.5, ...)

## S3 method for class 'aov_b'
WAIC(object, ...)

Arguments

Details

AIC and BIC are constructed using the posterior mean. DIC and WAIC are computed via independent posterior sampling, ensuring that the final computed numbers is within mc_error of the actual DIC/WAIC with high probability.

Value

Numeric (or in the case of DIC, a numeric vector)

Examples

set.seed(2025)
N = 500
test_data <-
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome <-
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )
fit1 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data)
AIC(fit1)
BIC(fit1)
DIC(fit1)
WAIC(fit1)

Create a Survival Object

Description

Create a survival object, usually used as a response variable in a model formula. Argument matching is special for this function, see Details under Surv. This is a restricted wrapper around Surv and currently supports only right-censored data.

Usage

Surv(...)

Arguments

Value

An object of class "Surv".

References

Therneau T (2024). A Package for Survival Analysis in R. R package version 3.8-3, https://CRAN.R-project.org/package=survival.

Examples

set.seed(2025)
N = 300
test_data = 
  data.frame(outcome = 
               rweibull(N,2,5))
test_data$observed = 
  ifelse(test_data$outcome >= 7, 0, 1)
test_data$outcome =
  ifelse(dplyr::near(test_data$observed,1), test_data$outcome, 7)
Surv(test_data$outcome,
     test_data$observed)

Analysis of Variance Using Bayesian Methods

Description

Analysis of Variance Using Bayesian Methods

Usage

aov_b(
  formula,
  data,
  heteroscedastic = TRUE,
  prior_mean_mu,
  prior_mean_nu = 0.001,
  prior_var_shape = 0.001,
  prior_var_rate = 0.001,
  CI_level = 0.95,
  ROPE = 0.1,
  contrasts,
  improper = FALSE,
  seed = 1,
  mc_error = 0.002,
  compute_bayes_factor = TRUE
)

Arguments

Details

MODEL: The likelihood model is given by

y_{gi} \overset{iid}{\sim} N(\mu_g,\sigma^2_g),

(although if heteroscedastic is set to FALSE, \sigma^2_g=\sigma^2_h \forall g,h).

The prior is given by

\mu_g|\sigma^2_g \overset{iid}{\sim} N\left(\mu,\frac{\sigma^2_g}{\nu}\right), \\ \sigma^2_g \overset{iid}{\sim} \Gamma^{-1}(a/2,b/2),

where mu is set by prior_mean_mu, nu is set by prior_mean_nu, a is set by prior_var_shape, and b is set by prior_var_rate.

The posterior is

\mu_g|y,\sigma^2_g \overset{iid}{\sim} N\left(\hat\mu_g,\frac{\sigma^2_g}{\nu_g}\right), \\ \sigma^2_g|y \overset{iid}{\sim} \Gamma^{-1}(a_g/2,b_g/2),

where \hat\mu_g, \nu_g, a_g, and b_g are all returned by aov_b in the named element posterior_parameters.

ROPE:

If missing, the ROPE bounds will be given under the principle of "half of a small effect size." Using Cohen's D of 0.2 as a small effect size, the ROPE is defined in terms of -0.1 < Cohen's D < 0.1.

Value

Object of class aov_b and lm_b.

References

Charles R. Doss, James M. Flegal, Galin L. Jones, Ronald C. Neath "Markov chain Monte Carlo estimation of quantiles," Electronic Journal of Statistics, Electron. J. Statist. 8(2), 2448-2478, (2014)

Examples

# Create data
set.seed(2025)
N = 500
test_data = 
  data.frame(x1 = rep(letters[1:5],N/5))
test_data$outcome = 
  rnorm(N,-1 + 2 * (test_data$x1 %in% c("d","e")) )

# Fit 1-way ANOVA model
fit1 <-
  aov_b(outcome ~ x1,
        test_data,
        prior_mean_mu = 2,
        prior_mean_nu = 0.5,
        prior_var_shape = 0.01,
        prior_var_rate = 0.01)
fit1
summary(fit1)
plot(fit1)
coef(fit1)
credint(fit1)
credint(fit1,
        CI_level = 0.99)
vcov(fit1)
fit1_predictions <- 
  predict(fit1,
          CI_level = 0.99,
          PI_level = 0.9)
AIC(fit1)
BIC(fit1)
DIC(fit1)
WAIC(fit1)

# Implement contrasts
## One contrast
fit2 <-
  aov_b(outcome ~ x1,
        test_data,
        mc_error = 0.01,
        contrasts = c(-1/3,-1/3,-1/3,1/2,1/2))
fit2$contrasts
summary(fit2)
## Multiple contrasts
fit3 <-
  aov_b(outcome ~ x1,
        test_data,
        mc_error = 0.01,
        contrasts = rbind(c(-1/3,-1/3,-1/3,1/2,1/2),
                          c(-1/3,-1/3,-1/3,1,0)))
fit3$contrasts
summary(fit3)

b_procedure Objects

Description

Objects of class b_procedure represent the result of a Bayesian procedure, including the data, prior specification, posterior summaries, and optional plotting output.

Details

A b_procedure object is a named list with the following components (not all bayesian procedures will yield objects with all of these entries):

name

Character string giving the name of the procedure.

data

A tibble containing the data used in the analysis.

print_data

Logical; whether the data should be printed by print.b_procedure.

CI_level

Numeric scalar giving the credible interval level as provided by the user.

sampling_design

Character; Used by independence_b.

prior

Character string describing the prior used.

results

A tibble containing posterior summaries with columns:

• Quantity: character

• Post Mean: numeric

• Lower: numeric

• Upper: numeric

• ROPE: optional numeric

• ROPE_lower_bound,ROPE_upper_bound: optional numeric

PDir

List containing:

• description: Character describing the probability of direction.

• pdir: Numeric scalar giving the probability of direction.

BF

List containing:

• description: character

• BF: numeric scalar giving the Bayes factor

• interpretation: character

overall_ROPE

List containing:

• description: character

• Pr_in_ROPE: numeric

plot

A ggplot object associated with the procedure.

object_fit

If applicable, the underlying fitted model object (e.g., aov_b).

notes

character vector

S3 methods

The following methods are available for lm_b class objects: print() and plot().

See Also

print.b_procedure, plot.b_procedure

Examples

cc_fit <- case_control_b(matrix(c(8,47,1,26),2,2))
cc_fit
plot(cc_fit)

Bayes Factors for lm_b, glm_b, and survfit_b Objects

Description

Bayes factors for Bayesian regression objects using the Savage-Dickey ratio

Usage

bayes_factors(object, ...)

## S3 method for class 'lm_b'
bayes_factors(object, by = "coefficient", ...)

## S3 method for class 'glm_b'
bayes_factors(object, by = "coefficient", ...)

## S3 method for class 'survfit_b'
bayes_factors(object, object2, ...)

Arguments

Details

Bayes factors are given in terms of favoring the two-tailed alternative hypothesis vs. the null hypothesis that the regression coefficient equals zero. Currently implemented for lm_b or glm_b objects. Note that for glm_b objects, if importance sampling was used, the model will be refit using fixed form variational Bayes to get the multivariate posterior density. The Bayes factor is then computed using the Savage-Dickey ratio.

Interpretation is taken from Kass and Raftery.

Value

A tibble with Bayes factors and interpretations.

References

Chib, S. (1995). Marginal Likelihood from the Gibbs Output. Journal of the American Statistical Association, 90(432), 1313–1321. https://doi.org/10.1080/01621459.1995.10476635

James M. Dickey. "The Weighted Likelihood Ratio, Linear Hypotheses on Normal Location Parameters." Ann. Math. Statist. 42 (1) 204 - 223, February, 1971. https://doi.org/10.1214/aoms/1177693507

Kass, R. E., & Raftery, A. E. (1995). Bayes Factors. Journal of the American Statistical Association, 90(430), 773–795.

Examples

# Generate some binomial data
set.seed(2025)
N = 500
test_data = 
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome = 
  rbinom(N,1,1.0 / (1.0 + exp(-(-2 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) ))))


# Fit a GLM
fit <-
  glm_b(outcome ~ x1 + x2 + x3,
        data = test_data,
        family = binomial(),
        seed = 2025)

# Compute the BF for each coefficient
bayes_factors(fit)

# Compute the BF for each variable
bayes_factors(fit,
              by = "variable")

Bayesian P-values for Regression Models

Description

Bayesian P-values for Regression Models

Usage

bayes_pvalue(object, statistic, mc_error, seed, ...)

## S3 method for class 'lm_b'
bayes_pvalue(object, statistic, mc_error = 0.005, seed = 1, ...)

## S3 method for class 'aov_b'
bayes_pvalue(object, statistic, mc_error = 0.005, seed = 1, ...)

Arguments

Details

Overview:

Bayesian p-values are measures of how well the model match the data, as evaluated through the predictive posterior distribution. While they can be extremely flexible - testing very specific aspects of the model being fitted-, the default setting for bayes_pvalue is the deviance, acting to do an overall goodness-of-fit test. More generally, a Bayesian p-value takes a test statistic of the data and the model parameters T(y,\theta) and compares the posterior probability that the test statistic evaluated at the observed data compared to data randomly generated according to the model. I.e., the Bayesian p-value is given by

\Pr(T(y_{obs},\theta) > T(y_{pred},\theta) | y_{pred}).

MC error:

The number of posterior samples is determined by the fact that if the true Bayesian p-value is p, the Monte Carlo estimate of the Bayesian p-value will have 99% probability of being with \approx 2.3 \sqrt{p(1-p)/L}, where L is the number of posterior draws. The worst case scenario, in terms of MC variance is when the Bayesian p-value is p=0.5. However, in such a case, there will be no question that the model fit is adequate and we can afford a much larger MC error. It is nearer thresholds (typically Bayesian p-values less than 0.05 or greater than 0.95 are cause for alarm) where we need more precise estimates of what the Bayesian p-value actually is, but at near these boundaries the MC error is much smaller than the worst case scenario of p=0.5. Hence we compute the number of posterior samples to be within the user-specified mc_error at p(1-p)=(0.15)(0.85).

The Monte Carlo error is implemented in such a way so as to obtain as few posterior samples as necessary to ensure a small probability of incorrectly determining an inadequate fit when the fit is good and vice versa.

Value

named list with:

• bpvalue - numeric between 0 and 1, giving the Bayesian p-value

• statistic_posterior_draws - tibble with two columns, one for T(y_{obs},\theta) and the other for T(y_{pred},\theta). See details.

Examples

# Create some data
set.seed(2026)
N = 500
test_data = 
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome = 
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )

# Fit a linear regression model
fit = 
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data)

# Compute the Bayesian p-value based on the deviance
bp_deviance = 
  bayes_pvalue(fit)
bp_deviance$bpvalue # We want a value near 0.5, say in (0.05,0.95).
plot(T_y_observed ~ T_y_predicted,
     data = bp_deviance$statistic_posterior_draws,
     xlab = expression(T(y[pred],theta)),
     ylab = expression(T(y[obs],theta)),
     pch = 16,
     cex = 0.1,
     col = gray(0.5,0.25))
abline(0,1,
       lwd = 2)

# Use a custom test statistic
bp_sw = 
  bayes_pvalue(fit,
               statistic = 
                 function(y,mu,dispersion){
                   shapiro.test((y - mu)/sqrt(dispersion))$statistic
                 }
  )
bp_sw$bpvalue
 

Bayesian Model Averaging

Description

Estimates and CIs from BMA

Usage

bma_inference(
  formula,
  data,
  zellner_g = nrow(data),
  CI_level = 0.95,
  ROPE,
  mcmc_draws = 10000,
  n_models = 500,
  mc_error = 0.001,
  seed = 1,
  compute_residuals = TRUE,
  ...
)

Arguments

Details

bma_inference leverages the bms function from its eponymous R package, and then uses lm_b to obtain inference on the regression coefficients for Bayesian model averaging.

Value

Object of class lm_b_bma and lm_b.

References

Feldkircher, M. and S. Zeugner (2015): Bayesian Model Averaging Employing Fixed and Flexible Priors: The BMS Package for R, Journal of Statistical Software 68(4).

Examples

# Create data
set.seed(2025)
N = 500
test_data = 
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5],
             x4 = rnorm(N),
             x5 = rnorm(N),
             x6 = rnorm(N),
             x7 = rnorm(N),
             x8 = rnorm(N),
             x9 = rnorm(N),
             x10 = rnorm(N))
test_data$outcome = 
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )

# Fit linear model using Bayesian model averaging
fit <-
  bma_inference(outcome ~ .,
                test_data,
                user.int = FALSE)
summary(fit)
coef(fit)
credint(fit)
plot(fit)

Case-Control Analysis

Description

Bayesian analysis of a case-control study (without covariates).

Usage

case_control_b(
  cases,
  controls,
  x,
  large_sample_approx,
  ROPE,
  prior_mean = 0,
  prior_sd = log(10)/1.96,
  plot = TRUE,
  CI_level = 0.95,
  seed = 1,
  mc_error = 0.005
)

Arguments

Details

If large_sample_approx = TRUE (the default if left missing and all cell counts are at least 5), then the likelihood is

\log(\hat\omega) \sim N\left(\log(\omega),\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} \right),

where \omega is the odds ratio, \hat\omega is the empirical odds ratio, n_{ij}, i,j = 1,2 are the cells of the 2x2 contingency table. The prior on \log\omega is

\log\omega \sim N(\texttt{prior\_mean},\texttt{prior\_sd}^2).

If the large sample approximation is not used, then inference is made on the odds ratio by instead putting uniform priors on \Pr(exposure|outcome).

Value

An object of class b_procedure-class.

Examples

case_control_b(matrix(c(8,47,1,26),2,2))

case_control_b(c(8,47),
               c(1,26))

Test of Independence for 2-way Contingency Tables

Description

Test of Independence for 2-way Contingency Tables

Usage

independence_b(
  x,
  sampling_design = c("multinomial", "fixed rows", "fixed columns"),
  ROPE,
  prior = c("jeffreys", "uniform"),
  prior_shapes,
  CI_level = 0.95,
  seed = 1,
  mc_error = 0.002
)

Arguments

Details

For a 2-way contingency table with R rows and C columns, evaluate the probability that

• the joint probabilities p_{ij} are all within the ROPE of p_{i\cdot}\times p_{\cdot j} for sampling_design = "multinomial"

• the probabilities p_{j|i} are all within the ROPE of p_{\cdot j} if sampling_design = "fixed rows" or "fixed columns"

Value

An object of class b_procedure-class.

References

Gunel, Erdogan & Dickey, James (1974). Bayes factors for independence in contingency tables, Biometrika, 61(3), Pages 545–557, https://doi.org/10.1093/biomet/61.3.545

Kass, R. E., & Raftery, A. E. (1995). Bayes Factors. Journal of the American Statistical Association, 90(430), 773–795.

Examples

# Generate data
set.seed(2025)
N = 500
nR = 5
nC = 3
dep_probs = 
  extraDistr::rdirichlet(1,rep(2,nR*nC)) |> 
  matrix(nR,nC)

# Multinomial sampling
## Test independence
independence_b(round(N * dep_probs))

## Use other priors
independence_b(round(N * dep_probs),
               prior = "uniform")
independence_b(round(N * dep_probs),
               prior_shapes = 2)
independence_b(round(N * dep_probs),
               prior_shapes = matrix(1:(nR*nC),nR,nC))

# Fixed marginals
independence_b(round(N * dep_probs),
               sampling_design = "fixed rows")
independence_b(round(N * dep_probs),
               sampling_design = "fixed col")

Coefficient Extraction for bayesics Objects

Description

Coefficient Extraction for bayesics Objects

Usage

## S3 method for class 'lm_b'
coef(object, ...)

## S3 method for class 'aov_b'
coef(object, ...)

Arguments

Value

vector of coefficients

Examples

set.seed(2025)
N = 500
test_data = 
  data.frame(x1 = rep(letters[1:5],N/5))
test_data$outcome = 
  rnorm(N,-1 + 2 * (test_data$x1 %in% c("d","e")) )

# Fit 1-way ANOVA model
fit1 <-
  aov_b(outcome ~ x1,
        test_data,
        prior_mean_mu = 2,
        prior_mean_nu = 0.5,
        prior_var_shape = 0.01,
        prior_var_rate = 0.01)
coef(fit1)

Test for Association/Correlation Between Paired Samples via Kendall's tau

Description

Test for Association/Correlation Between Paired Samples via Kendall's tau

Usage

cor_test_b(x, ...)

## Default S3 method:
cor_test_b(
  x,
  y,
  tau = 0,
  ROPE,
  prior = c("centered", "uniform", "positive", "negative"),
  prior_shapes,
  CI_level = 0.95,
  plot = TRUE,
  ...
)

## S3 method for class 'formula'
cor_test_b(
  formula,
  data,
  tau = 0,
  ROPE,
  prior = "centered",
  prior_shapes,
  CI_level = 0.95,
  plot = TRUE,
  ...
)

Arguments

Details

cor_test_b relies on the robust Kendall's tau, defined to be

\tau := \frac{(\# \text{concordant pairs}) - (\# \text{discordant pairs})}{(\# \text{concordant pairs}) + (\# \text{discordant pairs})},

where a concordant pair is a pair of points such that if the rank of the x values is higher for the first (second) point of the pair, so too the rank of the y value is higher for the first (second) point of the pair.

The Bayesian approach of Chechile (2020) puts a Beta prior on phi, the proportion of concordance, i.e.,

\phi := \frac{(\# \text{concordant pairs})}{(\# \text{concordant pairs}) + (\# \text{discordant pairs})}.

The relationship between the two, then, is \tau = 2\phi - 1, or equivalently \phi = (\tau + 1)/2.

For more information, see dfba_bivariate_concordance and vignette("dfba_bivariate_concordance",package = "DFBA").

Value

An object of class b_procedure-class.

References

Chechile, R.A. (2020). Bayesian Statistics for Experimental Scientists: A General Introduction Using Distribution_Free Statistics. Cambridge: MIT Press.

Chechile, R.A., & Barch, D.H. (2021). A distribution-free, Bayesian goodness-of-fit method for assessing similar scientific prediction equations. Journal of Mathematical Psychology. https://doi.org/10.1016/j.jmp.2021.102638

Lindley, D. V., & Phillips, L. D. (1976). Inference for a Bernoulli process (a Bayesian view). The American Statistician, 30, 112-119.

Barch DH, Chechile RA (2023). DFBA: Distribution-Free Bayesian Analysis. doi:10.32614/CRAN.package.DFBA

Examples

# Generate data
set.seed(2025)
N = 50
x = rnorm(N)
y = x + 4 * rnorm(N)

# Test for non-zero correlation
cor_test_b(x,y)

# Input can be in the form of formula and data
cor_test_b(~ asdf + qwer,
           data = data.frame(asdf = x,
                             qwer = y))

# Other priors can be used, also.  See help for details.
cor_test_b(x,y,
           prior = "uniform")
cor_test_b(x,y,
           prior = "negative")
cor_test_b(x,y,
           prior = "positive")
cor_test_b(x,y,
           prior_shapes = c(10,10))

Credible Intervals for Model Parameters

Description

Computes credible intervals for one or more parameters in a fitted model.

Usage

credint(object, ...)

## S3 method for class 'lm_b'
credint(object, CI_level = 0.95, ...)

## S3 method for class 'aov_b'
credint(object, CI_level = 0.95, which = c("means", "pairwise"), ...)

Arguments

Value

Matrix of credible intervals

Examples

set.seed(2025)
N = 500
test_data = 
  data.frame(x1 = rep(letters[1:5],N/5))
test_data$outcome = 
  rnorm(N,-1 + 2 * (test_data$x1 %in% c("d","e")) )

# Fit 1-way ANOVA model
fit1 <-
  aov_b(outcome ~ x1,
        test_data,
        prior_mean_mu = 2,
        prior_mean_nu = 0.5,
        prior_var_shape = 0.01,
        prior_var_rate = 0.01)
credint(fit1)

Find Parameters for Beta Prior Based on Prior Mean and One Quantile

Description

Find Parameters for Beta Prior Based on Prior Mean and One Quantile

Usage

find_beta_parms(
  mean,
  quantile,
  left_tail_prob,
  plot_results = TRUE,
  search_bounds = c(0.001, 100)
)

Arguments

Value

Vector of beta shape parameters

Examples

find_beta_parms(2/5,0.68,0.9)
2/ (2 + 3)
qbeta(0.9,2,3)

Find Parameters for Inverse Gamma Prior Based on Prior Mean and One Quantile

Description

Find Parameters for Inverse Gamma Prior Based on Prior Mean and One Quantile

Usage

find_invgamma_parms(
  lower_quantile,
  upper_quantile,
  response_variance,
  lower_R2,
  upper_R2,
  probability,
  plot_results = TRUE
)

Arguments

Details

Either provide the lower and upper quantiles that contain probability of the inverse gamma distribution, or if this is for linear regression, you can specify that you are a priori probability sure that the coefficient of determination (R^2) falls within the two bounds provided, assuming that the residual variance is 1-R^2 times the total variance.

Value

twice the shape and rate of the inverse gamma distribution.

Examples

# When aimed at linear regression via coefficient of determination...
hypothetical_s2_y = 2.0
lower_R2 = 0.05
upper_R2 = 0.85
find_invgamma_parms(response_variance = hypothetical_s2_y,
                    lower_R2 = lower_R2,
                    upper_R2 = upper_R2,
                    probability = 0.8)

# More arbitrary task...
find_invgamma_parms(0.3, # hypothetical_s2_y * (1.0 - upper_R2)
                    1.9, #hypothetical_s2_y * (1.0 - lower_R2)
                    probability = 0.8)

Fractional Bayes Factors

Description

Compute fractional Bayes factors for lm_b objects

Usage

frac_bayes_factors(object1, object2, fractional_proportion)

Arguments

Details

Fractional Bayes factors, devised by O'Hagan, are a way to use flat, even improper, priors to obtain valid Bayes factors. The idea is built on the notion of partial Bayes factors, where a part of the data is used to determine the prior, and the remaining is used to compare the models.

References

O’Hagan, Anthony. “Fractional Bayes Factors for Model Comparison.” Journal of the Royal Statistical Society. Series B (Methodological), vol. 57, no. 1, 1995, pp. 99–138. https://doi.org/10.1111/j.2517-6161.1995.tb02017.x

Examples

set.seed(2026)
N = 500
test_data <-
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome <-
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )
fit_full <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data)
fit_no_x1 <-
  lm_b(outcome ~ x2 + x3,
       data = test_data)
fit_no_x2 <-
  lm_b(outcome ~ x1 + x3,
       data = test_data)

frac_bayes_factors(fit_full,
                   fit_no_x1)
frac_bayes_factors(fit_full,
                   fit_no_x2)

Get Posterior Samples from lm_b Object

Description

Get Posterior Samples from lm_b Object

Usage

get_posterior_draws(object, n_draws, seed, ...)

## S3 method for class 'lm_b'
get_posterior_draws(object, n_draws = 10000, seed = 1, ...)

## S3 method for class 'aov_b'
get_posterior_draws(object, n_draws = 10000, seed = 1, ...)

Arguments

Value

matrix of posterior draws

Examples

set.seed(2025)
N = 500
test_data <-
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome <-
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )
fit1 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data)
pdraws <-
  get_posterior_draws(fit1)
  

Bayesian Generalized Linear Models

Description

glm_b is used to fit linear models. It can be used to carry out regression for gaussian, binomial, and poisson data. Note that if the family is gaussian, this is just a wrapper for lm_b.

Usage

glm_b(
  formula,
  data,
  family,
  trials,
  prior = c("zellner", "normal", "improper"),
  zellner_g,
  prior_beta_mean,
  prior_beta_precision,
  prior_phi_mean = 1,
  ROPE,
  CI_level = 0.95,
  vb_maximum_iterations = 1000,
  algorithm = c("VB", "IS", "LSA"),
  proposal_df = 5,
  seed = 1,
  mc_error = 0.01,
  save_memory = FALSE
)

Arguments

Value

Object of class glm_b and lm_b.

Importance sampling:

glm_b will, unless use_importance_sampling = FALSE, perform importance sampling. The proposal will use a multivariate t distribution, centered at the posterior mode, with the negative hessian as its precision matrix. Do NOT treat the proposal_draws as posterior draws.

Priors:

If the prior is set to be either "zellner" or "normal", a normal distribution will be used as the prior of \beta, specifically

\beta \sim N(\mu, V)

where \mu is the prior_beta_mean and V is the prior_beta_precision (not covariance) matrix.

• zellner: glm_b sets \mu=0 and V = \frac{1}{g} X^{\top} X.

• normal: If missing prior_beta_mean, glm_b sets \mu=0, and if missing prior_beta_precision V will be a diagonal matrix. The first element, corresponding to the intercept, will be (2.5\times \max{\tilde{s}_y,1})^{-2}, where \tilde{s}_y is max of 1 and the standard deviation of y. Remaining diagonal elements will equal (2.5 s_y/s_x)^{-2}, where s_x is the standard deviations of the covariates. This equates to being 95% certain a priori that a change in x by one standard deviation (s_x) would not lead to a change in the linear predictor of more than 5 standard deviations (5s_y). This imposes weak regularization that adapts to the scale of the data elements.

ROPE:

If missing, the ROPE bounds will be given under the principle of "half of a small effect size."

• Gaussian. Using Cohen's D of 0.2 as a small effect size, the ROPE is built under the principle that moving the full range of X (i.e., \pm 2 s_x) will not move the mean of y by more than the overall mean of y minus 0.1s_y to the overall mean of y plus 0.1s_y. The result is a ROPE equal to |\beta_j| < 0.05s_y/s_j. If the covariate is binary, then this is simply |\beta_j| < 0.2s_y.

• Poisson. FDA guidance suggests a small effect is a rate ratio less than 1.25. We use half this effect: 1.125, and consider ROPE to indicate that a moving the full range of X (\pm 2s_x will not change the rate ratio by more than this amount. Thus the ROPE for the regression coefficient equals |\beta| < \frac{\log(1.125)}{4s_x}. For binary covariates, this is simply |\beta| < \log(1.125).

References

Kruschke JK. Rejecting or Accepting Parameter Values in Bayesian Estimation. Advances in Methods and Practices in Psychological Science. 2018;1(2):270-280. doi:10.1177/2515245918771304

Tim Salimans. David A. Knowles. "Fixed-Form Variational Posterior Approximation through Stochastic Linear Regression." Bayesian Anal. 8 (4) 837 - 882, December 2013. https://doi.org/10.1214/13-BA858

Examples

# Generate some negative-binomial data
set.seed(2025)
N = 500
test_data =
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5],
             time = rexp(N))
test_data$outcome =
  rnbinom(N,
          mu = exp(-2 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e"))) * test_data$time,
          size = 0.7)

# Fit using variational Bayes (default)
fit_vb1 <-
  glm_b(outcome ~ x1 + x2 + x3 + offset(log(time)),
        data = test_data,
        family = negbinom(),
        seed = 2025)
fit_vb1
summary(fit_vb1,
        CI_level = 0.9)
plot(fit_vb1)
coef(fit_vb1)
credint(fit_vb1,
        CI_level = 0.99)
bayes_factors(fit_vb1,
              by = "v")
preds =
  predict(fit_vb1)

# Try different priors
fit_vb2 <-
  glm_b(outcome ~ x1 + x2 + x3 + offset(log(time)),
        data = test_data,
        family = negbinom(),
        seed = 2025,
        prior = "normal")
fit_vb2
fit_vb3 <-
  glm_b(outcome ~ x1 + x2 + x3 + offset(log(time)),
        data = test_data,
        family = negbinom(),
        seed = 2025,
        prior = "improper")
fit_vb3

# Use Importance sampling instead of VB
fit_is <-
  glm_b(outcome ~ x1 + x2 + x3 + offset(log(time)),
        data = test_data,
        family = negbinom(),
        algorithm = "IS",
        seed = 2025)
summary(fit_is)

# Use large sample approximation instead of VB
fit_lsa <-
  glm_b(outcome ~ x1 + x2 + x3 + offset(log(time)),
        data = test_data,
        family = negbinom(),
        algorithm = "LSA",
        seed = 2025)
summary(fit_lsa)

Test for Heteroscedasticity in AOV Models

Description

Use Chib's method to compute the Bayes factor to test for heteroscedasticity in analysis of variance models.

Usage

heteroscedasticity_test(hetero_model, homo_model)

Arguments

Value

(returned invisible) A tibble with Bayes factors and interpretations.

References

Kass, R. E., & Raftery, A. E. (1995). Bayes Factors. Journal of the American Statistical Association, 90(430), 773–795.

Examples

# Test homoscedastic case
## Generate some data
set.seed(2025)
N = 200
test_data = 
  data.frame(x1 = rep(letters[1:5],N/5))
test_data$outcome = 
  rnorm(N,-1 + 2 * (test_data$x1 %in% c("d","e")) )

## Fit the anova models
hetero_model = 
  aov_b(outcome ~ x1,
        test_data)
homo_model = 
  aov_b(outcome ~ x1,
        test_data,
        heteroscedastic = FALSE)

## Perform test for heteroscedasticity using Bayes factors
heteroscedasticity_test(hetero_model,
                        homo_model)

# Test heteroscedastic case
## Generate some data
set.seed(2025)
N = 200
test_data = 
  data.frame(x1 = rep(letters[1:5],N/5))
test_data$outcome = 
  rnorm(N,
        -1 + 2 * (test_data$x1 %in% c("d","e")),
        sd = 3 - 2 * (test_data$x1 %in% c("d","e")))

## Fit the anova models
hetero_model = 
  aov_b(outcome ~ x1,
        test_data)
homo_model = 
  aov_b(outcome ~ x1,
        test_data,
        heteroscedastic = FALSE)

## Perform test for heteroscedasticity using Bayes factors
heteroscedasticity_test(hetero_model,
                        homo_model)

Bayesian Linear Models

Description

lm_b is used to fit linear models. It can be used to carry out regression, single stratum analysis of variance and analysis of covariance (although aov_b may provide a more convenient interface for ANOVA.)

Usage

lm_b(
  formula,
  data,
  weights,
  prior = c("zellner", "conjugate", "improper"),
  zellner_g,
  prior_beta_mean,
  prior_beta_precision,
  prior_var_shape,
  prior_var_rate,
  ROPE,
  CI_level = 0.95
)

Arguments

Details

MODEL:

The likelihood is given by

y_i \overset{ind}{\sim} N(x_i'\beta,\sigma^2).

The prior is given by

\beta|\sigma^2 \sim N\left( \mu, \sigma^2 V^{-1} \right) \\ \sigma^2 \sim \Gamma^{-1}(a/2,b/2).

• For Zellner's g prior, \frac{1}{g}X'X.

• The default for the conjugate prior is based on arguments from standardized regression. The default V is set according to the following: "a priori, we are 95% certain that a standard deviation increase in X will not lead to more than a 5 standard deviation in the mean of y." This leads to

  V^{1/2} = diag(v^{1/2},s_{x_1},\ldots,s_{x_p}) / 2.5,

  v^{-1} := \sum \frac{\bar x_j^2}{s_j^2}

  where s_y is the standard deviation of y, and s_{x_j} is the standard deviation of the j^{th} covariate.

• Unless prior_var_shape AND prior_var_rate are provided, the inverse gamma prior on the residual variance will place 50% prior probability that the correlation between the response and the fitted values is between 0.1 and 0.9.

• If prior = "improper", then the prior is

  \pi(\beta,\sigma^2) \propto \frac{1}{\sigma^2}.

ROPE:

If missing, the ROPE bounds will be given under the principle of "half of a small effect size." Using Cohen's D of 0.2 as a small effect size, the ROPE is built under the principle that moving the full range of X (i.e., \pm 2 s_x) will not move the mean of y by more than the overall mean of y minus 0.1s_y to the overall mean of y plus 0.1s_y. The result is a ROPE equal to |\beta_j| < 0.05s_y/s_j. If the covariate is binary, then this is simply |\beta_j| < 0.2s_y.

Value

Object of class lm_b.

Examples

# Generate some data
set.seed(2025)
N = 500
test_data = 
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome = 
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )

# Find good hyperparameters for the residual variance
s2_hyperparameters = 
  find_invgamma_parms(lower_R2 = 0.05,
                      upper_R2 = 0.95,
                      probability = 0.8,
                      response_variance = var(test_data$outcome))
## Check (should equal ~ 0.8)
extraDistr::pinvgamma((1.0 - 0.05) * var(test_data$outcome),
                      0.5 * s2_hyperparameters[1],
                      0.5 * s2_hyperparameters[2]) -
  extraDistr::pinvgamma((1.0 - 0.95) * var(test_data$outcome),
                        0.5 * s2_hyperparameters[1],
                        0.5 * s2_hyperparameters[2])

# Fit the linear model using a conjugate prior
fit1 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data,
       prior = "conj",
       prior_var_shape = s2_hyperparameters["shape"],
       prior_var_rate = s2_hyperparameters["rate"])
fit1
summary(fit1,
        CI_level = 0.99)
plot(fit1)
coef(fit1)
credint(fit1,
        CI_level = 0.9)
bayes_factors(fit1)
bayes_factors(fit1,
              by = "var")
AIC(fit1)
BIC(fit1)
DIC(fit1)
WAIC(fit1)
vcov(fit1)
preds = predict(fit1)

# Try other priors
fit2 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data) # Default is prior = "zellner"
fit3 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data,
       prior = "improper")

lm_b Objects

Description

Objects of Class lm_b

Details

An object of class lm_b contains at least the following:

summary

Tibble giving the summary of the model parameters of having a minimum of:

Variable

character

Post Mean

Numeric

Lower

Numeric

Upper

Numeric

Prob Dir

Numeric

formula

data

A tibble containing the data used in the analysis.

CI_level

Numeric scalar giving the credible interval level as provided by the user.

fitted

Vector of fitted values

residuals

Vector of Pearson residuals

family

xlevels

Named list, giving the levels for each factor covariate.

model_type

Character, either "parametric" or "nonparametric".

Further, one and only one of the following must be contained:

posterior_covariance

matrix

importance_sampling_weights,proposal_draws

Numeric and matrix, respectively

posterior_draws

matrix

S3 methods

Methods are available for print() and plot(), depending on which components are present.

See Also

print.b_procedure, plot.b_procedure

Examples

## Not run: 
cc_fit <- case_control_b(matrix(c(8,47,1,26),2,2))
cc_fit
plot(cc_fit)

## End(Not run)

Extract Log-Likelihood

Description

Computes the log-likelihood for fitted model objects.

Usage

## S3 method for class 'lm_b'
logLik(object, ...)

Arguments

Value

An object of class "logLik" with attributes "df" and "nobs".

See Also

logLik

Mediation using Bayesian Methods

Description

Mediation analysis done in the framework of Imai et al. (2010).

Usage

mediate_b(
  model_m,
  model_y,
  treat,
  control_value,
  treat_value,
  n_draws = 500,
  ask_before_full_sampling = TRUE,
  CI_level = 0.95,
  seed = 1,
  mc_error = ifelse("glm_b" %in% model_y, 0.01, 0.002),
  batch_size = 500
)

Arguments

Details

The model is the same as that of Imai et al. (2010):

M_i(X) = w_i'\alpha_m + X\beta_m + \epsilon_{m,i}, \\ y_i(X, M(\tilde X)) = w_i'\alpha_y + X\beta_y + M(\tilde X)\gamma + \epsilon_{y,i}, \\ \epsilon_{m,i} \overset{iid}{\sim} N(0,\sigma^2_m), \\ \epsilon_{y,i} \overset{iid}{\sim} N(0,\sigma^2_y), \\

where M_i(X) is the mediator as a function of the treatment variable X, and w_i are confounder covariates.

Unlike the mediation R package, the estimation in mediate_b is fully Bayesian (as opposed to "quasi-Bayesian").

Value

A list with the following elements:

• summary - tibble giving results for causal mediation quantities

• posterior_draws (of counterfactual expectations)

• mc_error absolute error used, including any rescaling to match the scale of the outcome

• other inputs to mediate_b

References

Imai, Kosuke, et al. “A General Approach to Causal Mediation Analysis.” Psychological Methods, vol. 15, no. 4, 2010, pp. 309–34, https://doi.org/10.1037/a0020761.

Examples

# Simplest case
## Generate some data
set.seed(2025)
N = 500
test_data = 
  data.frame(tr = rnorm(N),
             x1 = rnorm(N))
test_data$m = 
  rnorm(N, 0.4 * test_data$tr - 0.25 * test_data$x1)
test_data$outcome = 
  rnorm(N,-1 + 0.6 * test_data$tr + 1.5 * test_data$m + 0.25 * test_data$x1)

## Fit the mediator and outcome models
m1 = 
  lm_b(m ~ tr + x1,
       data = test_data)
m2 = 
  lm_b(outcome ~ m + tr + x1,
       data = test_data)
## Estimate the causal mediation quantities
m3 <-
  mediate_b(m1,m2,
            treat = "tr",
            control_value = -2,
            treat_value = 2,
            n_draws = 500,
            mc_error = 0.05,
            ask_before_full_sampling = FALSE)
m3
summary(m3,
        CI_level = 0.9)

# More complicated scenario
## Generate some data
set.seed(2025)
N = 500
test_data = 
  data.frame(tr = rep(0:1,N/2),
             x1 = rnorm(N))
test_data$m = 
  rnorm(N, 0.4 * test_data$tr - 0.25 * test_data$x1)
test_data$outcome = 
  rpois(N,exp(-1 + 0.6 * test_data$tr + 1.5 * test_data$m + 0.25 * test_data$x1))

## Fit the mediator and outcome models
m1 = 
  lm_b(m ~ tr + x1,
       data = test_data)
m2 = 
  glm_b(outcome ~ m + tr + x1,
        data = test_data,
        family = poisson())

##  Estimate the causal mediation quantities
m3 <-
  mediate_b(m1,m2,
            treat = "tr",
            control_value = 0,
            treat_value = 1,
            n_draws = 500,
            mc_error = 0.05,
            ask_before_full_sampling = FALSE)
summary(m3)

Negative-Binomial Family

Description

The negbinom() is an additional family to be considered alongside others documented under stats::family.

Usage

negbinom()

Value

an object of class "family". See stats::family.

Examples

negbinom()

Non-parametric Linear Models

Description

np_glm_b uses general Bayesian inference with loss-likelihood bootstrap. This is, as implemented here, a Bayesian non-parametric linear models inferential engine. Applicable data types are continuous (use family = gaussian()), count (use family = poisson()), or binomial (use family = binomial()).

Usage

np_glm_b(
  formula,
  data,
  family,
  loss = "selfinformation",
  loss_gradient,
  trials,
  n_draws,
  ask_before_full_sampling = TRUE,
  CI_level = 0.95,
  ROPE,
  seed = 1,
  mc_error = 0.01
)

Arguments

Details

Consider a population parameter of interest defined in terms of minimizing a loss function \ell wrt the population distribution:

\theta(F_y) := \underset{\theta\in\Theta}{\text{argmax}} \int \ell(\theta,y)dF_y

If we use a non-parametric Dirichlet process prior on the distribution of y, F_y, and let the concentration parameter go to zero, we have the Bayesian bootstrap applied to a general Bayesian updating framework dictated by the loss function.

By default, the loss function is the self-information loss, i.e., the negative log likelihood. This then resembles a typical glm_b implementation, but is more robust to model misspecification.

Value

Object of class np_glm_b and lm_b.

References

S P Lyddon, C C Holmes, S G Walker, General Bayesian updating and the loss-likelihood bootstrap, Biometrika, Volume 106, Issue 2, June 2019, Pages 465–478, https://doi.org/10.1093/biomet/asz006

Examples

# Generate some data
set.seed(2025)
N = 500
test_data = 
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome = 
  rbinom(N,1,1.0 / (1.0 + exp(-(-2 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) ))))

# Fit the GLM via the (non-parametric) loss-likelihood bootstrap.
fit1 <-
  np_glm_b(outcome ~ x1 + x2 + x3,
           data = test_data,
           family = binomial())
fit1
summary(fit1,
        CI_level = 0.99)
plot(fit1)
coef(fit1)
credint(fit1)
predict(fit1,
        newdata = fit1$data[1,])
vcov(fit1)

Plots bayesics Objects.

Description

Plots bayesics Objects.

Usage

## S3 method for class 'lm_b'
plot(
  x,
  type = c("diagnostics", "cred band", "pred band"),
  statistic,
  mc_error = 0.005,
  seed = 1,
  variable,
  exemplar_covariates,
  combine_pred_cred = TRUE,
  variable_seq_length = 30,
  return_as_list = FALSE,
  CI_level = 0.95,
  PI_level = 0.95,
  backtransformation = function(x) {
     x
 },
  ...
)

## S3 method for class 'mediate_b'
plot(
  x,
  type = c("diagnostics", "acme", "ade"),
  statistic = list(m = NULL, y = NULL),
  return_as_list = FALSE,
  seed = 1,
  mc_error = 0.005,
  ...
)

## S3 method for class 'survfit_b'
plot(x, n_draws = 10000, seed = 1, CI_level = 0.95, ...)

## S3 method for class 'b_procedure'
plot(x, ...)

Arguments

Value

If return_as_list=TRUE, a list of requested ggplots.

Examples

set.seed(2025)
N = 500
test_data <-
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome <-
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )
fit1 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data)
plot(fit1)

Plot credible and prediction bands

Description

Plot credible and prediction bands

Usage

plot_bands(
  x,
  type,
  combine_pred_cred,
  CI_level,
  PI_level,
  backtransformation,
  return_as_list,
  ...
)

## S3 method for class 'lm_b'
plot_bands(
  x,
  type = c("cred band", "pred band"),
  combine_pred_cred = TRUE,
  CI_level = 0.95,
  PI_level = 0.95,
  backtransformation = function(x) {
     x
 },
  return_as_list = TRUE,
  variable,
  variable_seq_length = 30,
  exemplar_covariates,
  ...
)

## S3 method for class 'aov_b'
plot_bands(
  x,
  type = c("cred band", "pred band"),
  combine_pred_cred = TRUE,
  CI_level = 0.95,
  PI_level = 0.95,
  backtransformation = function(x) {
     x
 },
  return_as_list = TRUE,
  ...
)

Arguments

Diagnostic Plots for Bayesian Regression Objects

Description

Diagnostic Plots for Bayesian Regression Objects

Usage

plot_dx(x, statistic, mc_error, seed, return_as_list, ...)

## S3 method for class 'lm_b'
plot_dx(x, statistic, mc_error = 0.005, seed = 1, return_as_list = TRUE, ...)

## S3 method for class 'aov_b'
plot_dx(x, statistic, mc_error = 0.005, seed = 1, return_as_list = TRUE, ...)

## S3 method for class 'mediate_b'
plot_dx(
  x,
  statistic = list(m = NULL, y = NULL),
  mc_error = 0.005,
  seed = 1,
  return_as_list = TRUE,
  ...
)

Arguments

Poisson Procedures

Description

Make inference on one or two populations using Poisson distributed count data

Usage

poisson_test_b(
  x,
  offset,
  r,
  ROPE,
  prior = c("jeffreys", "flat"),
  prior_shape_rate,
  CI_level = 0.95,
  plot = TRUE,
  seed = 1,
  mc_error = 0.002
)

Arguments

Details

The likelihood is

y \sim Poi(\lambda t),

where \lambda is the rate, and t is the time or area observed and is given by the argument offset.

The prior is given by

\lambda \sim \Gamma(a,b),

where a and b are given by the argument prior_shape_rate. If prior_shape_rate is missing and prior = "jeffreys", then a Jeffrey's prior will be used, i.e., \Gamma(0.5,0) (improper), while if prior = "flat", \Gamma(0.001,0.001) will be used.

Value

An object of class b_procedure-class.

Examples

# One sample
poisson_test_b(x = 12)
## You can compute the posterior probability that the rate is less than r
poisson_test_b(x = 12,
               r = 8)

# Two samples
poisson_test_b(x = c(12,20))

# Offsets can be included:
poisson_test_b(x = c(12,20),
               offset = c(10,9))

# Different priors can be used
poisson_test_b(x = c(12,20),
               prior = "flat")
poisson_test_b(x = c(12,20),
               prior_shape_rate = c(20,1.5))

Predict Method for lm_b Model Fits

Description

Predict Method for lm_b Model Fits

Usage

## S3 method for class 'lm_b'
predict(
  object,
  newdata,
  trials,
  CI_level = 0.95,
  PI_level = 0.95,
  seed = 1,
  n_draws = 5000,
  ...
)

## S3 method for class 'aov_b'
predict(object, CI_level = 0.95, PI_level = 0.95, ...)

Arguments

Value

tibble with estimate (posterior mean), prediction intervals, and credible intervals for the mean.

Examples

# lm_b
## Create data
N = 500
test_data <-
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome <-
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )

## Fit linear model
fit <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data)
predict(fit)


# glm_b
## Generate some negative binomial data
set.seed(2025)
N = 500
test_data =
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5],
             time = rexp(N))
test_data$outcome =
  rnbinom(N,
          mu = exp(-2 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e"))) * test_data$time,
          size = 0.7)

## Fit using variational Bayes (default)
fit_vb1 <-
  glm_b(outcome ~ x1 + x2 + x3 + offset(log(time)),
        data = test_data,
        family = negbinom(),
        seed = 2025)
# Predict
predict(fit_vb1)

## Fit the GLM via the (non-parametric) loss-likelihood bootstrap.
fit_np <-
  np_glm_b(outcome ~ x1 + x2 + x3 + offset(log(time)),
           data = test_data,
           family = negbinom())
predict(fit_np)


# bma_inference
## Create data
set.seed(2025)
N = 500
test_data = 
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5],
             x4 = rnorm(N),
             x5 = rnorm(N),
             x6 = rnorm(N),
             x7 = rnorm(N),
             x8 = rnorm(N),
             x9 = rnorm(N),
             x10 = rnorm(N))
test_data$outcome = 
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )

## Fit linear model using Bayesian model averaging
fit <-
  bma_inference(outcome ~ .,
                test_data,
                user.int = FALSE)
predict(fit)

Print bayesics Objects.

Description

Print bayesics Objects.

Usage

## S3 method for class 'aov_b'
print(x, ...)

## S3 method for class 'lm_b'
print(x, ...)

## S3 method for class 'mediate_b'
print(x, ...)

## S3 method for class 'survfit_b'
print(x, ...)

## S3 method for class 'b_procedure'
print(x, ...)

Arguments

Value

None

Examples

set.seed(2025)
N = 500
test_data <-
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome <-
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )
fit1 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data)
print(fit1)

Binomial Procedures

Description

prop_test_b either makes inference on a single population proportion, or else compares two population proportions. binom_test_b is the same as prop_test_b.

Usage

prop_test_b(
  n_successes,
  n_failures,
  n_total,
  p,
  ROPE,
  prior = c("jeffreys", "uniform"),
  prior_shapes,
  CI_level = 0.95,
  plot = TRUE,
  seed = 1,
  mc_error = 0.002
)

Arguments

Details

The likelihood is given by

y \sim \text{Binom}(n,p),

and the prior on p is

p \sim Beta(a,b),

where a and b are given by the argument prior_shapes. If prior_shapes is missing and prior = "jeffreys", then a Jeffreys prior will be used (Beta(1/2,1/2)), and if prior = "uniform", then a uniform prior will be used (Beta(1,1)).

Value

An object of class b_procedure-class.

Examples

# Single population
prop_test_b(14,
            19)
# or another way of the same thing;
prop_test_b(14,
            n_total = 14 + 19)

# A null value compared against can be added:
prop_test_b(14,
            19,
            p = 0.5)

# Two populations
prop_test_b(c(14,22),
            c(19,45))
# or equivalently
prop_test_b(c(14,22),
            n_total = c(14,22) + c(19,45))

Paired Sign Test

Description

Sign test for paired data.

Usage

sign_test_b(
  x,
  y,
  p0 = 0.5,
  prior = c("jeffreys", "uniform"),
  prior_shapes,
  ROPE,
  CI_level = 0.95,
  plot = TRUE
)

Arguments

Details

The sign test looks at z_i:= 1_{[x_i > y_i]} rather than trying to model the distribution of (x_i,y_i). sign_test_b then uses the fact that

z_i \overset{iid}{\sim} Bernoulli(p)

and then makes inference on p. The prior on p is

p \sim Beta(a,b),

where a and b are given by the argument prior_shapes. If prior_shapes is missing and prior = "jeffreys", then a Jeffreys prior will be used (Beta(1/2,1/2)), and if prior = "uniform", then a uniform prior will be used (Beta(1,1)).

Value

An object of class b_procedure-class.

Examples

# Single population
sign_test_b(x = rnorm(50))

## Change ROPE
sign_test_b(x = rnorm(50),
            ROPE = 0.1)

## Change the prior
sign_test_b(x = rnorm(50),
            prior = "uniform")
sign_test_b(x = rnorm(50),
            prior_shapes = c(2,2))

## Two populations
sign_test_b(x = rnorm(50,1),
            y = rnorm(50,0))

## Change reference probability
sign_test_b(x = rnorm(50,1),
            y = rnorm(50,0),
            p0 = 0.7)

Summary functions for bayesics objects

Description

Summary functions for bayesics objects

Usage

## S3 method for class 'lm_b'
summary(
  object,
  CI_level = 0.95,
  interpretable_scale = TRUE,
  print_results = TRUE,
  ...
)

## S3 method for class 'aov_b'
summary(object, CI_level = 0.95, print_results = TRUE, ...)

## S3 method for class 'mediate_b'
summary(object, CI_level = 0.95, print_results = TRUE, ...)

Arguments

Value

tibble with summary values

Examples

set.seed(2025)
N = 500
test_data <-
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome <-
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )
fit1 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data)
summary(fit1)

Create Survival Curves

Description

Use the semi-parametric piecewise exponential survival model to fit a survival curve to one or more samples

Usage

survfit_b(formula, data, prior_shape, prior_rate, max_n_time_bins, n_time_bins)

Arguments

Details

The approach proposed by Qing et al. (2023) models the survival curve by way of piecewise exponential curves. That is, the hazard function is a piecewise function. The prior on the hazard within each "piece", or equivalently the rate of the exponential distribution, is a conjugate gamma distribution. Unless specified, the prior shape and rate for each piece is the posterior under the assumption that the data follow a single exponential distribution.

Unless prespecified by the user, the number of breaks in the hazard function is determined by Bayes factors, which can be quickly computed analytically.

If more than one population is being compared, then as before Bayes factors will be used to determine the number of breaks in each group's hazard function, and then Bayes factors will be used to compare the hypothesis that each group has a separate survival function vs. the null hypothesis that all groups share the same survival function.

Value

Object of class survfit_b with the following:

• posterior_parameters An n_time_binsx2 matrix whose columns provide shapes and rates of the gamma posterior distribution of each of the piecewise hazard rates.

• intervals An n_time_binsx2 matrix whose columns provide the start and endpoints of each time bin. If comparing multiple samples, a list of such matrices will be provided.

• marginal_likelihood

• data

If comparing multiple samples, each group will have a list of posterior_parameters and intervals.

References

Qing Y, Thall PF, Yuan Y. A Bayesian piecewise exponential phase II design for monitoring a time-to-event endpoint. Pharm Stat. 2023 Jan;22(1):34-44. doi: 10.1002/pst.2256. Epub 2022

Examples

# Single population
set.seed(2025)
N = 300
test_data = 
  data.frame(outcome = 
               rweibull(N,2,5))
test_data$observed = 
  ifelse(test_data$outcome >= 7, 0, 1)
test_data$outcome =
  ifelse(dplyr::near(test_data$observed,1), test_data$outcome, 7)
fit1 = 
  survfit_b(Surv(test_data$outcome,
                 test_data$observed) ~ 1)
fit1
plot(fit1)

# Multiple populations
set.seed(2025)
N = 300
test_data = 
  data.frame(outcome = 
               c(rweibull(2*N/3,2,5),
                 rweibull(N/3,2,10)),
             x1 = rep(letters[1:3],each = N/3))
test_data$observed = 
  ifelse(test_data$outcome >= 9, 0, 1)
test_data$outcome =
  ifelse(dplyr::near(test_data$observed,1), test_data$outcome, 9)
fit2 =
  survfit_b(Surv(outcome,
                 observed) ~ x1,
            data = test_data)
fit2
plot(fit2)

t-test

Description

One and two sample t-tests on vectors of data

Usage

t_test_b(
  x,
  y,
  mu,
  paired = FALSE,
  data,
  heteroscedastic = TRUE,
  prior_mean_mu,
  prior_mean_nu = 0.001,
  prior_var_shape = 0.001,
  prior_var_rate = 0.001,
  CI_level = 0.95,
  ROPE = 0.1,
  improper = FALSE,
  plot = TRUE,
  seed = 1,
  mc_error = 0.002
)

Arguments

Details

A one and two sample t-test is nothing more than a special case of one-way anova. See aov_b for details.

Value

An object of class b_procedure-class.

Examples

# Single population
t_test_b(rnorm(50))
# or an alternative input format
t_test_b(outcome ~ 1,
         data = data.frame(outcome = rnorm(50)))

# Two populations
t_test_b(rnorm(50),
         rnorm(15,1))

# or an alternative input format
t_test_b(outcome ~ group_variable,
         data = 
           data.frame(outcome = c(rnorm(50),
                                  rnorm(15,1)),
                      group_variable = rep(c("a","b"),
                                           c(50,15))))

Calculate Posterior Variance-Covariance Matrix for a Bayesian Fitted Model Object

Description

Calculate Posterior Variance-Covariance Matrix for a Bayesian Fitted Model Object

Usage

## S3 method for class 'lm_b'
vcov(object, ...)

Arguments

Value

A matrix of the covariance matrix for the regression coefficients. If the posterior is a multivariate t distribution (or consists of independent t's in the case of heteroscedastic 1-way ANOVA), the degrees of freedom are returned as the df attribute of the matrix. Note that for lm_b and aov_b objects, this function already takes into account the uncertainty around the residual variance.

Examples

set.seed(2025)
N = 500
test_data <-
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome <-
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )
fit1 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data)
vcov(fit1)

Bayesian Wilcoxon Rank Sum (aka Mann-Whitney U) and Signed Rank Analyses

Description

Bayesian Wilcoxon Rank Sum (aka Mann-Whitney U) and Signed Rank Analyses

Usage

wilcoxon_test_b(
  x,
  y,
  paired = FALSE,
  p = 0.5,
  ROPE,
  prior = c("centered", "uniform"),
  prior_shapes,
  CI_level = 0.95,
  plot = TRUE,
  seed = 1
)

Arguments

Details

Bayesian Wilcoxon signed rank analysis For a single input vector or paired data, the Bayesian signed rank analysis will be performed. The estimand is the proportion of (differenced) values that are positive. For more information, see dfba_wilcoxon and vignette("dfba_wilcoxon",package = "DFBA").

Bayesian Wilcoxon rank sum/Mann-Whitney analysis For unpaired x and y inputs, the Bayesian rank sum analysis will be performed. The estimand is \Omega_x:=\lim_{n\to\infty} \frac{U_x}{U_x + U_y}, where U_x is the number of pairs (i,j) such that x_i > y_j, and vice versa for U_y. That is, it is the population proportion of all untied pairs for which x > y. Larger values imply that x is stochastically larger than y. For more information, see dfba_mann_whitney and vignette("dfba_mann_whitney",package = "DFBA").

Value

An object of class b_procedure-class.

References

Chechile, R.A. (2020). Bayesian Statistics for Experimental Scientists: A General Introduction to Distribution-Free Methods. Cambridge: MIT Press.

Chechile, R. A. (2018) A Bayesian analysis for the Wilcoxon signed-rank statistic. Communications in Statistics - Theory and Methods, https://doi.org/10.1080/03610926.2017.1388402

Chechile, R.A. (2020). A Bayesian analysis for the Mann-Whitney statistic. Communications in Statistics – Theory and Methods 49(3): 670-696. https://doi.org/10.1080/03610926.2018.1549247.

Barch DH, Chechile RA (2023). DFBA: Distribution-Free Bayesian Analysis. doi:10.32614/CRAN.package.DFBA

Examples

# Signed rank analysis
## Generate data
N = 150
set.seed(2025)
test_data = 
  data.frame(x = rbeta(N,2,10),
             y = rbeta(N,5,10))

## input differenced data
wilcoxon_test_b(test_data$x - test_data$y)
## input paired data vectors individually
wilcoxon_test_b(test_data$x,
                test_data$y,
                paired = TRUE)

## Use different priors
wilcoxon_test_b(test_data$x - test_data$y,
                prior = "uniform")
wilcoxon_test_b(test_data$x - test_data$y,
                prior_shapes = c(5,5))

## Change ROPE bounds
wilcoxon_test_b(test_data$x - test_data$y,
                ROPE = 0.1)

# Rank sum analysis
## Generate data
set.seed(2025)
N = 150
x = rbeta(N,2,10)
y = rbeta(N + 1,5,10)

## Perform analysis
wilcoxon_test_b(x,y)