--- title: "Bayesian regression tables" description: > Publication-ready tables for Bayesian regression in R with `table_regression()`: posterior medians and credible intervals, the deliberate absence of p-values, convergence checks before reading, prior disclosure, random effects from the posterior draws, and the `rstanarm` and `brms` engines. output: rmarkdown::html_vignette: toc: true toc_depth: 2 vignette: > %\VignetteIndexEntry{Bayesian regression tables} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} have_rstanarm <- requireNamespace("rstanarm", quietly = TRUE) && requireNamespace("posterior", quietly = TRUE) knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = have_rstanarm ) build_rich_tables <- identical(Sys.getenv("IN_PKGDOWN"), "true") # The brms chunk compiles C++ (minutes, and a working toolchain) -- # too fragile for CI site builds. Opt in locally with # Sys.setenv(SPICY_BUILD_BRMS = "true") before knitting. run_brms <- identical(Sys.getenv("SPICY_BUILD_BRMS"), "true") && requireNamespace("brms", quietly = TRUE) source("_pkgdown-helpers.R") ``` ```{r setup, message = FALSE} library(spicy) library(rstanarm) ``` This vignette covers **Bayesian regression tables** — the same `table_regression()` call on a posterior instead of a maximum-likelihood fit. The companion vignette [*Publication-ready regression tables*](table-regression.html) covers the shared mechanics (the class-by-class map is [*Supported models*](table-regression-supported-models.html)); here we focus on what genuinely changes when the estimates are draws from a posterior distribution: what the columns mean, which familiar columns are *deliberately absent*, and what to check before reading anything. `table_regression()` supports two engines: * **`rstanarm`** (`stan_glm()`, `stan_lmer()`, `stan_glmer()`) — ships precompiled models, so fits run in seconds and every chunk in this vignette executes live; * **`brms`** (`brm()`) — a wider model space, at the cost of compiling each model; shown at the end (the code is identical in spirit). The running example is the logistic regression of `smoking` from the main vignette's glm section — same data (`sochealth`), same predictors — so the Bayesian table can be read against its frequentist twin. ```{r fit} sh <- sochealth sh$education <- factor(sh$education, ordered = FALSE) fit <- stan_glm(smoking ~ sex + age + education, data = sh, family = binomial(), iter = 1000, chains = 2, seed = 123, refresh = 0) ``` The sampler settings are trimmed so this vignette builds in seconds — for real work, keep rstanarm's defaults (four chains of 2,000 iterations): R-hat compares chains, so it needs several to compare, and credible bounds are tail quantiles, estimated exactly where draws are scarcest (Vehtari et al. 2021). ## What changes in a Bayesian table ```{r basic} table_regression(fit) ``` Reading the columns against their frequentist counterparts: * **B is the posterior median** — the 50% quantile of each coefficient's draws, robust to skewed posteriors (the `rstanarm` convention). It is not a maximum-likelihood point estimate, although under the default weakly informative priors and this sample size the two are nearly identical. * **SE is the posterior MAD SD** — the median absolute deviation of the draws, rescaled so it equals the standard deviation for a normal posterior — and the footer says so (`Std. errors: posterior MAD SD (scaled median absolute deviation).`) rather than borrowing the frequentist term. Median + MAD SD is the robust pairing of *Regression and Other Stories* (§5.3) and of `rstanarm`'s own print: both members resist skewed posteriors, where a mean + SD pair can be dominated by one heavy tail. As Gelman, Hill & Vehtari put it, "for convenience we sometimes call the posterior mad sd the 'standard error'". * **The interval is a credible interval**, and the header says so: `95% CrI`, the equal-tailed 2.5% and 97.5% posterior quantiles. Its reading is the direct probability statement that confidence intervals are so often misread as making: *given model, priors, and data, the coefficient lies in [-1.31, -0.56] with 95% probability* (the `Tertiary` row). No appeal to repeated sampling is required. The equal-tailed interval is the default because it is transformation-invariant (which is what makes the exact odds-ratio bounds below possible) and matches `rstanarm`; the Kruschke tradition's **highest-density interval** is one argument away — `ci_method = "hdi"` computes the shortest interval containing 95% of the draws (Kruschke 2015) and relabels the header `95% HDI`. For these near-symmetric posteriors the two differ only in the second decimal; they part ways on skewed posteriors, where the HDI shifts toward the mode. * **There is no p column, by design.** A posterior has no p-value, and the Bayesian Analysis Reporting Guidelines (Kruschke 2021) ask that any decision rule be stated explicitly and kept separate from the posterior summary itself — so the default table reports the posterior and leaves decisions to the reader (requesting `"p"` explicitly is refused; in a mixed frequentist–Bayesian table the shared p column simply dashes the Bayesian rows). The closest posterior summary — the **probability of direction** (`pd`, the share of draws on the dominant side of zero; Makowski et al. 2019) — is opt-in: ```{r pd} table_regression(fit, show_columns = c("b", "ci", "pd")) ``` Read against the table: `Male`, at .617, is barely better than a coin flip on sign; `age`, at .886, leans positive without settling the direction; the education contrasts show 1.000 and .999 — essentially every retained draw is negative (all 1,000 for `Tertiary`; 999 of 1,000 for `Upper secondary` — the three-decimal display keeps the two cases distinguishable). A `pd` of .500 would say the draws split evenly, with zero at the posterior median. It reads as evidence of *direction*, never of practical size — and one disclosure keeps this section honest: under weak priors `pd` maps almost one-to-one onto the frequentist one-sided p-value (two-sided p ≈ 2(1 − pd), so `age`'s 0.89 corresponds to p ≈ 0.23; Makowski et al. 2019). `pd` is not a route around the p-value's pathologies; it is the same directional information stated as a posterior probability. For a claim of practical size, the Kruschke tradition asks how much of the posterior falls inside a **region of practical equivalence** (ROPE) around zero. spicy deliberately renders no ROPE column: the reporting guidelines require the equivalence limits to be *justified* (BARG step 4.C), and limits are substantive, per-coefficient claims — ±0.05 log-odds may be negligible for one predictor and consequential for another — that a package cannot default on the analyst's behalf. When you can defend limits, compute the shares outside the table with `bayestestR::rope(fit, range = c(-0.18, 0.18))` (or `bayestestR::describe_posterior()`), state the limits and their justification in the text, and keep the table's posterior summaries as they are. ## Check convergence before reading anything A frequentist table is wrong only if the model is wrong; a Bayesian table can additionally be wrong because the *sampler* has not converged — the draws then misrepresent the very posterior the columns summarize. spicy therefore checks every Bayesian fit before rendering: when any sampled parameter misses its target — R-hat at or above 1.01, effective sample size below 100 per chain (never less than 400, so fewer chains do not weaken the bar), any divergent transition, E-BFMI below 0.2 (Vehtari et al. 2021) — the table gains a `Sampler diagnostics:` footer line naming the failure and a warning fires (classed `spicy_bayes_diagnostics`, so scripts can catch or mute this guard selectively with `withCallingHandlers()`). A clean fit prints nothing: the publication table stays lean, and silence means the checks passed. The guard reads *all* sampled parameters — including the group-level ones a multilevel table summarizes — not just the displayed coefficients. For parameter-level triage the same three diagnostics are available as opt-in columns in all-Bayesian tables, together with a fourth, `mcse` — the Monte Carlo standard error of the displayed posterior median: ```{r diag-cols} table_regression(fit, show_columns = c("b", "ci", "rhat", "mcse")) ``` The `mcse` column answers a question the other diagnostics do not: *how many of the printed digits are real?* MCMC output is stochastic — rerun the sampler with another seed and the third significant digit of most estimates will move. The working rule (Gelman, Vehtari, McElreath et al. 2026, §11.6): a displayed digit is Monte-Carlo stable when **twice the MCSE stays below one unit of that digit**. Here twice the MCSE of the coefficients is on the order of their second decimal: that digit sits near the stability edge — honest to report, and a reminder that a two-decimal table is close to the resolution limit of this deliberately short chain. Most of the time two *significant* digits are all a posterior summary can support, halving the MCSE costs four times the iterations, and tail quantities (the CrI bounds) are noisier than the median. When the displayed digits matter — a boundary against a clinical threshold, say — check `mcse` before trusting the last one. Under the hood these are three numbers per parameter (Vehtari et al. 2021): R-hat, which compares chains (at convergence, 1.00; worry above roughly 1.01), and two effective sample sizes that discount autocorrelation — bulk ESS, which governs point estimates such as the posterior median, and tail ESS, which governs the credible-interval bounds: ```{r convergence} posterior::summarise_draws( posterior::subset_draws(posterior::as_draws_array(fit), variable = names(fit$coefficients)) )[, c("variable", "rhat", "ess_bulk", "ess_tail")] ``` Read the output rather than glossing it. R-hat sits at 1.00 for every slope but prints 1.01 for the intercept — exactly the alert level, and exactly what a deliberately short two-chain run invites; a manuscript fit would simply use the defaults. Bulk effective sizes exceed the 1,000 retained draws — legitimate, since Hamiltonian Monte Carlo can produce anticorrelated draws that beat independence — while the tail effective sizes, near 700–800, are what guard the 2.5% and 97.5% bounds the CI column reports: comfortable for two-decimal reporting. When these diagnostics genuinely fail, the remedy is more iterations or a reparameterized model — never a more optimistic reading of the table. ## Priors are part of the model — disclose them The posterior blends data with priors, so a reported Bayesian model is incompletely specified until the priors are stated. `rstanarm`'s defaults are weakly informative (centred at zero, scaled to the data), which is why the table above sits so close to its frequentist twin — but that is a property to *verify and report*, not assume: ```{r priors} prior_summary(fit) ``` The *Adjusted prior* line is the disclosure that matters: `rstanarm` rescales the default `normal(0, 2.5)` by each predictor's spread (scale 2.5/sd(x)), so the priors actually applied here are normal(0, 5.00) on the sex contrast and normal(0, 0.17) on age — equivalent to placing normal(0, 2.5) on the coefficients of *standardized* predictors. For a manuscript, one sentence suffices ("rstanarm defaults: independent normal(0, 2.5) priors, autoscaled by each predictor's SD — the *Adjusted prior* scales above"), and a sensitivity re-fit with wider or substantive priors belongs in the appendix when the priors could plausibly matter. ## Odds ratios: `exponentiate = TRUE` The ratio scale works exactly as for `glm`, with one pleasant difference: every displayed quantity is computed from the exponentiated draws themselves, so all of them are exact posterior statements about the odds ratio — the credible bounds are the exponentiated quantiles, and the SE is the posterior MAD SD of `exp(draws)` rather than the delta-method approximation (SE_OR = OR × SE_log-odds), which understates the ratio-scale spread severalfold when the posterior is wide. Ratio-scale posteriors are right-tail-heavy, so the interval remains the primary uncertainty summary — never reconstruct one as B ± 2 SE: ```{r or} table_regression(fit, exponentiate = TRUE) ``` Reading the `Tertiary` row with the factor-change template: holding sex and age constant, the odds of smoking for a tertiary-educated respondent are 0.40 times the odds for one with lower secondary education — and given model, priors, and data, that odds ratio lies between 0.27 and 0.57 with 95% probability. ## Marginal effects, from the draws Odds ratios answer on the odds scale; the **average marginal effect** answers on the probability scale — here, percentage points of P(smoking). For a posterior, the AME is computed *per draw* (`marginaleffects::avg_slopes()`), and the table summarizes those draws under the same conventions as every other column: posterior median, MAD SD, and a credible interval that follows `ci_method` (equal-tailed by default, HDI on request): ```{r ame} table_regression(fit, show_columns = c("b", "ame", "ame_ci")) ``` Reading the `Tertiary` row: averaged over the sample, tertiary education shifts the probability of smoking by the displayed amount — given model, priors, and data, that shift lies in its credible interval with 95% probability, a statement no delta-method AME can make exactly. There is no `ame_p`, for the same reason there is no p column: the preset `"all_ame"` expands without it, the atomic token is refused, and in a mixed frequentist–Bayesian table the shared `ame_p` column dashes the Bayesian rows. ## The frequentist twin, side by side Because the layout machinery is shared, the comparison table is one `list()` away — and with weak priors and n = 1,175 it mostly documents agreement, which is itself worth showing to a skeptical reader: ```{r side-by-side} gf <- glm(smoking ~ sex + age + education, data = sh, family = binomial) table_regression(list("ML (glm)" = gf, "Bayes (stan_glm)" = fit), show_columns = c("b", "ci")) ``` The intervals differ in *reading*, not much in width: the frequentist column's CI is a coverage statement about the procedure, the Bayesian column's CrI a probability statement about the parameter. The shared header stays `95% CI` here — the honest common label — and the footer discloses per model that the Bayesian column is an equal-tailed posterior credible interval (the all-Bayesian table above relabels the header itself, `95% CrI`). In the fit-statistics block, `n` and `R² (Bayes)` fill per model, but the likelihood-based rows — the two pseudo-R² and AIC — fill only the frequentist one, because a posterior has none of them. Which brings us to what the table refuses. ## Fit statistics for a posterior The default fit block reports `R² (Bayes)` — the posterior median of the draws-based R² (Gelman et al. 2019), the one R² a posterior genuinely has. Predictive accuracy is opt-in, because it costs a few seconds of PSIS-LOO per model: ```{r bayes-fit-stats} table_regression(fit, show_fit_stats = c("nobs", "r2_bayes", "elpd_loo")) ``` The `ELPD (LOO)` row is the expected log predictive density (Vehtari, Gelman & Gabry 2017), with its standard error disclosed in the footer — compare models with `loo::loo_compare()` outside the table, and read differences against the SE of the *difference*. Prefer the log scale: `"looic"` (its −2 twin) persists for historical reasons only (Gelman, Vehtari, McElreath et al. 2026, App. B), and `"waic"` exists with the same disclosures though PSIS-LOO is generally preferred. ## Standardized coefficients from the draws The algebraic flavors run natively on the draws — an exact affine rescale, so the median, MAD SD and credible bounds all move by the same SD ratio and nothing is re-summarized: ```{r bayes-beta} table_regression(fit, standardized = "posthoc", show_columns = c("b", "beta")) ``` On this logistic model the betas are x-standardized on the link scale (the glm convention): `age`, the only continuous predictor, moves from a per-year log-odds of 0.01 to 0.09 per SD of age — the dummies are left unscaled under `"posthoc"`, so their rows match `B`. The same call works on standard-formula `brm()` fits (the design matrix is recovered through `insight`), and the scale factors are engine-invariant across rstanarm and brms. The refusals — multilevel and special-term formulas, `"refit"`, `"pseudo"` — are itemized in the next section. ## What is refused, and why Several familiar requests are meaningless for a posterior, and spicy refuses them with an explanatory error rather than rendering something silently wrong: * `p_adjust` — there are no p-values to adjust; * likelihood-based fit statistics (`"aic"`, `"bic"`, pseudo-R²) — a posterior has no likelihood-based information criteria; the Bayesian block of the *Fit statistics* section above replaces them. The predictive estimates come with their own reliability diagnostics, and spicy never mutes them: when PSIS-LOO flags observations with Pareto k above the sample-size-specific threshold — min(1 − 1/log10(S), 0.7), the same bound `loo::loo()` itself prints (Vehtari et al. 2024) — or WAIC flags p_waic above 0.4, the footer says the estimate is unreliable and a `spicy_bayes_diagnostics` warning fires. The remediation ladder (Gelman, Vehtari, McElreath et al. 2026, ch. 24): `loo::loo_moment_match()` first, then refitting the flagged folds (`k_threshold = 0.7`, brms `reloo = TRUE`), then K-fold cross-validation. Model *comparison* belongs to `loo::loo_compare()` outside the table: the comparison uncertainty is the standard error of the elpd *difference* (not the per-model SEs shown here), and differences smaller than about 4 carry no practical meaning (Gelman, Vehtari, McElreath et al. 2026, ch. 9). The canonical fit *check* is the graphical posterior predictive check (`pp_check(fit)`; Gelman et al. 2013, ch. 6) — Bayes factors stay out by design (Gelman et al. 2013, §7.4); * robust / cluster-robust `vcov` — a sandwich estimator corrects a misspecified likelihood's standard errors; nothing standard plays that role for a posterior (misspecification-robust Bayesian procedures exist but remain research-grade), so model the clustering instead (next section); * `ci_method = "profile"` / `"boot_percentile"` — profile likelihood and bootstrap replicates are frequentist constructions with no posterior analogue (the credible interval is computed from the draws; `"hdi"` is the one alternative flavor); * `standardized = "refit"` and `"pseudo"` — refit would re-run the MCMC sampler on z-scored data inside a table call (minutes per model), and pseudo's latent variance varies per draw (its draws-native version is planned). The algebraic flavors are demonstrated in *Standardized coefficients from the draws* above; their conventions in brief: `"posthoc"` / `"basic"` match `effectsize::standardize_parameters()` to numerical precision on `stanreg` fits (`"basic"` on `brmsfit` too), and `"smart"` follows Gelman's (2008) convention of leaving binary inputs and factor dummies raw (`effectsize`'s `two_sd = TRUE` variant 2-SD-scales binary numerics). Multilevel fits (`stan_glmer`, `brm()` with group-level terms) are refused — a standardized β would need an explicit sd(Y) decomposition, the choice the frequentist mixed engines make by refitting on z-scored data — and so are `stan_polr` / `stan_betareg` (no validated convention) and brms formulas with distributional, multivariate or special terms (`mo()`, `s()`, ...). Pre-standardizing predictors before fitting (Gelman, Hill & Vehtari 2020, ch. 12) remains the cleanest route when the whole analysis should live on the standardized scale (and the route for every refused case); * variational and optimizing fits (`stan_glm(..., algorithm = "meanfield")`, `"optimizing"`) — they carry an *approximate* posterior on which the MCMC diagnostics above would be vacuous; refit with `algorithm = "sampling"`; * a ROPE column — deliberately absent, as discussed in the `pd` section: equivalence limits are substantive claims the analyst must justify per coefficient (`bayestestR::rope()` renders them once justified). ```{r refused, error = TRUE} table_regression(fit, p_adjust = "holm") ``` ## Multilevel models: random effects from the draws Clustered data is where the Bayesian machinery earns its keep — the priors regularize the group-level variance exactly where a likelihood-based fit struggles (few groups, boundary estimates). The random effects render as the same block as for `lmer`, with every number taken directly from the posterior draws: the estimate is the posterior median of each SD, the interval its posterior quantiles — and there is no boundary-corrected test because nothing is being tested — the chi-bar-squared machinery exists only because the frequentist null σ = 0 sits on the edge of the parameter space, a situation a posterior quantile never encounters: ```{r multilevel} mfit <- stan_lmer(Reaction ~ Days + (Days | Subject), data = lme4::sleepstudy, iter = 1000, chains = 2, seed = 7, refresh = 0) table_regression(mfit) ``` Note the guard from the convergence section doing its job: this deliberately trimmed run (two short chains, a correlated random slope) misses the R-hat and ESS targets, so the table says so in its footer and a `spicy_bayes_diagnostics` warning fires — a manuscript fit with rstanarm's default sampler settings renders this same table clean. Compare this block with the `lmer` table of the [mixed-effects vignette](table-regression-mixed.html): same layout, but the credible interval on the `Days` slope SD — the `σ Subject Days` row — is a genuine posterior statement — no chi-bar-squared correction, no profile refit — and the footer reads `Random effects (MCMC)`. The per-group deviations themselves (one per subject) are deliberately *not* listed as coefficient rows; they are shrinkage estimates, and their summary is the block's SD. ## `brms` `brms` reaches the same table through the same code path; the only practical difference is that each model compiles before sampling (minutes, and a C++ toolchain), which is why this chunk is shown without being evaluated. Run it locally to see the table — it matches the `rstanarm` one above row for row, as the package's shared frame tests verify: ```{r brms, eval = run_brms} bfit <- brms::brm(Reaction ~ Days + (Days | Subject), data = lme4::sleepstudy, chains = 2, iter = 1000, seed = 7, refresh = 0, silent = 2) table_regression(bfit) ``` Coefficients are posterior medians with equal-tailed credible intervals, the random-effect block is built from the `sd_*` / `cor_*` draws, and the refusals above apply identically. What `brms` buys is model space: distributional regression, non-linear formulas, and families neither `glm` nor `rstanarm` covers — all rendered through the same table. ## Output formats The default console table is one of several targets. The `output` argument also produces a raw data frame, a long broom-style tibble, and — with the corresponding Suggests package — rich `gt`, `flextable`, `tinytable`, Excel, or Word tables. The Bayesian structure carries through: an all-Bayesian table carries no p column at all (in a mixed frequentist–Bayesian table the shared p column dashes on the Bayesian side), and the posterior summaries keep full precision in the tidy frame. In [`broom::tidy()`](https://broom.tidymodels.org) the `estimate` column is the posterior median, `std.error` the posterior MAD SD, and `conf.low` / `conf.high` the credible bounds; `p.value`, `statistic`, and `df` are `NA` — pipelines that filter on `p.value` must switch to the interval (or `pd`): ```{r broom} td <- broom::tidy(table_regression(fit)) td[, c("term", "estimate", "std.error", "conf.low", "conf.high", "p.value")] ``` ```{r outputs-gt, eval = build_rich_tables} table_regression(fit, exponentiate = TRUE, output = "gt") ``` ## See also - `vignette("table-regression", package = "spicy")` for the shared mechanics (columns, layouts, filtering). - `vignette("table-regression-mixed", package = "spicy")` for the frequentist multilevel counterpart of the random-effects block. ## References * Gelman, A. (2008). Scaling regression inputs by dividing by two standard deviations. *Statistics in Medicine*, 27(15), 2865–2873. * Gelman, A., Hill, J., & Vehtari, A. (2020). *Regression and Other Stories*. Cambridge University Press. * Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). *Bayesian Data Analysis* (3rd ed.). CRC Press. * Gelman, A., Goodrich, B., Gabry, J., & Vehtari, A. (2019). R-squared for Bayesian regression models. *The American Statistician*, 73(3), 307–309. * Gelman, A., Vehtari, A., McElreath, R., et al. (2026). *Bayesian Workflow*. Chapman & Hall / CRC Press. * Kruschke, J. K. (2015). *Doing Bayesian Data Analysis* (2nd ed.). Academic Press. * Kruschke, J. K. (2021). Bayesian analysis reporting guidelines. *Nature Human Behaviour*, 5(10), 1282–1291. * Makowski, D., Ben-Shachar, M. S., Chen, S. H. A., & Lüdecke, D. (2019). Indices of effect existence and significance in the Bayesian framework. *Frontiers in Psychology*, 10, 2767. * Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. *Statistics and Computing*, 27(5), 1413–1432. * Vehtari, A., Simpson, D., Gelman, A., Yao, Y., & Gabry, J. (2024). Pareto smoothed importance sampling. *Journal of Machine Learning Research*, 25(72), 1–58. * Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC. *Bayesian Analysis*, 16(2), 667–718.