Package {smoothbp}

Title:	Hierarchical Piecewise Regression with Smoothed Change-Points
Version:	0.2.1
Description:	Fits Bayesian hierarchical piecewise regression models with multiple logistic-smoothed change-points. Non-linear parameters (change-point locations and transition sharpness) and linear parameters can each be conditioned on covariates and factors via flexible design matrices. A random-intercept structure is supported for any parameter. Spike-and-slab regularization is supported for selecting the number of breakpoints. Posterior inference uses a Metropolis-within-Gibbs sampler implemented in 'Rust' for speed. Methods are based on the smooth transition piecewise regression model of Bacon and Watts (1971) <doi:10.2307/2334389> and variable selection spike-and-slab priors of Kuo and Mallick (1998) https://www.jstor.org/stable/25053023.
License:	MIT + file LICENSE
Encoding:	UTF-8
RoxygenNote:	7.3.3
Config/rextendr/version:	0.5.0
SystemRequirements:	Cargo (Rust's package manager), rustc >= 1.65.0
Depends:	R (≥ 4.2)
Imports:	posterior, ggplot2, stats, bridgesampling, loo, bayesplot
Suggests:	testthat (≥ 3.0.0), withr, knitr, rmarkdown, brms, mcp, dplyr, gt, tidyr, scales, rjags
Config/testthat/edition:	3
VignetteBuilder:	knitr
URL:	https://github.com/ABindoff/smoothbp
BugReports:	https://github.com/ABindoff/smoothbp/issues
NeedsCompilation:	yes
Packaged:	2026-05-27 21:44:34 UTC; bindoffa
Author:	Aidan D Bindoff [aut, cre]
Maintainer:	Aidan D Bindoff <aidan.bindoff@utas.edu.au>
Repository:	CRAN
Date/Publication:	2026-06-01 08:30:02 UTC

Convert draws to a data frame

Description

Convert draws to a data frame

Usage

## S3 method for class 'smoothbp_fit'
as.data.frame(x, ...)

Arguments

Value

A data.frame containing the posterior draws of the model parameters.

Bayes Factor for smoothbp_fit

Description

Bayes Factor for smoothbp_fit

Usage

## S3 method for class 'smoothbp_fit'
bayes_factor(x1, x2, log = FALSE, ...)

Arguments

Value

A numeric value representing the Bayes Factor (or log Bayes Factor if log = TRUE) comparing x1 to x2.

Bridge Sampler for smoothbp_fit

Description

Bridge Sampler for smoothbp_fit

Usage

## S3 method for class 'smoothbp_fit'
bridge_sampler(
  samples,
  method = c("auto", "rust", "bridgesampling"),
  seed = 42,
  ...
)

Arguments

Value

An object of class "bridge" or "bridge_list" containing the log marginal likelihood estimate.

Fitted values for smoothbp_fit objects

Description

Fitted values for smoothbp_fit objects

Usage

## S3 method for class 'smoothbp_fit'
fitted(object, newdata = NULL, summary = TRUE, ...)

Arguments

Value

If summary = TRUE, a data.frame containing the observation index, mean fitted value, and 95% CI bounds. If summary = FALSE, a matrix of dimension S x N where S is the number of posterior draws and N is the number of observations, containing the posterior draws of the fitted values.

Fix a parameter at a specific value

Description

Used within omega or rho lists in smoothbp() to specify that a parameter is fixed and should not be estimated.

Usage

fixed(value)

Arguments

Value

A smoothbp_fixed object.

Test hypotheses and compute evidence ratios from posterior draws

Description

Evaluates one or more hypothesis expressions against the posterior draws of a smoothbp_fit object, returning posterior probabilities and evidence ratios.

Usage

hypothesis(object, hypotheses, ci = 0.95, ...)

Arguments

Value

An object of class c("smoothbp_hypothesis", "data.frame") with one row per hypothesis and columns:

Hypothesis

The original hypothesis string.

Estimate

Posterior mean of the contrast.

Est.Error

Posterior SD of the contrast.

CI.lower, CI.upper

Credible interval bounds.

P(H)

Posterior probability of the hypothesis.

Evid.Ratio

Evidence ratio P(H)/(1-P(H)).

Star

Informal star coding based on the evidence ratio.

Hypothesis syntax

Write hypotheses as character strings using exact parameter names as they appear in fit$param_names (e.g. "b2_(Intercept)", "omega_(Intercept)"). No backtick-quoting is needed; the function handles special characters internally.

Two forms are accepted:

Directional

An expression containing >, <, >=, or <=. The hypothesis is evaluated as a contrast: the left-hand side minus the right-hand side (direction-adjusted), and P(H \mid \text{data}) is the proportion of posterior draws satisfying the condition. Examples: "b2_(Intercept) > 0", "omega_(Intercept) < 4", "b2_(Intercept) - b1_(Intercept) > 0".

Contrast

A numeric expression without a comparison operator. The function summarises the posterior distribution of the derived quantity and reports P(\text{expression} > 0) as the directional probability. Example: "b2_(Intercept) - b1_(Intercept)".

Point-null hypotheses (==) are not supported because they require the Savage-Dickey density ratio; use bayestestR::rope() for interval-based equivalence testing instead.

Evidence ratio interpretation

ER = \frac{P(H \mid \text{data})}{1 - P(H \mid \text{data})}

An ER of 19 corresponds to P(H) = 0.95; ER = 1 means the posterior is equally split. Star codes: *** ER > 99, ** ER > 19, * ER > 3.

Examples

## Not run: 
# Is the change in slope positive?
hypothesis(fit, "b2_(Intercept) > 0")

# Does the change-point fall before time 4?
hypothesis(fit, "omega_(Intercept) < 4")

# Multiple hypotheses at once
hypothesis(fit, c(
  "b2_(Intercept) > 0",
  "omega_(Intercept) < 4",
  "b2_(Intercept) - b1_(Intercept) > 0"
))

# Posterior summary of a contrast (no comparison operator)
hypothesis(fit, "b2_(Intercept) - b1_(Intercept)")

## End(Not run)

Pointwise log-likelihood matrix

Description

Pointwise log-likelihood matrix

Usage

log_lik(object, ...)

## S3 method for class 'smoothbp_fit'
log_lik(object, ...)

Arguments

Value

A matrix of pointwise log-likelihood values of dimension S x N, where S is the number of posterior draws and N is the number of observations.

Posterior inclusion probabilities from a spike-and-slab fit

Description

Posterior inclusion probabilities from a spike-and-slab fit

Usage

pip(object, ...)

Arguments

Value

A named numeric vector of posterior inclusion probabilities (PIPs) for each b2 coefficient that had a spike-and-slab prior.

Trace and density plots for a smoothbp_fit

Description

A thin wrapper around trace_plot for the standard plot() interface.

Usage

## S3 method for class 'smoothbp_fit'
plot(x, type = "trace", pars = NULL, ...)

Arguments

Value

A ggplot object, or a named list of two when type = "both".

Plot posterior inclusion probabilities

Description

Plot posterior inclusion probabilities

Usage

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

Arguments

Value

A ggplot object.

Print a smoothbp_fit

Description

Print a smoothbp_fit

Usage

## S3 method for class 'smoothbp_fit'
print(x, digits = 3, effects = c("fixed", "ran_pars"), ...)

Arguments

Value

The input object x (invisibly), called for its printing side effects.

Specify a gamma prior for a parameter

Description

Specify a gamma prior for a parameter

Usage

prior_gamma(shape = 1, scale = 1)

Arguments

Value

A smoothbp_prior object.

Specify an inverse-gamma prior for a variance component

Description

Specify an inverse-gamma prior for a variance component

Usage

prior_invgamma(shape = 1, scale = 1)

Arguments

Value

A smoothbp_prior object.

Specify a normal (or truncated normal) prior for a regression coefficient

Description

Specify a normal (or truncated normal) prior for a regression coefficient

Usage

prior_normal(mean = 0, sd = 1, lb = -Inf, ub = Inf)

Arguments

Value

A smoothbp_prior object.

Specify a spike-and-slab prior for variable selection

Description

Used with smoothbp_ss() to place a point-mass spike at zero on selected coefficients.

Usage

prior_spike_slab(
  pi = 0.5,
  slab = prior_normal(0, 2),
  learn_pi = FALSE,
  a = 1,
  b = 1
)

Arguments

Value

A smoothbp_spike_slab object.

Parameter recovery plot

Description

Compares posterior estimates against the known data-generating values stored in attr(dat, "true_params") (as produced by simulate_smoothbp). For each matched parameter the plot shows the posterior mean, a credible interval, and the true value, coloured by whether the interval contains the truth.

Usage

recovery_plot(fit, dat, level = 0.95)

Arguments

Details

The function looks for population-level intercept parameters only. If the fitted model has covariates in a given component (e.g. omega = ~ 1 + group) the function still extracts the (Intercept) term for comparison against the scalar true value from the simulation.

Value

A ggplot object.

Examples

## Not run: 
dat <- simulate_smoothbp(n_subj = 20, n_obs = 8, seed = 42)
fit <- smoothbp(y ~ tau, b0 = ~ 1 + (1 | subject), data = dat,
                priors = smoothbp_priors(omega = prior_normal(3, 2, lb = 0)),
                chains = 4L, iter = 2000L, warmup = 1000L, seed = 42L)
recovery_plot(fit, dat)

## End(Not run)

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

bayesplot

pp_check

loo

loo, waic

Robustify a smoothbp_fit object using a Bayesian Sandwich approach

Description

Robustify a smoothbp_fit object using a Bayesian Sandwich approach

Usage

robustify(object, cluster, ...)

Arguments

Value

A new smoothbp_fit object where the MCMC draws have been affinely transformed to match the robust clustered covariance matrix.

Round numeric columns in a data frame

Description

Round numeric columns in a data frame

Usage

round_df(df, digits = 3)

Simulate data from the smooth change-point model

Description

Generates synthetic data from the model used by smoothbp, including optional between-subject random intercepts. True parameter values are stored as the "true_params" attribute so they can be compared against posterior estimates.

Usage

simulate_smoothbp(
  n_subj = 20L,
  n_obs = 8L,
  b0 = 5,
  b1 = -0.3,
  b2 = 1.2,
  omega = 3,
  rho = 4,
  sigma = 0.4,
  sigma_u = 0.5,
  tau_range = c(0, 6),
  seed = NULL
)

Arguments

Details

The data-generating model is:

y_{ij} = (b0 + u_j) + b1 \cdot d_{ij} + b2 \cdot d_{ij} \cdot \sigma(d_{ij} \cdot \rho) + \varepsilon_{ij}

where d_{ij} = \tau_{ij} - \omega and \sigma(\cdot) is the logistic sigmoid.

Value

A data.frame with columns:

subject

Subject identifier (factor).

tau

Time variable.

mu

Noise-free conditional mean \mu_{ij}.

y

Observed response.

The attribute "true_params" is a named list containing the data-generating values of b0, b1, b2, omega, rho, sigma, sigma_u, the vector of subject-level deviations u, and the seed used.

Examples

dat <- simulate_smoothbp(
  n_subj = 20, n_obs = 8,
  b0 = 5, b1 = -0.3, b2 = 1.2,
  omega = 3, rho = 4, sigma = 0.4, sigma_u = 0.5,
  seed = 42
)
head(dat)
attr(dat, "true_params")

Fit a hierarchical piecewise regression model with smoothed change-points

Description

Fit a hierarchical piecewise regression model with smoothed change-points

Usage

smoothbp(
  formula,
  b0 = ~1,
  b1 = ~1,
  deltas = list(~1),
  omega = list(~1),
  rho = list(~1),
  data,
  priors = smoothbp_priors(),
  chains = 4L,
  iter = 2000L,
  warmup = 1000L,
  seed = NULL,
  step_om = 0.3,
  step_rho = 0.3,
  target_accept = 0.8,
  cores = getOption("smoothbp.cores", 1L),
  hierarchical = NULL,
  .verbose = TRUE
)

Arguments

Value

A smoothbp_fit object.

Collect priors for all model parameters

Description

Each argument accepts either:

• A single prior_normal() applied to all coefficients of that parameter, or

• A named list mapping coefficient names (matching column names of the design matrix) to individual prior_normal() objects.

Usage

smoothbp_priors(
  b0 = prior_normal(0, 10),
  b1 = prior_normal(0, 2),
  deltas = prior_normal(0, 2),
  omega = prior_normal(3, 2, lb = 0),
  rho = prior_normal(3, 2, lb = 0),
  sigma = prior_invgamma(1, 1),
  sigma_u = prior_invgamma(1, 1),
  sigma_re_om = prior_invgamma(1, 1)
)

Arguments

Details

For multi-breakpoint models, deltas, omega, and rho can also be lists of prior specifications (one per breakpoint slot). If a single specification is provided, it is applied to all slots.

Value

A smoothbp_priors list.

Fit a smooth change-point model with spike-and-slab variable selection

Description

Fit a smooth change-point model with spike-and-slab variable selection

Usage

smoothbp_ss(
  formula,
  b0 = ~1,
  b1 = ~1,
  deltas = list(~1),
  omega = list(~1),
  rho = list(~1),
  data,
  priors = smoothbp_priors(),
  spike = prior_spike_slab(),
  b1_spike = FALSE,
  hierarchical = NULL,
  chains = 4L,
  iter = 2000L,
  warmup = 1000L,
  seed = NULL,
  step_om = 0.3,
  step_rho = 0.3,
  target_accept = 0.65,
  cores = getOption("smoothbp.cores", 1L),
  .verbose = TRUE
)

Arguments

Value

A smoothbp_fit object.

Generate evenly spaced priors for candidate breakpoints

Description

This helper function generates a list of prior_normal objects for omega (breakpoint locations) that are evenly spaced across the range of your time variable tau. This is highly recommended when using smoothbp_ss() to ensure the candidate breakpoints cover the entire domain without clumping.

Usage

space_omega_priors(K, tau_min, tau_max)

Arguments

Details

For hierarchical models where omega has random effects (e.g., ~ 1 + (1 | group)), this function automatically names the prior (Intercept) so it applies correctly to the global market mean, while the random effects are handled automatically by the sigma_re_om shrinkage variance.

Value

A list of length K containing prior specifications for omega.

Summarise a smoothbp_fit

Description

Summarise a smoothbp_fit

Usage

## S3 method for class 'smoothbp_fit'
summary(object, effects = c("fixed", "ran_pars"), digits = 3, ...)

Arguments

Value

A data.frame containing summary statistics (mean, SD, 2.5% and 97.5% quantiles, Rhat, and bulk/tail ESS) for the requested effects.

Summarise a smoothbp_ss_fit

Description

Returns a data frame of posterior summaries for selected parameters.

Usage

## S3 method for class 'smoothbp_ss_fit'
summary(object, effects = "all", digits = 3, ...)

Arguments

Value

A data frame with one row per selected parameter and columns variable, mean, sd, Q2.5, Q97.5, rhat, ess_bulk, ess_tail. For gamma (spike-and-slab) parameters, the mean column contains the posterior inclusion probability (PIP).

Fixed-effects table for smoothbp_fit objects

Description

Collects posterior summaries from one or more smoothbp_fit objects and displays them in a single table with parameters on rows and models in columns. All parameters present in any model are shown; models that do not include a given parameter display a dash.

Usage

tab_smoothbp(
  ...,
  labels = NULL,
  digits = 2,
  fmt = c("mean [CI]", "mean (SD)"),
  show_rhat = FALSE
)

Arguments

Value

A gt_tbl object (or a knitr_kable if gt is not installed).

Examples

## Not run: 
tab_smoothbp(m.ep.pw1, m.wo.pw1, m.la.pw1,
             labels = c("Episodic", "Working", "Language"))

## End(Not run)

Trace plots with automatic poor-mixing highlighting

Description

Produces per-parameter trace plots from a smoothbp_fit object. Parameters with \hat{R} > 1.05 are flagged with a light-red background and their panel labels include the \hat{R} value and a warning symbol. Parameters with low bulk-ESS (< 100) are further annotated.

Usage

trace_plot(
  fit,
  pars = NULL,
  type = "trace",
  rhat_thresh = 1.05,
  ess_thresh = 100
)

Arguments

Value

A ggplot object (or a named list of two when type = "both").

Print true parameters from a simulated dataset

Description

Convenience function to display the data-generating parameters stored in the "true_params" attribute of a dataset returned by simulate_smoothbp.

Usage

true_params(dat)

Arguments

Value

The true_params list, invisibly.

Update a fitted smoothbp model

Description

Re-fits the model, replacing any arguments supplied here with the corresponding values stored in the original fit for anything left unspecified.

Usage

## S3 method for class 'smoothbp_fit'
update(
  object,
  formula,
  b0,
  b1,
  omega,
  rho,
  data,
  priors,
  chains,
  iter,
  warmup,
  seed,
  step_om,
  step_rho,
  target_accept,
  cores,
  .verbose = TRUE,
  ...
)

Arguments

Value

A new smoothbp_fit object.