Bayesian regression tables
Source:vignettes/articles/table-regression-bayesian.Rmd
table-regression-bayesian.RmdThis article covers Bayesian regression tables — the
same table_regression() call on a posterior instead of a
maximum-likelihood fit. The companion article Publication-ready regression
tables covers the shared mechanics (the class-by-class map is
Supported
models); 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 everyrstanarmchunk in this article 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 article’s glm section — same data
(sochealth), same predictors — so the Bayesian table can be
read against its frequentist twin.
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 article 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
table_regression(fit)
#> Bayesian logistic regression (stanreg): smoking
#>
#> Variable │ B SE 95% CrI
#> ──────────────────────────┼───────────────────────────────
#> (Intercept) │ -1.11 0.29 [-1.66, -0.55]
#> sex: │
#> Female (ref.) │ – – –
#> Male │ -0.04 0.15 [-0.34, 0.25]
#> age │ 0.01 0.01 [-0.00, 0.02]
#> education: │
#> Lower secondary (ref.) │ – – –
#> Upper secondary │ -0.49 0.16 [-0.81, -0.16]
#> Tertiary │ -0.91 0.20 [-1.31, -0.56]
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Bayes) │ 0.02
#>
#> Note. Bayesian logistic regression (stanreg).
#> Std. errors: posterior MAD SD (scaled median absolute deviation).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
rstanarmconvention). 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 ofrstanarm’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 (theTertiaryrow). 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 matchesrstanarm; 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 header95% 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:
table_regression(fit, show_columns = c("b", "ci", "pd"))
#> Bayesian logistic regression (stanreg): smoking
#>
#> Variable │ B 95% CrI pd
#> ──────────────────────────┼────────────────────────────────
#> (Intercept) │ -1.11 [-1.66, -0.55] 1.000
#> sex: │
#> Female (ref.) │ – – –
#> Male │ -0.04 [-0.34, 0.25] .617
#> age │ 0.01 [-0.00, 0.02] .886
#> education: │
#> Lower secondary (ref.) │ – – –
#> Upper secondary │ -0.49 [-0.81, -0.16] .999
#> Tertiary │ -0.91 [-1.31, -0.56] 1.000
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Bayes) │ 0.02
#>
#> Note. Bayesian logistic regression (stanreg).
#> Std. errors: posterior MAD SD (scaled median absolute deviation).
#> pd = probability of direction (share of the posterior on the dominant side of zero; Makowski et al. 2019).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 or
effective sample size below 100 per chain, never less than 400, so fewer
chains do not weaken the bar (Vehtari et al. 2021); any divergent
transition (Stan guidance); E-BFMI below 0.2 (Betancourt 2017) — 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:
table_regression(fit, show_columns = c("b", "ci", "rhat", "mcse"))
#> Bayesian logistic regression (stanreg): smoking
#>
#> Variable │ B 95% CrI R-hat MCSE
#> ──────────────────────────┼─────────────────────────────────────────
#> (Intercept) │ -1.11 [-1.66, -0.55] 1.007 0.014
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ -0.04 [-0.34, 0.25] 1.000 0.0050
#> age │ 0.01 [-0.00, 0.02] 1.002 0.00015
#> education: │
#> Lower secondary (ref.) │ – – – –
#> Upper secondary │ -0.49 [-0.81, -0.16] 1.002 0.0071
#> Tertiary │ -0.91 [-1.31, -0.56] 1.002 0.0097
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Bayes) │ 0.02
#>
#> Note. Bayesian logistic regression (stanreg).
#> Std. errors: posterior MAD SD (scaled median absolute deviation).
#> MCSE = Monte Carlo standard error of the posterior median (Vehtari et al. 2021).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:
posterior::summarise_draws(
posterior::subset_draws(posterior::as_draws_array(fit),
variable = names(fit$coefficients))
)[, c("variable", "rhat", "ess_bulk", "ess_tail")]
#> # A tibble: 5 × 4
#> variable rhat ess_bulk ess_tail
#> <chr> <dbl> <dbl> <dbl>
#> 1 (Intercept) 1.01 1301. 689.
#> 2 sexMale 1.00 1292. 743.
#> 3 age 1.00 1270. 788.
#> 4 educationUpper secondary 1.00 1064. 816.
#> 5 educationTertiary 1.00 1133. 797.Read the output rather than glossing it. R-hat sits at 1.00 for every
slope but prints as 1.01 for the intercept in this
three-significant-digit display — brushing the alert level, exactly what
a deliberately short two-chain run invites. The exact value, 1.007
(visible in the R-hat column of the spicy table above),
still clears the guard’s at-or-above-1.01 bar — which is why no
Sampler diagnostics: footer appears on any table of this
fit; 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:
prior_summary(fit)
#> Priors for model 'fit'
#> ------
#> Intercept (after predictors centered)
#> ~ normal(location = 0, scale = 2.5)
#>
#> Coefficients
#> Specified prior:
#> ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
#> Adjusted prior:
#> ~ normal(location = [0,0,0,...], scale = [5.00,0.17,5.02,...])
#> ------
#> See help('prior_summary.stanreg') for more detailsThe 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) — the same robust pairing as the log-odds table,
median + MAD SD, now applied to the OR posterior — rather than the
delta-method approximation (SE_OR = OR × SE_log-odds), which is accurate
only while the posterior is narrow (on this fit the two agree to the
displayed precision). Ratio-scale posteriors are right-tail-heavy, and
no single scale number summarizes a heavy right tail well: on a wide
posterior the delta approximation understates the ratio-scale
SD severalfold, while the MAD SD deliberately discounts the
tail. So the interval remains the primary uncertainty summary — never
reconstruct one as B ± 2 SE:
table_regression(fit, exponentiate = TRUE)
#> Bayesian logistic regression (stanreg): smoking
#>
#> Variable │ OR SE 95% CrI
#> ──────────────────────────┼─────────────────────────────
#> (Intercept) │ 0.33 0.09 [0.19, 0.57]
#> sex: │
#> Female (ref.) │ – – –
#> Male │ 0.96 0.15 [0.71, 1.28]
#> age │ 1.01 0.01 [1.00, 1.02]
#> education: │
#> Lower secondary (ref.) │ – – –
#> Upper secondary │ 0.61 0.10 [0.44, 0.85]
#> Tertiary │ 0.40 0.08 [0.27, 0.57]
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Bayes) │ 0.02
#>
#> Note. Bayesian logistic regression (stanreg).
#> Std. errors: posterior MAD SD (scaled median absolute deviation).
#> OR = odds ratio.
#> Coefficients exponentiated and displayed as OR; SE on the OR scale (posterior MAD SD of the exponentiated draws); CI bounds exponentiated (asymmetric).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, changes in
P(smoking), displayed as proportions. 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):
table_regression(fit, show_columns = c("b", "ame", "ame_ci"))
#> Bayesian logistic regression (stanreg): smoking
#>
#> Variable │ B AME 95% CrI
#> ──────────────────────────┼────────────────────────────────
#> (Intercept) │ -1.11
#> sex: │
#> Female (ref.) │ – – –
#> Male │ -0.04 -0.01 [-0.06, 0.04]
#> age │ 0.01 0.00 [-0.00, 0.00]
#> education: │
#> Lower secondary (ref.) │ – – –
#> Upper secondary │ -0.49 -0.09 [-0.16, -0.03]
#> Tertiary │ -0.91 -0.15 [-0.22, -0.09]
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Bayes) │ 0.02
#>
#> Note. Bayesian logistic regression (stanreg).
#> Std. errors: posterior MAD SD (scaled median absolute deviation).
#> AME = average marginal effect.Reading the Tertiary row: averaged over the sample,
tertiary education shifts the probability of smoking by −0.15 — 15
percentage points lower — and, given model, priors, and data, that shift
lies in [−0.22, −0.09] 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:
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"))
#> Regression comparison: smoking
#>
#> ML (glm) Bayes (stan_glm)
#> ─────────────────────── ───────────────────────
#> Variable │ B 95% CI B 95% CI
#> ──────────────────────────┼──────────────────────────────────────────────────
#> (Intercept) │ -1.11 [-1.67, -0.55] -1.11 [-1.66, -0.55]
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ -0.04 [-0.32, 0.25] -0.04 [-0.34, 0.25]
#> age │ 0.01 [-0.00, 0.02] 0.01 [-0.00, 0.02]
#> education: │
#> Lower secondary (ref.) │ – – – –
#> Upper secondary │ -0.48 [-0.82, -0.14] -0.49 [-0.81, -0.16]
#> Tertiary │ -0.91 [-1.29, -0.52] -0.91 [-1.31, -0.56]
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175 1175
#> R² (McFadden) │ 0.02 –
#> R² (Nagelkerke) │ 0.03 –
#> AIC │ 1200.9 –
#> R² (Bayes) │ – 0.02
#>
#> Note. ML (glm): logistic regression; Bayes (stan_glm): Bayesian logistic regression (stanreg).
#> Std. errors:
#> ML (glm): classical (Fisher information)
#> Bayes (stan_glm): posterior MAD SD (scaled median absolute deviation)
#> Model 2: 95% CI is an equal-tailed posterior credible interval.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 fills for both models; the remaining rows split by
paradigm — R² (Bayes) fills only the Bayesian column, and
the likelihood-based rows (the two pseudo-R² and AIC) only the
frequentist one, because a posterior has none of them. Which brings us
to the fit statistics a posterior does have — and to those it lacks.
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:
table_regression(fit, show_fit_stats = c("nobs", "r2_bayes",
"elpd_loo"))
#> Bayesian logistic regression (stanreg): smoking
#>
#> Variable │ B SE 95% CrI
#> ──────────────────────────┼───────────────────────────────
#> (Intercept) │ -1.11 0.29 [-1.66, -0.55]
#> sex: │
#> Female (ref.) │ – – –
#> Male │ -0.04 0.15 [-0.34, 0.25]
#> age │ 0.01 0.01 [-0.00, 0.02]
#> education: │
#> Lower secondary (ref.) │ – – –
#> Upper secondary │ -0.49 0.16 [-0.81, -0.16]
#> Tertiary │ -0.91 0.20 [-1.31, -0.56]
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Bayes) │ 0.02
#> ELPD (LOO) │ -600.6
#>
#> Note. Bayesian logistic regression (stanreg).
#> Std. errors: posterior MAD SD (scaled median absolute deviation).
#> Predictive accuracy by PSIS-LOO; SE(ELPD) = 18.7.The ELPD (LOO) row is the expected log predictive
density (Vehtari, Gelman & Gabry 2017), with its standard error
disclosed in the footer. 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.
These 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).
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:
table_regression(fit, standardized = "posthoc",
show_columns = c("b", "beta"))
#> Bayesian logistic regression (stanreg): smoking
#>
#> Variable │ B β
#> ──────────────────────────┼────────────────
#> (Intercept) │ -1.11 –
#> sex: │
#> Female (ref.) │ – –
#> Male │ -0.04 -0.04
#> age │ 0.01 0.09
#> education: │
#> Lower secondary (ref.) │ – –
#> Upper secondary │ -0.49 -0.49
#> Tertiary │ -0.91 -0.91
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Bayes) │ 0.02
#>
#> Note. Bayesian logistic regression (stanreg).
#> Std. errors: posterior MAD SD (scaled median absolute deviation).
#> β = standardised coefficient ("posthoc": B × SD(X)/SD(Y) for numeric predictors, B/SD(Y) for factor dummies).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 Fit statistics section above shows what replaces them (and the reliability guards that come with the replacements); - 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"matcheffectsize::standardize_parameters()to numerical precision onstanregfits ("basic"onbrmsfittoo), and"smart"follows Gelman’s (2008) convention of leaving binary inputs and factor dummies raw (effectsize’stwo_sd = TRUEvariant 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 arestan_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 withalgorithm = "sampling"; - a ROPE column — deliberately absent, as discussed in the
pdsection: equivalence limits are substantive claims the analyst must justify per coefficient (bayestestR::rope()renders them once justified).
table_regression(fit, p_adjust = "holm")
#> Error in `table_regression()`:
#> ! `p_adjust` is not available for Bayesian fits: there are no p-values to adjust.
#> ℹ Bayesian tables report posterior medians and credible intervals; the probability-of-direction column (`show_columns = "pd"`) is the closest posterior summary.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:
mfit <- stan_lmer(Reaction ~ Days + (Days | Subject),
data = lme4::sleepstudy,
iter = 1000, chains = 2, seed = 7, refresh = 0)
table_regression(mfit)
#> Warning: Bayesian fit (outcome: Reaction) shows sampler problems -- Sampler
#> diagnostics: max R-hat = 1.010 (target < 1.01); min ESS = 253 (target > 400).
#> Do not report as-is; run longer or reparameterize (Vehtari et al. 2021).
#> Bayesian linear regression (stanreg): Reaction
#>
#> Variable │ B SE 95% CrI
#> ──────────────────────────────────┼────────────────────────────────
#> (Intercept) │ 251.23 6.53 [237.54, 264.18]
#> Days │ 10.38 1.82 [ 6.58, 13.95]
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> Random effects: │
#> σ Subject (Intercept) │ 21.99 5.52 [ 12.68, 35.63]
#> σ Subject Days │ 6.53 1.33 [ 4.37, 10.18]
#> ρ Subject (Days × (Intercept)) │ 0.11 0.30 [ -0.42, 0.70]
#> σ (Residual) │ 25.99 1.63 [ 23.14, 29.28]
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 180
#> R² (Bayes) │ 0.79
#>
#> Note. Bayesian linear regression (stanreg).
#> Std. errors: posterior MAD SD (scaled median absolute deviation).
#> Random effects (MCMC).
#> Sampler diagnostics: max R-hat = 1.010 (target < 1.01); min ESS = 253 (target > 400). Do not report as-is; run longer or reparameterize (Vehtari et al. 2021).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 article: 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:
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() 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):
td <- broom::tidy(table_regression(fit))
td[, c("term", "estimate", "std.error", "conf.low", "conf.high",
"p.value")]
#> # A tibble: 5 × 6
#> term estimate std.error conf.low conf.high p.value
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) -1.11 0.287 -1.66 -0.554 NA
#> 2 sexMale -0.0380 0.153 -0.337 0.249 NA
#> 3 age 0.00614 0.00504 -0.00326 0.0152 NA
#> 4 educationUpper secondary -0.487 0.165 -0.811 -0.164 NA
#> 5 educationTertiary -0.909 0.197 -1.31 -0.559 NA
pkgdown_dark_gt(table_regression(fit, exponentiate = TRUE, output = "gt"))| Bayesian logistic regression (stanreg): smoking | ||||
|
Variable
|
OR
|
SE
|
95% CrI
|
|
|---|---|---|---|---|
| LL | UL | |||
| (Intercept) | 0.33 | 0.09 | 0.19 | 0.57 |
| sex: | ||||
| Female (ref.) | – | – | – | – |
| Male | 0.96 | 0.15 | 0.71 | 1.28 |
| age | 1.01 | 0.01 | 1.00 | 1.02 |
| education: | ||||
| Lower secondary (ref.) | – | – | – | – |
| Upper secondary | 0.61 | 0.10 | 0.44 | 0.85 |
| Tertiary | 0.40 | 0.08 | 0.27 | 0.57 |
| n | 1175 | |||
| R² (Bayes) | 0.02 | |||
See also
- Publication-ready regression tables for the shared mechanics (columns, layouts, filtering).
- Mixed-effects regression tables for the frequentist multilevel counterpart of the random-effects block.
References
- Betancourt, M. (2017). A conceptual introduction to Hamiltonian Monte Carlo. arXiv preprint, arXiv:1701.02434.
- 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.