inst/doc/brms-comparison.R

In smoothbp: Hierarchical Piecewise Regression with Smoothed Change-Points

## ----include = FALSE----------------------------------------------------------
NOT_CRAN <- identical(tolower(Sys.getenv("NOT_CRAN")), "true")
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  eval     = NOT_CRAN,
  fig.width  = 7,
  fig.height = 4,
  dpi        = 96
)

## ----setup, message = FALSE, warning = FALSE, eval = TRUE---------------------
library(smoothbp)
library(posterior)
library(ggplot2)
library(dplyr)
library(tidyr)

have_brms <- requireNamespace("brms", quietly = TRUE) && NOT_CRAN
if (!have_brms && NOT_CRAN) {
  message("brms not installed -- install.packages('brms') to render the comparison.")
}

have_mcp <- requireNamespace("mcp",   quietly = TRUE) &&
            requireNamespace("rjags", quietly = TRUE) && NOT_CRAN
if (!have_mcp && NOT_CRAN) {
  message("mcp or rjags not installed -- install JAGS then install.packages('mcp') ",
          "to render the mcp comparison.")
}

## ----simulate-----------------------------------------------------------------
# set.seed(31)
# dat <- simulate_smoothbp(
#   n_subj = 25, n_obs = 10,
#   b0 = 5.0, b1 = -0.4, delta = 1.4,
#   omega = 3.2, rho = 4.0,
#   sigma = 0.4, sigma_u = 0.7,
#   seed = 31L
# )
# true_params(dat)

## ----fit-smoothbp-------------------------------------------------------------
# t0 <- Sys.time()
# fit_sbp <- smoothbp(
#   formula = y ~ tau,
#   b0      = ~ 1 + (1 | subject),
#   b1      = ~ 1,
#   deltas  = list(~ 1),
#   omega   = list(~ 1),
#   rho     = list(~ 1),
#   data    = dat,
#   priors  = smoothbp_priors(
#     b0      = prior_normal(0, 10),
#     b1      = prior_normal(0, 2),
#     omega   = list(prior_normal(3, 2, lb = 0, ub = max(dat$tau))),
#     rho     = list(prior_normal(3, 2, lb = 0)),
#     # Diffuse variance priors, matched (as closely as the two parameterisations
#     # allow) to the brms half-t below. smoothbp's variance update is conjugate,
#     # so the family must be inverse-gamma; (0.001, 0.001) makes it effectively
#     # non-informative for the well-identified residual SD.
#     sigma   = prior_invgamma(0.001, 0.001),
#     sigma_u = prior_invgamma(0.001, 0.001)
#   ),
#   chains  = 4L, iter = 2000L, warmup = 1000L,
#   seed    = 31L, .verbose = FALSE
# )
# time_sbp <- as.numeric(difftime(Sys.time(), t0, units = "secs"))

## ----fit-brms, eval = have_brms, message = FALSE, warning = FALSE, results = "hide"----
# suppressPackageStartupMessages(library(brms))
# 
# # Note: brms requires named non-linear parameters in its formula.  Here `b2`
# # is the brms-side name for the slope-change parameter; it corresponds to
# # `delta` (delta_1) in smoothbp's parameterisation.
# bf_smoothed <- brms::bf(
#   y ~ b0 + b1 * (tau - omega) +
#         b2 * (tau - omega) / (1 + exp(-(tau - omega) * rho)),
#   b0 ~ 1 + (1 | subject),
#   b1 ~ 1, b2 ~ 1, omega ~ 1, rho ~ 1,  # b2 = smoothbp's delta_1
#   nl = TRUE
# )
# 
# # brms::prior() captures arguments by non-standard evaluation, so
# # `ub = max(dat$tau)` would be embedded literally in the generated Stan
# # code. Pre-compute the value and pass it via prior_string() which uses
# # standard evaluation and accepts numeric bounds.
# ub_omega <- max(dat$tau)
# 
# priors_brms <- c(
#   brms::prior(normal(0, 10), nlpar = "b0"),
#   brms::prior(normal(0, 2),  nlpar = "b1"),
#   brms::prior(normal(0, 2),  nlpar = "b2"),
#   brms::prior_string("normal(3, 2)", nlpar = "omega",
#                      lb = 0, ub = ub_omega),
#   brms::prior(normal(3, 2),  nlpar = "rho",   lb = 0),
#   # Diffuse variance priors, matched to smoothbp's near-non-informative
#   # inverse-gamma above. A wide half-t leaves the residual SD likelihood-
#   # dominated; the group SD remains weakly identified with only 25 subjects.
#   # In a non-linear model the group SD belongs to its nlpar, so `nlpar = "b0"`
#   # is required (a bare `class = "sd"` matches no parameter).
#   brms::prior(student_t(3, 0, 10), class = "sigma"),
#   brms::prior(student_t(3, 0, 10), class = "sd", nlpar = "b0")
# )
# 
# # Initialise every chain near the prior mean.  Without this, Stan's default
# # uniform-on-(-2, 2) initialisation lands some chains in the unidentifiable
# # omega > max(tau) trap and they never escape.
# init_fun <- function() list(
#   b_b0    = array(rnorm(1, 5, 1)),
#   b_b1    = array(rnorm(1, 0, 0.3)),
#   b_b2    = array(rnorm(1, 0, 0.3)),
#   b_omega = array(rnorm(1, 3, 0.3)),
#   b_rho   = array(rnorm(1, 3, 0.5))
# )
# 
# t0 <- Sys.time()
# fit_brms <- brms::brm(
#   bf_smoothed,
#   data    = dat,
#   prior   = priors_brms,
#   chains  = 4, iter = 2000, warmup = 1000,
#   seed    = 31, refresh = 0,
#   init    = init_fun,
#   control = list(adapt_delta = 0.95)
# )
# time_brms <- as.numeric(difftime(Sys.time(), t0, units = "secs"))

## ----brms-convergence, eval = have_brms---------------------------------------
# sbp_names <- c(
#   "b0_(Intercept)", "b1_(Intercept)", "delta1_(Intercept)",
#   "omega1_(Intercept)", "rho1_(Intercept)",
#   "sigma", "sigma_u"
# )
# brms_names <- c(
#   "b_b0_Intercept", "b_b1_Intercept", "b_b2_Intercept",
#   "b_omega_Intercept", "b_rho_Intercept",
#   "sigma", "sd_subject__b0_Intercept"
# )
# 
# posterior::summarise_draws(
#   posterior::as_draws_df(fit_brms), rhat, ess_bulk
# ) %>%
#   dplyr::filter(variable %in% brms_names) %>%
#   dplyr::transmute(parameter = sbp_names[match(variable, brms_names)],
#                    rhat, ess_bulk) %>%
#   knitr::kable(digits = 3,
#                caption = "brms convergence on population-level parameters.")

## ----compare-summaries, eval = have_brms--------------------------------------
# sbp_draws  <- posterior::as_draws_df(fit_sbp$draws)[, sbp_names]
# brms_draws <- as.data.frame(posterior::as_draws_df(fit_brms))[, brms_names]
# colnames(brms_draws) <- sbp_names
# 
# draws_long <- dplyr::bind_rows(
#   tidyr::pivot_longer(sbp_draws,  cols = dplyr::everything(),
#                       names_to = "parameter", values_to = "value") %>%
#     dplyr::mutate(package = "smoothbp"),
#   tidyr::pivot_longer(brms_draws, cols = dplyr::everything(),
#                       names_to = "parameter", values_to = "value") %>%
#     dplyr::mutate(package = "brms")
# )
# 
# draws_long %>%
#   dplyr::group_by(parameter, package) %>%
#   dplyr::summarise(
#     mean = mean(value),
#     sd   = sd(value),
#     q025 = quantile(value, 0.025),
#     q975 = quantile(value, 0.975),
#     .groups = "drop"
#   ) %>%
#   tidyr::pivot_wider(
#     names_from  = package,
#     values_from = c(mean, sd, q025, q975),
#     names_glue  = "{package}_{.value}"
#   ) %>%
#   dplyr::mutate(
#     delta_mean_pct = 100 * abs(smoothbp_mean - brms_mean) /
#                        (abs(brms_mean) + 1e-8)
#   ) %>%
#   knitr::kable(digits = 3,
#                caption = "Posterior summaries for both packages.")

## ----overlay, eval = have_brms, fig.height = 5--------------------------------
# p_overlay <- ggplot(draws_long, aes(x = value, fill = package, colour = package)) +
#   geom_density(alpha = 0.35) +
#   facet_wrap(~ parameter, scales = "free", ncol = 3) +
#   scale_fill_manual(values  = c(smoothbp = "#1f77b4", brms = "#d62728")) +
#   scale_colour_manual(values = c(smoothbp = "#1f77b4", brms = "#d62728")) +
#   labs(
#     title    = "Marginal posteriors: smoothbp vs brms",
#     subtitle = sprintf(
#       "Same simulated data, same priors. smoothbp: %.1fs   brms: %.1fs",
#       time_sbp, time_brms
#     ),
#     x = NULL, y = "density"
#   ) +
#   theme_minimal(base_size = 11) +
#   theme(legend.position = "bottom")
# p_overlay

## ----w1, eval = have_brms-----------------------------------------------------
# # standardised 1-D 1-Wasserstein distance between two samples
# w1_std <- function(x, y) {
#   p  <- ppoints(1000)
#   qx <- quantile(x, p, names = FALSE, type = 8)
#   qy <- quantile(y, p, names = FALSE, type = 8)
#   mean(abs(qx - qy)) / stats::sd(c(x, y))
# }
# 
# # second smoothbp run (different seed) as a within-sampler Monte Carlo reference
# fit_sbp_b <- smoothbp(
#   formula = y ~ tau,
#   b0      = ~ 1 + (1 | subject),
#   b1      = ~ 1,
#   deltas  = list(~ 1),
#   omega   = list(~ 1),
#   rho     = list(~ 1),
#   data    = dat,
#   priors  = smoothbp_priors(
#     b0      = prior_normal(0, 10),
#     b1      = prior_normal(0, 2),
#     omega   = list(prior_normal(3, 2, lb = 0, ub = max(dat$tau))),
#     rho     = list(prior_normal(3, 2, lb = 0)),
#     sigma   = prior_invgamma(0.001, 0.001),
#     sigma_u = prior_invgamma(0.001, 0.001)
#   ),
#   chains  = 4L, iter = 2000L, warmup = 1000L,
#   seed    = 32L, .verbose = FALSE
# )
# 
# sbp_a  <- posterior::as_draws_df(fit_sbp$draws)[, sbp_names]
# sbp_b  <- posterior::as_draws_df(fit_sbp_b$draws)[, sbp_names]
# brms_w <- as.data.frame(posterior::as_draws_df(fit_brms))[, brms_names]
# colnames(brms_w) <- sbp_names
# 
# data.frame(
#   parameter = sbp_names,
#   W1_cross  = vapply(sbp_names, function(nm) w1_std(sbp_a[[nm]], brms_w[[nm]]), numeric(1)),
#   W1_null   = vapply(sbp_names, function(nm) w1_std(sbp_a[[nm]], sbp_b[[nm]]),  numeric(1)),
#   row.names = NULL
# ) |>
#   knitr::kable(digits = 3,
#     caption = "Standardised 1-Wasserstein distance (pooled-SD units): smoothbp vs
#                brms (W1_cross) and a within-sampler reference from two smoothbp
#                seeds (W1_null). Cross comparable to null => indistinguishable up
#                to Monte Carlo error.")

## ----overlay-export, eval = FALSE---------------------------------------------
# ggplot2::ggsave("figure_overlay_s1.png", p_overlay,
#                 width = 7, height = 5, dpi = 150)

## ----ess, eval = have_brms----------------------------------------------------
# ess_sbp  <- posterior::summarise_draws(fit_sbp$draws, ess_bulk) %>%
#   dplyr::filter(variable %in% sbp_names) %>%
#   dplyr::transmute(parameter = variable,
#                    smoothbp_ess_per_sec = ess_bulk / time_sbp)
# 
# ess_brms <- posterior::summarise_draws(
#   posterior::as_draws_df(fit_brms), ess_bulk
# ) %>%
#   dplyr::filter(variable %in% brms_names) %>%
#   dplyr::mutate(parameter = sbp_names[match(variable, brms_names)]) %>%
#   dplyr::transmute(parameter, brms_ess_per_sec = ess_bulk / time_brms)
# 
# dplyr::full_join(ess_sbp, ess_brms, by = "parameter") %>%
#   knitr::kable(digits = 1,
#                caption = "Bulk ESS / second by package and parameter.")

## ----timing-s1, eval = have_brms----------------------------------------------
# data.frame(
#   package     = c("smoothbp", "brms (Stan/NUTS)"),
#   wall_clock_s = round(c(time_sbp, time_brms), 1),
#   speedup      = c(sprintf("%.0fx faster", time_brms / time_sbp), "reference")
# ) |>
#   knitr::kable(col.names = c("Package", "Wall-clock (s)", "Relative speed"),
#                caption   = "Scenario 1 wall-clock time. Single run on the same machine.")

## ----simulate-cov-------------------------------------------------------------
# set.seed(8147)
# 
# n_subj   <- 30
# n_obs    <- 10
# tau_range <- c(0, 6)
# 
# b0_true      <- 5.0
# b1_true      <- -0.4
# delta_true   <- 1.4
# omega_int    <- 3.2
# omega_trt    <- -0.8
# rho_true     <- 4.0
# sigma_true   <- 0.4
# sigma_u_true <- 0.7
# 
# treatment <- rep(c(0, 1), each = n_subj / 2)
# u_j <- rnorm(n_subj, 0, sigma_u_true)
# .sigmoid <- function(x) 1 / (1 + exp(-x))
# 
# rows <- vector("list", n_subj)
# for (j in seq_len(n_subj)) {
#   tau_j   <- seq(tau_range[1], tau_range[2], length.out = n_obs)
#   omega_j <- omega_int + omega_trt * treatment[j]
#   d_j     <- tau_j - omega_j
#   s_j     <- .sigmoid(d_j * rho_true)
#   mu_j    <- (b0_true + u_j[j]) + b1_true * d_j + delta_true * d_j * s_j
#   rows[[j]] <- data.frame(
#     subject   = factor(j),
#     tau       = tau_j,
#     treatment = treatment[j],
#     y         = mu_j + rnorm(n_obs, 0, sigma_true)
#   )
# }
# dat_cov <- do.call(rbind, rows)
# 
# cat(sprintf("True omega: intercept = %.1f, treatment effect = %.1f\n",
#             omega_int, omega_trt))
# cat(sprintf("  Control  (treatment = 0): omega = %.1f\n", omega_int))
# cat(sprintf("  Treated  (treatment = 1): omega = %.1f\n", omega_int + omega_trt))

## ----fit-smoothbp-cov---------------------------------------------------------
# t0 <- Sys.time()
# fit_sbp_cov <- smoothbp(
#   formula = y ~ tau,
#   b0      = ~ 1 + (1 | subject),
#   b1      = ~ 1,
#   deltas  = list(~ 1),
#   omega   = list(~ 1 + treatment),
#   rho     = list(~ 1),
#   data    = dat_cov,
#   priors  = smoothbp_priors(
#     b0    = prior_normal(0, 10),
#     b1    = prior_normal(0, 2),
#     omega = list(list(
#       "(Intercept)" = prior_normal(3, 2, lb = 0, ub = max(dat_cov$tau)),
#       "treatment"   = prior_normal(0, 2)
#     )),
#     rho   = list(prior_normal(3, 2, lb = 0))
#   ),
#   chains  = 4L, iter = 2000L, warmup = 1000L,
#   seed    = 8147L, .verbose = FALSE
# )
# time_sbp_cov <- as.numeric(difftime(Sys.time(), t0, units = "secs"))

## ----fit-brms-cov, eval = have_brms, message = FALSE, warning = FALSE, results = "hide"----
# bf_cov <- brms::bf(
#   y ~ b0 + b1 * (tau - omega) +
#         b2 * (tau - omega) / (1 + exp(-(tau - omega) * rho)),
#   b0    ~ 1 + (1 | subject),
#   b1    ~ 1,
#   b2    ~ 1,
#   omega ~ 1 + treatment,
#   rho   ~ 1,
#   nl = TRUE
# )
# 
# priors_brms_cov <- c(
#   brms::prior(normal(0, 10), nlpar = "b0"),
#   brms::prior(normal(0, 2),  nlpar = "b1"),
#   brms::prior(normal(0, 2),  nlpar = "b2"),
#   brms::prior_string("normal(3, 2)", nlpar = "omega", coef = "Intercept"),
#   brms::prior(normal(0, 2),  nlpar = "omega", coef = "treatment"),
#   brms::prior(normal(3, 2),  nlpar = "rho", lb = 0)
# )
# 
# init_fun_cov <- function() list(
#   b_b0    = array(rnorm(1, 5, 1)),
#   b_b1    = array(rnorm(1, 0, 0.3)),
#   b_b2    = array(rnorm(1, 0, 0.3)),
#   b_omega = array(c(rnorm(1, 3, 0.3), rnorm(1, -0.5, 0.3))),
#   b_rho   = array(rnorm(1, 3, 0.5))
# )
# 
# t0 <- Sys.time()
# fit_brms_cov <- brms::brm(
#   bf_cov,
#   data    = dat_cov,
#   prior   = priors_brms_cov,
#   chains  = 4, iter = 2000, warmup = 1000,
#   seed    = 8147, refresh = 0,
#   init    = init_fun_cov,
#   control = list(adapt_delta = 0.95)
# )
# time_brms_cov <- as.numeric(difftime(Sys.time(), t0, units = "secs"))

## ----compare-cov, eval = have_brms--------------------------------------------
# sbp_names_cov <- c(
#   "b0_(Intercept)", "b1_(Intercept)", "delta1_(Intercept)",
#   "omega1_(Intercept)", "omega1_treatment",
#   "rho1_(Intercept)", "sigma", "sigma_u"
# )
# brms_names_cov <- c(
#   "b_b0_Intercept", "b_b1_Intercept", "b_b2_Intercept",
#   "b_omega_Intercept", "b_omega_treatment",
#   "b_rho_Intercept", "sigma", "sd_subject__b0_Intercept"
# )
# 
# sbp_draws_cov  <- posterior::as_draws_df(fit_sbp_cov$draws)[, sbp_names_cov]
# brms_draws_cov <- as.data.frame(
#   posterior::as_draws_df(fit_brms_cov)
# )[, brms_names_cov]
# colnames(brms_draws_cov) <- sbp_names_cov
# 
# draws_long_cov <- dplyr::bind_rows(
#   tidyr::pivot_longer(sbp_draws_cov, cols = dplyr::everything(),
#                       names_to = "parameter", values_to = "value") %>%
#     dplyr::mutate(package = "smoothbp"),
#   tidyr::pivot_longer(brms_draws_cov, cols = dplyr::everything(),
#                       names_to = "parameter", values_to = "value") %>%
#     dplyr::mutate(package = "brms")
# )
# 
# draws_long_cov %>%
#   dplyr::group_by(parameter, package) %>%
#   dplyr::summarise(
#     mean = mean(value),
#     sd   = sd(value),
#     q025 = quantile(value, 0.025),
#     q975 = quantile(value, 0.975),
#     .groups = "drop"
#   ) %>%
#   tidyr::pivot_wider(
#     names_from  = package,
#     values_from = c(mean, sd, q025, q975),
#     names_glue  = "{package}_{.value}"
#   ) %>%
#   dplyr::mutate(
#     delta_mean_pct = 100 * abs(smoothbp_mean - brms_mean) /
#                        (abs(brms_mean) + 1e-8)
#   ) %>%
#   knitr::kable(digits = 3,
#                caption = "Posterior summaries: omega ~ 1 + treatment.")

## ----overlay-cov, eval = have_brms, fig.height = 5----------------------------
# true_vals_cov <- data.frame(
#   parameter = c("b0_(Intercept)", "b1_(Intercept)", "delta1_(Intercept)",
#                  "omega1_(Intercept)", "omega1_treatment",
#                  "rho1_(Intercept)", "sigma", "sigma_u"),
#   true_value = c(b0_true, b1_true, delta_true,
#                   omega_int, omega_trt, rho_true, sigma_true, sigma_u_true)
# )
# 
# ggplot(draws_long_cov, aes(x = value, fill = package, colour = package)) +
#   geom_density(alpha = 0.35) +
#   geom_vline(data = true_vals_cov, aes(xintercept = true_value),
#              linetype = "dashed", linewidth = 0.5) +
#   facet_wrap(~ parameter, scales = "free", ncol = 3) +
#   scale_fill_manual(values  = c(smoothbp = "#1f77b4", brms = "#d62728")) +
#   scale_colour_manual(values = c(smoothbp = "#1f77b4", brms = "#d62728")) +
#   labs(
#     title    = "Marginal posteriors: smoothbp vs brms (omega ~ 1 + treatment)",
#     subtitle = sprintf(
#       "smoothbp: %.1fs   brms: %.1fs.  Dashed lines = true values.",
#       time_sbp_cov, time_brms_cov
#     ),
#     x = NULL, y = "density"
#   ) +
#   theme_minimal(base_size = 11) +
#   theme(legend.position = "bottom")

## ----ess-cov, eval = have_brms------------------------------------------------
# ess_sbp_cov  <- posterior::summarise_draws(fit_sbp_cov$draws, ess_bulk) %>%
#   dplyr::filter(variable %in% sbp_names_cov) %>%
#   dplyr::transmute(parameter = variable,
#                    smoothbp_ess_per_sec = ess_bulk / time_sbp_cov)
# 
# ess_brms_cov <- posterior::summarise_draws(
#   posterior::as_draws_df(fit_brms_cov), ess_bulk
# ) %>%
#   dplyr::filter(variable %in% brms_names_cov) %>%
#   dplyr::mutate(parameter = sbp_names_cov[match(variable, brms_names_cov)]) %>%
#   dplyr::transmute(parameter, brms_ess_per_sec = ess_bulk / time_brms_cov)
# 
# dplyr::full_join(ess_sbp_cov, ess_brms_cov, by = "parameter") %>%
#   knitr::kable(digits = 1,
#                caption = "Bulk ESS / second by package and parameter (covariate model).")

## ----timing-s2, eval = have_brms----------------------------------------------
# data.frame(
#   package      = c("smoothbp", "brms (Stan/NUTS)"),
#   wall_clock_s = round(c(time_sbp_cov, time_brms_cov), 1),
#   speedup       = c(sprintf("%.0fx faster", time_brms_cov / time_sbp_cov), "reference")
# ) |>
#   knitr::kable(col.names = c("Package", "Wall-clock (s)", "Relative speed"),
#                caption   = "Scenario 2 wall-clock time. Covariate-on-omega model.")

## ----mcp-simulate-------------------------------------------------------------
# set.seed(5501)
# 
# dat_mcp <- simulate_smoothbp(
#   n_subj    = 1,
#   n_obs     = 150,
#   b0        = 5.0,
#   b1        = -0.4,
#   delta     =  1.4,
#   omega     =  3.2,
#   rho       =  5.0,   # moderately sharp — should agree well with mcp
#   sigma     =  0.4,
#   sigma_u   =  0.0,   # no random effects
#   tau_range = c(0, 6),
#   seed      = 5501L
# )
# 
# # mcp expects a plain data.frame without the factor column
# dat_mcp_plain <- data.frame(tau = dat_mcp$tau, y = dat_mcp$y)
# 
# cat(sprintf("n = %d  |  true omega = 3.2  |  true b1 = -0.4  |  true delta1 = 1.4\n",
#             nrow(dat_mcp_plain)))

## ----mcp-smoothbp-------------------------------------------------------------
# t0 <- Sys.time()
# fit_sbp_mcp <- smoothbp(
#   formula = y ~ tau,
#   b0      = ~ 1,
#   b1      = ~ 1,
#   deltas  = list(~ 1),
#   omega   = list(~ 1),
#   rho     = list(fixed(100)), # Use fixed() to specify a sharp, hard kink
#   data    = dat_mcp,
#   priors  = smoothbp_priors(
#     b0    = prior_normal(0, 10),
#     b1    = prior_normal(0, 2),
#     omega = list(prior_normal(3, 2, lb = 0, ub = max(dat_mcp$tau)))
#   ),
#   chains = 4L, iter = 2000L, warmup = 1000L,
#   seed = 5501L, .verbose = FALSE
# )
# time_sbp_mcp <- as.numeric(difftime(Sys.time(), t0, units = "secs"))
# cat(sprintf("smoothbp: %.1f s\n", time_sbp_mcp))

## ----mcp-fit, eval = have_mcp, message = FALSE, warning = FALSE, results = "hide"----
# # Two-segment joined-slope model: the natural mcp equivalent
# model_mcp <- list(
#   y ~ 1 + tau,   # segment 1: intercept + pre-change slope
#   ~ 0 + tau      # segment 2: joined at cp_1, post-change slope
# )
# 
# # Restrict cp_1 to the observed time range, matching the smoothbp bound
# prior_mcp <- list(cp_1 = "dunif(0, 6)")
# set.seed(5501)
# t0 <- Sys.time()
# fit_mcp <- mcp::mcp(
#   model_mcp,
#   data  = dat_mcp_plain,
#   prior = prior_mcp,
#   par_x = "tau",
#   chains = 4L, iter = 2000L, adapt = 1000L,
# )
# time_mcp <- as.numeric(difftime(Sys.time(), t0, units = "secs"))
# cat(sprintf("mcp: %.1f s\n", time_mcp))

## ----mcp-compare, eval = have_mcp---------------------------------------------
# sbp_sum_mcp <- posterior::summarise_draws(fit_sbp_mcp$draws) |>
#   dplyr::filter(variable %in% c("b1_(Intercept)", "delta1_(Intercept)", "omega1_(Intercept)",
#                                "rho1_(Intercept)", "sigma"))
# 
# mcp_sum <- summary(fit_mcp) |>
#   dplyr::filter(name %in% c("cp_1", "int_1", "tau_1", "tau_2", "sigma_1"))
# 
# # Compute post-change slope from smoothbp draws for direct comparison
# dm_sbp_mcp <- posterior::as_draws_matrix(fit_sbp_mcp$draws)
# post_slope <- as.numeric(dm_sbp_mcp[, "b1_(Intercept)"] +
#                            dm_sbp_mcp[, "delta1_(Intercept)"])
# 
# comparison_mcp <- data.frame(
#   parameter          = c("omega / cp_1",
#                           "b1 / time_1 (pre-change slope)",
#                           "b1+delta1 / time_2 (post-change slope)",
#                           "sigma"),
#   truth              = c(3.2, -0.4, -0.4 + 1.4, 0.4),
#   sbp_mean           = c(
#     sbp_sum_mcp$mean[sbp_sum_mcp$variable == "omega1_(Intercept)"],
#     sbp_sum_mcp$mean[sbp_sum_mcp$variable == "b1_(Intercept)"],
#     mean(post_slope),
#     sbp_sum_mcp$mean[sbp_sum_mcp$variable == "sigma"]
#   ),
#   sbp_lo             = c(
#     sbp_sum_mcp$q2.5[sbp_sum_mcp$variable == "omega1_(Intercept)"],
#     sbp_sum_mcp$q2.5[sbp_sum_mcp$variable == "b1_(Intercept)"],
#     quantile(post_slope, 0.025),
#     sbp_sum_mcp$q2.5[sbp_sum_mcp$variable == "sigma"]
#   ),
#   sbp_hi             = c(
#     sbp_sum_mcp$q97.5[sbp_sum_mcp$variable == "omega1_(Intercept)"],
#     sbp_sum_mcp$q97.5[sbp_sum_mcp$variable == "b1_(Intercept)"],
#     quantile(post_slope, 0.975),
#     sbp_sum_mcp$q97.5[sbp_sum_mcp$variable == "sigma"]
#   ),
#   mcp_mean           = mcp_sum$mean[match(c("cp_1","tau_1","tau_2","sigma_1"),
#                                           mcp_sum$name)],
#   mcp_lo             = mcp_sum$lower[match(c("cp_1","tau_1","tau_2","sigma_1"),
#                                            mcp_sum$name)],
#   mcp_hi             = mcp_sum$upper[match(c("cp_1","tau_1","tau_2","sigma_1"),
#                                            mcp_sum$name)]
# )
# 
# knitr::kable(comparison_mcp, digits = 3,
#   caption = "Scenario 3: smoothbp vs mcp on comparable parameters.
#              smoothbp post-change slope is b1 + delta1.")

## ----mcp-trajectories, eval = have_mcp, fig.height = 4.5----------------------
# tau_grid <- seq(0, 6, length.out = 200)
# 
# # smoothbp posterior mean trajectory
# draws_sbp_mcp  <- posterior::as_draws_matrix(fit_sbp_mcp$draws)
# b0_m   <- mean(draws_sbp_mcp[, "b0_(Intercept)"])
# b1_m   <- mean(draws_sbp_mcp[, "b1_(Intercept)"])
# delta1_m <- mean(draws_sbp_mcp[, "delta1_(Intercept)"])
# om_m   <- mean(draws_sbp_mcp[, "omega1_(Intercept)"])
# rho_m  <- mean(draws_sbp_mcp[, "rho1_(Intercept)"])
# 
# sigmoid  <- function(x) 1 / (1 + exp(-x))
# d_grid   <- tau_grid - om_m
# mu_sbp   <- b0_m + b1_m * d_grid + delta1_m * d_grid * sigmoid(d_grid * rho_m)
# 
# # mcp posterior mean trajectory
# cp_m    <- mcp_sum$mean[mcp_sum$name == "cp_1"]
# int_m   <- mcp_sum$mean[mcp_sum$name == "int_1"]
# t1_m    <- mcp_sum$mean[mcp_sum$name == "tau_1"]
# t2_m    <- mcp_sum$mean[mcp_sum$name == "tau_2"]
# mu_mcp  <- ifelse(tau_grid < cp_m,
#                   int_m + t1_m * tau_grid,
#                   int_m + t1_m * cp_m + t2_m * (tau_grid - cp_m))
# 
# traj_df <- data.frame(
#   tau    = rep(tau_grid, 2),
#   mu     = c(mu_sbp, mu_mcp),
#   model  = rep(c("smoothbp (logistic)", "mcp (hard kink)"), each = 200)
# )
# 
# ggplot() +
#   geom_point(aes(tau, y), data = dat_mcp_plain, alpha = 0.25, size = 0.8) +
#   geom_line(aes(tau, mu, colour = model, linetype = model),
#             data = traj_df, linewidth = 0.9) +
#   geom_vline(xintercept = 3.2, linetype = "dotted", colour = "grey40") +
#   scale_colour_manual(values = c("smoothbp (logistic)" = "#2166ac",
#                                  "mcp (hard kink)"     = "#d6604d")) +
#   labs(title    = "Posterior mean fitted trajectories",
#        subtitle = "Dotted line = true change-point (omega = 3.2)",
#        x = "tau", y = "y", colour = NULL, linetype = NULL) +
#   theme_bw() +
#   theme(legend.position = "bottom")

## ----mcp-ess, eval = have_mcp-------------------------------------------------
# # smoothbp ESS on omega
# sbp_ess_mcp <- posterior::summarise_draws(
#   fit_sbp_mcp$draws[,, c("omega1_(Intercept)", "b1_(Intercept)",
#                           "delta1_(Intercept)", "sigma")],
#   ess_bulk = posterior::ess_bulk
# ) |>
#   dplyr::mutate(ess_per_sec = ess_bulk / time_sbp_mcp,
#                 package = "smoothbp")
# 
# # mcp ESS — extract draws from mcp object
# # mcp stores posterior draws as a coda::mcmc.list in $mcmc_post
# mcp_draws <- posterior::as_draws_array(fit_mcp$mcmc_post)
# mcp_ess_draws <- posterior::summarise_draws(
#   mcp_draws[, , c("cp_1", "tau_1", "sigma_1")],
#   ess_bulk = posterior::ess_bulk
# )
# mcp_ess_raw <- setNames(mcp_ess_draws$ess_bulk, mcp_ess_draws$variable)
# 
# data.frame(
#   parameter   = c("omega / cp_1", "b1 / time_1", "sigma"),
#   sbp_ess_s   = round(sbp_ess_mcp$ess_per_sec[
#     match(c("omega1_(Intercept)", "b1_(Intercept)", "sigma"),
#           sbp_ess_mcp$variable)], 1),
#   mcp_ess_s   = round(c(mcp_ess_raw["cp_1"],
#                          mcp_ess_raw["tau_1"],
#                          mcp_ess_raw["sigma_1"]) / time_mcp, 1)
# ) |>
#   knitr::kable(caption = "ESS/second: smoothbp vs mcp.
#                           Higher is better.")

## ----timing-s3----------------------------------------------------------------
# data.frame(
#   package      = c("smoothbp", "mcp (JAGS)"),
#   wall_clock_s = round(c(time_sbp_mcp, time_mcp), 1),
#   speedup       = c(sprintf("%.0fx faster", time_mcp / time_sbp_mcp), "reference")
# ) |>
#   knitr::kable(col.names = c("Package", "Wall-clock (s)", "Relative speed"),
#                caption   = "Scenario 3 wall-clock time. Single-group, n = 150.")

## ----simulate-multi-----------------------------------------------------------
# set.seed(9801)
# n  <- 300
# tau <- seq(0, 12, length.out = n)
# 
# # True model: two breakpoints
# om1_true    <-  3.5;  rho1_true <- 5.0;  delta1_true <- -1.2
# om2_true    <-  8.0;  rho2_true <- 4.0;  delta2_true <-  1.8
# b0_true     <-  5.0;  b1_true   <-  0.3;  sigma_true  <-  0.3
# 
# d1  <- tau - om1_true;  s1 <- plogis(d1 * rho1_true)
# d2  <- tau - om2_true;  s2 <- plogis(d2 * rho2_true)
# mu  <- b0_true + b1_true * (tau - om1_true) +
#        delta1_true * d1 * s1 + delta2_true * d2 * s2
# y   <- mu + rnorm(n, sd = sigma_true)
# dat_multi <- data.frame(y = y, tau = tau)
# 
# cat(sprintf("True: b0=%.1f b1=%.1f delta1=%.1f omega1=%.1f delta2=%.1f omega2=%.1f sigma=%.1f\n",
#             b0_true, b1_true, delta1_true, om1_true, delta2_true, om2_true, sigma_true))

## ----fit-smoothbp-multi-------------------------------------------------------
# t0 <- Sys.time()
# fit_sbp_multi <- smoothbp(
#   formula = y ~ tau,
#   b0      = ~ 1,
#   b1      = ~ 1,
#   deltas  = list(~ 1, ~ 1),
#   omega   = list(~ 1, ~ 1),
#   rho     = list(~ 1, ~ 1),
#   data    = dat_multi,
#   priors  = smoothbp_priors(
#     b0     = prior_normal(0, 10),
#     b1     = prior_normal(0, 2),
#     omega  = list(
#       prior_normal(3.5, 2, lb = 0, ub = 6),
#       prior_normal(8.0, 2, lb = 4, ub = max(dat_multi$tau))
#     ),
#     rho    = list(prior_normal(4, 2, lb = 0),
#                   prior_normal(4, 2, lb = 0))
#   ),
#   chains = 4L, iter = 2000L, warmup = 1000L,
#   seed = 9801L, .verbose = FALSE
# )
# time_sbp_multi <- as.numeric(difftime(Sys.time(), t0, units = "secs"))
# cat(sprintf("smoothbp: %.1f s\n", time_sbp_multi))

## ----fit-brms-multi, eval = have_brms, message = FALSE, warning = FALSE, results = "hide"----
# suppressPackageStartupMessages(library(brms))
# 
# # Two-breakpoint non-linear formula for brms
# bf_multi <- brms::bf(
#   y ~ b0 + b1 * (tau - om1) +
#         d1 * (tau - om1) / (1 + exp(-(tau - om1) * rho1)) +
#         d2 * (tau - om2) / (1 + exp(-(tau - om2) * rho2)),
#   b0 ~ 1, b1 ~ 1,
#   d1 ~ 1, om1 ~ 1, rho1 ~ 1,
#   d2 ~ 1, om2 ~ 1, rho2 ~ 1,
#   nl = TRUE
# )
# 
# priors_brms_multi <- c(
#   brms::prior(normal(0, 10),  nlpar = "b0"),
#   brms::prior(normal(0, 2),   nlpar = "b1"),
#   brms::prior(normal(0, 2),   nlpar = "d1"),
#   brms::prior_string("normal(3.5, 2)", nlpar = "om1", lb = 0, ub = 6),
#   brms::prior(normal(4, 2),   nlpar = "rho1", lb = 0),
#   brms::prior(normal(0, 2),   nlpar = "d2"),
#   brms::prior_string("normal(8.0, 2)", nlpar = "om2", lb = 4, ub = 12),
#   brms::prior(normal(4, 2),   nlpar = "rho2", lb = 0)
# )
# 
# init_multi <- function() list(
#   b_b0   = array(rnorm(1, 5, 0.5)),
#   b_b1   = array(rnorm(1, 0.3, 0.1)),
#   b_d1   = array(rnorm(1, -1, 0.3)),
#   b_om1  = array(rnorm(1, 3.5, 0.3)),
#   b_rho1 = array(rnorm(1, 5, 0.5)),
#   b_d2   = array(rnorm(1, 1.5, 0.3)),
#   b_om2  = array(rnorm(1, 8, 0.3)),
#   b_rho2 = array(rnorm(1, 4, 0.5))
# )
# 
# t0 <- Sys.time()
# fit_brms_multi <- brms::brm(
#   bf_multi,
#   data    = dat_multi,
#   prior   = priors_brms_multi,
#   chains  = 4, iter = 2000, warmup = 1000,
#   seed    = 9801, refresh = 0,
#   init    = init_multi,
#   control = list(adapt_delta = 0.95)
# )
# time_brms_multi <- as.numeric(difftime(Sys.time(), t0, units = "secs"))

## ----compare-multi, eval = have_brms------------------------------------------
# sbp_names_multi <- c(
#   "b0_(Intercept)", "b1_(Intercept)",
#   "delta1_(Intercept)", "omega1_(Intercept)", "rho1_(Intercept)",
#   "delta2_(Intercept)", "omega2_(Intercept)", "rho2_(Intercept)",
#   "sigma"
# )
# brms_names_multi <- c(
#   "b_b0_Intercept", "b_b1_Intercept",
#   "b_d1_Intercept", "b_om1_Intercept", "b_rho1_Intercept",
#   "b_d2_Intercept", "b_om2_Intercept", "b_rho2_Intercept",
#   "sigma"
# )
# true_vals_multi <- c(b0_true, b1_true,
#                      delta1_true, om1_true, rho1_true,
#                      delta2_true, om2_true, rho2_true,
#                      sigma_true)
# 
# sbp_draws_multi  <- posterior::as_draws_df(fit_sbp_multi$draws)[, sbp_names_multi]
# brms_draws_multi <- as.data.frame(
#   posterior::as_draws_df(fit_brms_multi)
# )[, brms_names_multi]
# colnames(brms_draws_multi) <- sbp_names_multi
# 
# draws_long_multi <- dplyr::bind_rows(
#   tidyr::pivot_longer(sbp_draws_multi, cols = dplyr::everything(),
#                       names_to = "parameter", values_to = "value") |>
#     dplyr::mutate(package = "smoothbp"),
#   tidyr::pivot_longer(brms_draws_multi, cols = dplyr::everything(),
#                       names_to = "parameter", values_to = "value") |>
#     dplyr::mutate(package = "brms")
# )
# 
# draws_long_multi |>
#   dplyr::group_by(parameter, package) |>
#   dplyr::summarise(mean = mean(value), sd = sd(value),
#                    q025 = quantile(value, 0.025),
#                    q975 = quantile(value, 0.975), .groups = "drop") |>
#   tidyr::pivot_wider(names_from = package,
#                      values_from = c(mean, sd, q025, q975),
#                      names_glue = "{package}_{.value}") |>
#   dplyr::mutate(
#     truth          = true_vals_multi[match(parameter, sbp_names_multi)],
#     delta_mean_pct = 100 * abs(smoothbp_mean - brms_mean) /
#                        (abs(brms_mean) + 1e-8)
#   ) |>
#   knitr::kable(digits = 3,
#                caption = "Scenario 4: two-breakpoint smoothbp vs brms.")

## ----overlay-multi, eval = have_brms, fig.height = 6--------------------------
# true_df_multi <- data.frame(
#   parameter  = sbp_names_multi,
#   true_value = true_vals_multi
# )
# 
# ggplot(draws_long_multi, aes(x = value, fill = package, colour = package)) +
#   geom_density(alpha = 0.35, linewidth = 0.6) +
#   geom_vline(data = true_df_multi,
#              aes(xintercept = true_value),
#              linetype = "dashed", linewidth = 0.5, colour = "black") +
#   facet_wrap(~ parameter, scales = "free", ncol = 3) +
#   scale_fill_manual(values   = c(smoothbp = "#1f77b4", brms = "#d62728")) +
#   scale_colour_manual(values = c(smoothbp = "#1f77b4", brms = "#d62728")) +
#   labs(
#     title    = "Scenario 4: Multi-breakpoint posteriors — smoothbp vs brms",
#     subtitle = sprintf(
#       "Two true change-points. smoothbp: %.1fs   brms: %.1fs\nDashed = true value. Overlapping densities = agreement.",
#       time_sbp_multi, time_brms_multi
#     ),
#     x = NULL, y = "Density"
#   ) +
#   theme_minimal(base_size = 11) +
#   theme(legend.position = "bottom")

## ----ess-multi, eval = have_brms----------------------------------------------
# ess_sbp_multi <- posterior::summarise_draws(
#   fit_sbp_multi$draws[, , sbp_names_multi], ess_bulk = posterior::ess_bulk
# ) |>
#   dplyr::transmute(parameter = variable,
#                    smoothbp_ess_per_sec = ess_bulk / time_sbp_multi)
# 
# ess_brms_multi <- posterior::summarise_draws(
#   posterior::as_draws_df(fit_brms_multi), ess_bulk = posterior::ess_bulk
# ) |>
#   dplyr::filter(variable %in% brms_names_multi) |>
#   dplyr::mutate(parameter = sbp_names_multi[match(variable, brms_names_multi)]) |>
#   dplyr::transmute(parameter, brms_ess_per_sec = ess_bulk / time_brms_multi)
# 
# dplyr::full_join(ess_sbp_multi, ess_brms_multi, by = "parameter") |>
#   knitr::kable(digits = 1,
#                caption = "ESS/second: smoothbp vs brms (two-breakpoint model).")

## ----timing-s4, eval = have_brms----------------------------------------------
# data.frame(
#   package      = c("smoothbp", "brms (Stan/NUTS)"),
#   wall_clock_s = round(c(time_sbp_multi, time_brms_multi), 1),
#   speedup       = c(sprintf("%.0fx faster", time_brms_multi / time_sbp_multi), "reference")
# ) |>
#   knitr::kable(col.names = c("Package", "Wall-clock (s)", "Relative speed"),
#                caption   = "Scenario 4 wall-clock time. Two-breakpoint model (single group, fixed effects).")