Publication-ready regression tables
Source:vignettes/articles/table-regression.Rmd
table-regression.Rmdtable_regression() produces a coefficient summary table
from one or several fitted regression models. The output is
publication-ready by default and follows APA Manual 7 (American
Psychological Association 2020, Tables 7.13–7.15) formatting
conventions: paired estimate-and-CI columns, APA p-values without
leading zero, factor levels grouped under their parent variable, fit
statistics at the foot of the table, and a self-documenting note line
that names the variance estimator and any methodological choice that
affected the rendered values.
APA is the default, not a limit: style = "nejm" re-types
every p-value to the New England Journal of Medicine’s statistical
reporting guideline, style = "lancet" writes The Lancet’s
midline decimal point and en-dash intervals, and "jama",
"annals", and "aer" apply their journals’
rules the same way. Every rule a style encodes is sourced:
?spicy_style quotes the sentence of the author guideline it
comes from.
The function is fit-first: you pass already-fitted
models, not raw data and a formula. The same object exports cleanly to a
long broom-canonical frame (broom::tidy()) and to a
one-row-per-model glance summary (broom::glance()).
This article teaches the shared mechanics on lm() and
glm() fits. The Generalized linear models section
covers the glm-specific argument semantics; the Mixed-effects
models section introduces the Random effects rows that
mixed-effects fits add below the fixed effects (lmer,
glmer, glmmTMB, lme, and the
ordinal clmm). Each further model family has a dedicated
article that applies these mechanics to its own estimands: Mixed-effects regression
tables for mixed-effects models, GEE regression tables for
population-averaged (GEE) models, Count and two-part regression
tables for Poisson, negative-binomial, and two-part models, Ordinal regression
tables for ordinal models (polr, clm,
and the mixed clmm), Multinomial regression
tables for multinomial models (multinom,
mlogit), Survival regression
tables for survival models (coxph,
survreg), and Bayesian regression
tables for Bayesian fits (rstanarm,
brms). The mechanics covered here — confidence levels,
multi-model and nested layouts, p-value adjustment, coefficient
filtering, output formats, and the broom methods — apply across every
supported class; variance-estimator availability varies by family, and
the complete map — every supported class with its estimands, variance
estimators and fit statistics, including the families without a
dedicated article (quantile regression, fixest,
betareg, …) — is Supported
models. One decision cuts across all of them — how to code,
test, and present categorical predictors — and has its own
article: Categorical
predictors.
Basic usage
Pass a fitted lm() object. The default rendering returns
a single-model table with B, SE,
95% CI, and p columns and a fit-statistics
footer (n, R², Adj. R²):
fit <- lm(wellbeing_score ~ age + sex + smoking, data = sochealth)
table_regression(fit)
#> Linear regression: wellbeing_score
#>
#> Variable │ B SE 95% CI p
#> ─────────────────┼──────────────────────────────────────
#> (Intercept) │ 65.20 1.66 [61.95, 68.45] <.001
#> age │ 0.05 0.03 [-0.01, 0.11] .130
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 3.86 0.91 [ 2.08, 5.63] <.001
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ -1.72 1.11 [-3.89, 0.45] .121
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² │ 0.02
#> Adj. R² │ 0.02
#>
#> Note. Linear regression. Std. errors: classical (OLS). CI = confidence interval.Reading the table:
- Each predictor occupies one row; factor predictors are grouped under
their parent variable name with one row per level. The reference level
carries the
(ref.)annotation and an en dash across the statistic columns, so the substantive comparison is visible in a single glance (NEJM / BMJ clinical convention). - The
95% CIcolumn reports a symmetric Wald interval at the level set byci_level(default 0.95). The CI header tracks the chosen level — switching toci_level = 0.99re-labels the column accordingly. - The footer line names the variance estimator in plain English so the
reader can find the inferential regime without leaving the table. The
estimator switches with
vcov(covered in Robust variance below). - Read substantively, one sentence per coefficient, from the exact
numbers shown:
Male— B = 3.86, 95% CI [2.08, 5.63], p < .001 — holding age and smoking constant, men score on average 3.9 points higher than women on the 0–100 WHO-5 wellbeing index.age— B = 0.05, 95% CI [-0.01, 0.11], p = .130 — the interval spans effects from slightly negative to about a tenth of a point per year, so the data do not establish an adjusted age gradient.
One argument re-typesets the same table for another journal — here
NEJM’s statistical-reporting rules: leading zero kept, two decimals
above 0.01, floor <0.001:
table_regression(fit, style = "nejm")
#> Linear regression: wellbeing_score
#>
#> Variable │ B SE 95% CI p
#> ─────────────────┼───────────────────────────────────────
#> (Intercept) │ 65.20 1.66 [61.95, 68.45] <0.001
#> age │ 0.05 0.03 [-0.01, 0.11] 0.13
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 3.86 0.91 [ 2.08, 5.63] <0.001
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ -1.72 1.11 [-3.89, 0.45] 0.12
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² │ 0.02
#> Adj. R² │ 0.02
#>
#> Note. Linear regression. Std. errors: classical (OLS). CI = confidence interval.The complete tour of the named styles — what each journal specifies, with the guideline sentence behind every rule — is in the House styles section of Summary tables for reporting.
Standardized coefficients
Standardized coefficients (β) make predictors with
different natural scales comparable: a one-standard-deviation increase
in X predicts a β-standard-deviation change in
Y. APA Manual 7 §7.13 recommends reporting both
B and β so the unstandardized effect (natural
units, interpretable) stays alongside the standardized effect
(comparable across predictors).
standardized selects the method. Four apply to linear
models — a fifth, glm-specific "pseudo", is covered in the
Generalized linear models section — and the choice is
consequential and well-documented (Cohen, Cohen, West, and Aiken 2003
§3.4; Gelman 2008). The methods differ most visibly in how they treat
factor dummies:
-
"refit"— refit the model on the z-scored outcome and z-scored numeric predictors; factors stay factors, so their dummies keep the 0/1 scale (the Cohen et al. 2003 convention, matching theeffectsizepackage). A dummy’sβis its group contrast inSD(Y)units. -
"posthoc"— the same estimand as"refit", obtained by algebraic rescaling of the original fit:β = B × SD(X) / SD(Y)for numeric predictors,β = B / SD(Y)for factor dummies (which stay on 0/1). Numerically identical to"refit"for purely linear-additive Gaussian models; preferred when refitting is expensive or whenlm()was wrapped in a pipeline that resists re-execution. -
"basic"— algebraic, but every design column is scaled by its column SD, factor dummies included:β = B × SD(column) / SD(Y). This is the Beta definition of SPSSREGRESSIONand Stataregress, beta(spicy’s oracle tests pin it against PSPP): reach for it when you need numbers that reproduce those packages, at the price of a dummyβthat depends on the level’s base rate throughSD(dummy). -
"smart"— Gelman’s (2008) input scaling: numeric predictors divided by2 × SD(X); binary predictors and factor dummies kept on their 0/1 scale. The resultingβis the change inY, inSD(Y)units, for a two-standard-deviation change inX(from mean − SD to mean + SD), which puts a continuous predictor’s ±1 SD swing on the same footing as a binary’s 0 → 1 step. Where the paper leavesYraw, spicy divides bySD(Y)soβstays comparable across the four methods.
When standardized != "none", the "beta"
token is auto-injected into show_columns immediately after
"b":
fit <- lm(wellbeing_score ~ age + sex + smoking, data = sochealth)
table_regression(fit, standardized = "refit")
#> Linear regression: wellbeing_score
#>
#> Variable │ B β SE 95% CI p
#> ─────────────────┼─────────────────────────────────────────────
#> (Intercept) │ 65.20 -0.10 1.66 [61.95, 68.45] <.001
#> age │ 0.05 0.04 0.03 [-0.01, 0.11] .130
#> sex: │
#> Female (ref.) │ – – – – –
#> Male │ 3.86 0.25 0.91 [ 2.08, 5.63] <.001
#> smoking: │
#> No (ref.) │ – – – – –
#> Yes │ -1.72 -0.11 1.11 [-3.89, 0.45] .121
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² │ 0.02
#> Adj. R² │ 0.02
#>
#> Note. Linear regression. Std. errors: classical (OLS). β = standardized coefficient ("refit": outcome and numeric predictors z-scored, factor dummies on 0/1); CI = confidence interval.On the standardized scale the predictors become directly comparable:
being male predicts a 0.25-SD higher wellbeing score, the largest
standardized effect in the model, while a one-SD increase in age
predicts a change of only 0.04 SD. The raw B column could
not support that ranking — 3.86 versus 0.05 compares a factor contrast
with a per-year slope, two different units.
Caveat: standardized coefficients with interactions or
transformed terms. When the model contains a product term, an
I(), poly(), log(), or
splines::ns() wrapper, the standardized coefficient of the
non-additive term has no closed-form “one-SD change in X” reading (Aiken
and West 1991; Cohen et al. 2003 §7.7). The function emits a classed
spicy_caveat warning at runtime and prints a
method-specific caveat line in the table footer, so the limitation is
exposed at the point of use without blocking the table.
Per-coefficient effect sizes
For each predictor, you can request a partial effect-size column.
Three measures are available — Cohen’s f²
(partial_f2), Pearson’s partial η²
(partial_eta2), and the Olejnik–Algina bias-corrected
partial ω² (partial_omega2). Each carries a
confidence interval derived from noncentral-F inversion
(Smithson 2003; Steiger 2004), exposed as a separate
<token>_ci column. The shortcuts
"all_f2", "all_eta2", and
"all_omega2" each expand to the point estimate together
with its CI column. The partial F the three are built on is a column of
its own, partial_f, printed as value (q, df):
the term’s degrees of freedom, then the residual degrees of freedom of
the model.
sochealth$education is stored as an ordered factor,
which R codes with polynomial contrasts (.L = linear,
.Q = quadratic) — trend components across the ordered
levels, not per-level effects against a reference. For a coefficient
table read level-against-reference, convert it to a plain factor; we
reuse this copy throughout the article:
sh <- sochealth |>
dplyr::mutate(education = factor(education, ordered = FALSE))
fit <- lm(wellbeing_score ~ age + sex + smoking + education,
data = sh)
table_regression(
fit,
show_columns = c("b", "p", "all_eta2", "all_omega2")
)
#> Linear regression: wellbeing_score
#>
#> Variable │ B p η² η² 95% CI ω²
#> ──────────────────────────┼──────────────────────────────────────────
#> (Intercept) │ 53.96 <.001
#> age │ 0.03 .344 0.00 [0.00, 0.00] 0.00
#> sex: │
#> Female (ref.) │ – –
#> Male │ 3.57 <.001 0.02 [0.00, 0.03] 0.02
#> smoking: │
#> No (ref.) │ – –
#> Yes │ 0.68 .496 0.00 [0.00, 0.00] 0.00
#> education: │
#> Lower secondary (ref.) │ – –
#> Upper secondary │ 11.89 <.001 0.21 [0.17, 0.25] 0.21
#> Tertiary │ 19.71 <.001 0.21 [0.17, 0.25] 0.21
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² │ 0.22
#> Adj. R² │ 0.22
#>
#> Variable │ ω² 95% CI
#> ──────────────────────────┼──────────────
#> (Intercept) │
#> age │ [0.00, 0.00]
#> sex: │
#> Female (ref.) │
#> Male │ [0.00, 0.03]
#> smoking: │
#> No (ref.) │
#> Yes │ [0.00, 0.00]
#> education: │
#> Lower secondary (ref.) │
#> Upper secondary │ [0.17, 0.25]
#> Tertiary │ [0.17, 0.25]
#>
#> Note. Linear regression. Std. errors: classical (OLS). η² = partial eta squared; CI = confidence interval; ω² = partial omega squared.Methodology notes:
- The partial F-test is a Type-II test (marginality-respecting,
contrast-invariant; Fox 2016; Fox and Weisberg 2019): the focal term is
tested by comparing the two nested models that both exclude every
higher-order term containing it —
Ainy ~ A * BisSS(A | B), withA:Bleft out of both sides — over the full model’s mean squared error. This matchescar::Anova(type = 2). For additive models it reduces to the classic “each term after all the others” partial F, and the result does not depend on the factor coding (treatment, sum, Helmert, polynomial). Readers coming from SPSSUNIANOVAor SASPROC GLM, whose default is Type III, should note the difference is a deliberate convention choice, not an oversight; the two conventions coincide for models without interactions. - The
ω²point estimate is bias-corrected via the Olejnik and Algina (2003) formula((F − 1) × df1) / (F × df1 + N − df1), yielding a less-biased small-sample estimator than partialη². - The CI bounds come from inverting the noncentrality parameter of the
F-distribution at the lower and upper confidence levels (Steiger 2004
§4), anchored at the bias-corrected
Fimplied byω²rather than at the observedF. Population partialη²and partialω²are the same quantity —ω²is simply the less-biased small-sample estimator — so both columns share that single interval, and it always brackets the bias-corrected point estimate, even when the lower bound clips at zero (common for near-null terms). Expect small differences from theeffectsizepackage on both columns: it inverts the observedFfor itsη²bounds and its own impliedFforω². - For factor predictors with
klevels, the partial F-test is the joint(k − 1)df Wald test, so the same effect-size value is broadcast across all non-reference dummy rows; the reference row leaves the effect-size cells blank.
Average marginal effects (AME)
The B coefficient is the model’s structural
effect — the change in the linear predictor Xβ for a
one-unit change in x, holding the other predictors
constant. The average marginal effect (AME) is the
observable effect on the response. For a numeric predictor it
is, by default, the average partial derivative dE[Y|X] / dx
over the sample; for a factor level it is the average discrete
change — the difference in predicted response between that
level and the reference, averaged over the sample (Long and Freese 2014
§6.2 distinguish the two). For a purely linear-additive
lm() both coincide with B. For models with
interactions or transformed terms, and for any glm() with a
non-identity link, the AME is the quantity the substantive reader
actually needs. The average in the name is load-bearing: the
effect is computed for every observation and then averaged over the
sample, not evaluated once at the covariate means — the marginal effect
at the mean (MEM), a point at which no actual observation may sit. The
two generally differ in nonlinear models; Long and Freese (2014, §6.2.3)
discuss the choice in detail.
Add "ame" to show_columns for the point
estimate. Companion tokens display the SE, CI, and p-value:
"ame_se", "ame_ci", "ame_p". The
shortcut "all_ame" expands to the full set.
fit <- lm(wellbeing_score ~ age + sex + smoking, data = sochealth)
table_regression(
fit,
show_columns = c("b", "se", "p", "ame", "ame_ci", "ame_p")
)
#> Linear regression: wellbeing_score
#>
#> Variable │ B SE p AME 95% CI p
#> ─────────────────┼───────────────────────────────────────────────────
#> (Intercept) │ 65.20 1.66 <.001
#> age │ 0.05 0.03 .130 0.05 [-0.01, 0.11] .130
#> sex: │
#> Female (ref.) │ – – – – – –
#> Male │ 3.86 0.91 <.001 3.86 [ 2.08, 5.63] <.001
#> smoking: │
#> No (ref.) │ – – – – – –
#> Yes │ -1.72 1.11 .121 -1.72 [-3.89, 0.45] .121
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² │ 0.02
#> Adj. R² │ 0.02
#>
#> Note. Linear regression. Std. errors: classical (OLS). AME of a continuous predictor: average slope (dY/dX). AME = average marginal effect; CI = confidence interval.Reading conventions:
-
"p"always refers to theB(orβ) coefficient, never to the AME. Use"ame_p"for the AME-specific p-value. - Placing
"ame"after"p"keeps the attribution of each p-value unambiguous — the B-block (B / SE / p) is closed before the AME-block opens. - For non-linear models the AME is reported on the response scale, not the link scale. The Generalized linear models section below works out a logistic-regression example.
- The note names what the AME of a continuous predictor measures: here
the average slope (dY/dX), the default. In an
lm()without transformed terms it is also the change for one more unit of the predictor. Under a non-identity link the two differ, andame_typechooses between them (Slope or discrete change, below).
Inference is delegated to marginaleffects::avg_slopes(),
or to avg_comparisons() for the changes
ame_type asks for, and respects the chosen variance
estimator. Under cluster-robust variance (covered in Robust
variance), the AME inference shares the coefficient’s
t-distribution and Satterthwaite-corrected degrees of freedom via
clubSandwich::linear_contrast(), so B and AME
are reported on the same inferential footing in the same table.
marginaleffects is a Suggests dependency with two distinct
failure modes: requesting AME columns without the package
installed is a hard, classed error (spicy_missing_pkg,
with the install hint), while a computation failure inside
marginaleffects on an installed system warns
(spicy_fallback) and dashes the AME cells, leaving the rest
of the table intact.
Multiple models side by side
Pass a list of lm() fits. The default column layout
gives each model its own group of columns under a centered
spanner label showing the model name; sub-columns
(B / SE / p) are repeated under the spanner. When dependent
variables differ across models and the user did not supply labels
(model_labels = or named list), each outcome is lifted into
the spanner automatically: its variable name by default, its label
attribute with labels = TRUE (see Display
options):
m_wellbeing <- lm(wellbeing_score ~ age + sex + smoking,
data = sochealth)
m_bmi <- lm(bmi ~ age + sex + smoking,
data = sochealth)
table_regression(list(m_wellbeing, m_bmi))
#> Linear regression comparison
#>
#> wellbeing_score bmi
#> ──────────────────── ────────────────────
#> Variable │ B SE p B SE p
#> ─────────────────┼────────────────────────────────────────────
#> (Intercept) │ 65.20 1.66 <.001 23.98 0.40 <.001
#> age │ 0.05 0.03 .130 0.04 0.01 <.001
#> sex: │
#> Female (ref.) │ – – – – – –
#> Male │ 3.86 0.91 <.001 0.51 0.22 .018
#> smoking: │
#> No (ref.) │ – – – – – –
#> Yes │ -1.72 1.11 .121 -0.06 0.26 .822
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175 1163
#> R² │ 0.02 0.02
#> Adj. R² │ 0.02 0.02
#>
#> Note. Linear regression models. Std. errors: classical (OLS).When all models share the same DV, the DV appears in the title. Use a
named list — e.g. list(Crude = m1, Adjusted = m2) — to set
the spanner labels explicitly; pass model_labels = c(...)
to override the names from the list. When such labels replace the names
of different DVs, an “Outcome” row says which model has which
(outcome_labels = FALSE removes it).
Models one under the other
Side by side, models that share few predictors leave a wide and
sparse table, and the default drops the confidence intervals to save
width. layout = "stacked" sets the models one under the
other instead: each model gets a panel with the columns it would have on
its own, its interval included, opened by a row with its label. As side
by side, models with different outcomes, labeled here by their names in
the list, get an “Outcome” row under each label
(outcome_labels = FALSE removes it).
table_regression(
list(Wellbeing = m_wellbeing, BMI = m_bmi),
layout = "stacked"
)
#> Linear regression comparison
#>
#> Variable │ B SE 95% CI p
#> ─────────────────┼──────────────────────────────────────
#> Wellbeing: │
#> Outcome │ wellbeing_score
#> (Intercept) │ 65.20 1.66 [61.95, 68.45] <.001
#> age │ 0.05 0.03 [-0.01, 0.11] .130
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 3.86 0.91 [ 2.08, 5.63] <.001
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ -1.72 1.11 [-3.89, 0.45] .121
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² │ 0.02
#> Adj. R² │ 0.02
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> BMI: │
#> Outcome │ bmi
#> (Intercept) │ 23.98 0.40 [23.20, 24.75] <.001
#> age │ 0.04 0.01 [ 0.02, 0.05] <.001
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 0.51 0.22 [ 0.09, 0.94] .018
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ -0.06 0.26 [-0.58, 0.46] .822
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1163
#> R² │ 0.02
#> Adj. R² │ 0.02
#>
#> Note. Linear regression models. Std. errors: classical (OLS). CI = confidence interval.table_stack() stacks tables built separately, each with
its own arguments. A panel whose columns differ from the first one opens
with its own header, here the odds ratios of a logistic model:
m_dentist <- glm(dentist_12m ~ age + sex + smoking,
data = sochealth, family = binomial)
table_stack(
BMI = table_regression(m_bmi),
Dentist = table_regression(m_dentist, exponentiate = TRUE)
)
#> Variable │ B SE 95% CI p
#> ─────────────────┼──────────────────────────────────────
#> BMI: │
#> (Intercept) │ 23.98 0.40 [23.20, 24.75] <.001
#> age │ 0.04 0.01 [ 0.02, 0.05] <.001
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 0.51 0.22 [ 0.09, 0.94] .018
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ -0.06 0.26 [-0.58, 0.46] .822
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1163
#> R² │ 0.02
#> Adj. R² │ 0.02
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> Dentist: │
#> │ OR SE 95% CI p
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> (Intercept) │ 1.66 0.39 [ 1.05, 2.62] .030
#> age │ 1.01 0.00 [ 1.00, 1.02] .064
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 1.00 0.13 [ 0.77, 1.28] .975
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ 0.95 0.15 [ 0.70, 1.30] .767
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Tjur) │ 0.00
#> R² (McFadden) │ 0.00
#> AIC │ 1421.8
#>
#> Note. BMI: linear regression; Dentist: logistic regression. Dentist: coefficients exponentiated and displayed as OR; SE on the OR scale (delta method); CI bounds exponentiated (asymmetric). Dentist: the (Intercept) row estimates the baseline odds. Std. errors: BMI: classical (OLS), Dentist: classical (Fisher information). CI = confidence interval; AIC = Akaike information criterion.The notes of the panels are merged as the note of the same models
side by side: one sentence gives the type of each model, one its
standard errors, and a sentence that holds for one panel only opens with
its label. A table without a name takes its label from
panel_labels, else from a title you typed for it, else from
its outcome. table_stack() adds no “Outcome” row of its
own: a table built with outcome_labels
(table_regression(m_bmi, outcome_labels = "bmi")) brings
one.
Univariable screening
Epidemiological reporting often starts one step before the
side-by-side table: a univariable screen — one simple
regression per candidate predictor — set against the multivariable
model, the classic “Table 2” of clinical and epidemiological journals
(the layout of gtsummary::tbl_uvregression() and the
EpiRHandbook workflow). table_regression_uv() builds that
entire layout in one call: pass the data, the outcome, and the candidate
predictors — reusing sh, the
unordered-education copy of sochealth
introduced in Per-coefficient effect sizes:
table_regression_uv(sh, outcome = smoking,
predictors = c(sex, age, education),
family = binomial(), exponentiate = TRUE)
#> Univariable and multivariable logistic regression: smoking
#>
#> Univariable Multivariable
#> ─────────────────────────────── ─────────────
#> Variable │ N OR 95% CI p OR
#> ──────────────────────────┼────────────────────────────────────────────────
#> sex: │
#> Female (ref.) │ – – – – –
#> Male │ 1175 0.95 [0.72, 1.26] .713 0.96
#> age │ 1175 1.01 [1.00, 1.01] .292 1.01
#> education: │
#> Lower secondary (ref.) │ – – – – –
#> Upper secondary │ 1175 0.62 [0.44, 0.87] .005 0.62
#> Tertiary │ 0.41 [0.28, 0.60] <.001 0.40
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Tjur) │ 0.02
#> R² (McFadden) │ 0.02
#> AIC │ 1200.9
#>
#> Multivariable
#> ────────────────────
#> Variable │ 95% CI (OR) p (OR)
#> ──────────────────────────┼──────────────────────
#> sex: │
#> Female (ref.) │ – –
#> Male │ [0.73, 1.28] .800
#> age │ [1.00, 1.02] .214
#> education: │
#> Lower secondary (ref.) │ – –
#> Upper secondary │ [0.44, 0.87] .005
#> Tertiary │ [0.27, 0.59] <.001
#>
#> Note. Logistic regression models. Coefficients exponentiated and displayed as OR; CI bounds exponentiated. Std. errors: classical (Fisher information). CI = confidence interval; AIC = Akaike information criterion.Reading it: each Univariable row is its own simple logistic
regression — the N column is per predictor, because
listwise deletion can differ across candidates — while the
Multivariable column group is the single joint model on the
same candidates (multivariable = FALSE drops it,
complete_cases = TRUE aligns every fit on the common
complete cases). Comparing the two columns per row is the point of the
layout: an association that survives adjustment reads differently from
one that collapses. The multivariable n appears once, as a
fit statistic — it belongs to the joint model, not to any predictor. A
p_adjust argument treats the screen as one family of tests.
For a time-to-event screen the same call takes
method = "coxph" with a Surv() outcome — see
Survival regression
tables.
Hierarchical / nested regression
Set nested = TRUE to add in-table
change-statistic rows (APA Table 7.13 / Stata
esttab / SPSS Model Summary convention). Each adjacent pair
(M2 vs M1, M3 vs M2, …) contributes one column of change stats below
R² / Adj. R²; the first model column gets en-dashes (no
previous model to compare to).
Hierarchical regression requires every model in the sequence to share
the same analytic sample: identical nobs and an identical
response variable. R’s listwise deletion otherwise produces different
n per model as soon as a predictor introduced by M2 or M3
has missing values that M1 did not see — a silent bias on the change
statistics. We restrict to complete cases on the union of predictors up
front:
sochealth_cc <- sochealth |>
dplyr::select(wellbeing_score, age, sex, smoking, bmi, region) |>
na.omit()
m1 <- lm(wellbeing_score ~ age + sex, data = sochealth_cc)
m2 <- lm(wellbeing_score ~ age + sex + smoking, data = sochealth_cc)
m3 <- lm(wellbeing_score ~ age + sex + smoking + bmi, data = sochealth_cc)
table_regression(list(m1, m2, m3), nested = TRUE)
#> Hierarchical linear regression: wellbeing_score
#>
#> Model 1 Model 2
#> ──────────────────── ──────────────────────────────
#> Variable │ B SE p B SE p
#> ─────────────────┼──────────────────────────────────────────────────────
#> (Intercept) │ 64.70 1.66 <.001 65.00 1.67 <.001
#> age │ 0.05 0.03 .118 0.05 0.03 .109
#> sex: │
#> Female (ref.) │ – – – – – –
#> Male │ 3.89 0.91 <.001 3.88 0.91 <.001
#> smoking: │
#> No (ref.) │ – – –
#> Yes │ -1.68 1.11 .132
#> bmi │
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1163 1163
#> R² │ 0.02 0.02
#> Adj. R² │ 0.02 0.02
#> ΔR² │ – +0.00
#> F-change │ – 2.28 (1, 1159)
#> p (change) │ – .132
#>
#> Model 3
#> ──────────────────────────────
#> Variable │ B SE p
#> ─────────────────┼────────────────────────────────
#> (Intercept) │ 80.57 3.37 <.001
#> age │ 0.07 0.03 .019
#> sex: │
#> Female (ref.) │ – – –
#> Male │ 4.21 0.90 <.001
#> smoking: │
#> No (ref.) │ – – –
#> Yes │ -1.71 1.10 .119
#> bmi │ -0.65 0.12 <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1163
#> R² │ 0.04
#> Adj. R² │ 0.04
#> ΔR² │ +0.02
#> F-change │ 28.13 (1, 1158)
#> p (change) │ <.001
#>
#> Note. Linear regression models. Std. errors: classical (OLS).Reading the change rows: each test carries its degrees of freedom,
2.28 (1, 1159) for an F on 1 and 1159 df (the added term
and the residual df of the larger model), as Stata prints F(1, 1159).
Adding smoking in Model 2 does not improve fit over the age
+ sex baseline (ΔR² = +0.00, F-change = 2.28, p = .132) — with a single
added term, this test is equivalent to the test on the
smoking coefficient itself (p = .132 in the Model 2
column). Adding bmi in Model 3 does (ΔR² = +0.02, F-change
= 28.13, p < .001). The hierarchical sequence therefore supports
retaining bmi and provides no incremental evidence for
smoking.
Default change tokens auto-injected:
c("r2_change", "f_change", "p_change") for lm
(APA hierarchical-regression standard),
c("lrt_change", "p_change") for glm (Hosmer
& Lemeshow §3.5), c("f_change", "p_change") for a
glm whose dispersion is estimated (the F test
anova(test = "F") reports for gaussian, Gamma, inverse
Gaussian, Tweedie and quasi families),
c("aic_change", "bic_change", "lrt_change", "p_change") for
mixed-effects fits (lmer / glmer /
glmmTMB / nlme::lme — Pinheiro & Bates
2000 §2.4.1). Customize via show_fit_stats; the order of
tokens controls the order of the rows. Other change tokens are
available: "adj_r2_change", "f2_change",
"deviance_change", "aic_change" /
"aicc_change" / "bic_change".
Variance-explained change tokens (Δr², Δf²) are undefined for
mixed-effects pairs, and asking for one is a
spicy_invalid_input error naming "lrt_change"
+ "p_change" instead: the F-test framework that grounds
them doesn’t apply. The same refusal covers every hierarchy whose nested
comparison is a likelihood-ratio test (survreg,
polr, clm, coxph,
gls, betareg, multinom,
fixest), while a quantile hierarchy is pointed at
"f_change" + "p_change".
Validation is strict: identical nobs and identical
response across all models, otherwise a spicy_invalid_input
error explains the listwise-deletion trap and suggests refitting on the
common subset.
The change test follows vcov. Under any estimator other
than the classical one, each pair is tested by a Wald test of the
coefficients the current model adds, on that model’s robust matrix, so
the block test and the coefficients agree. The row then reads Wald
F-change (also for a glm whose dispersion is
estimated; Wald χ² (change) for the other glm
fits, coxph and the other likelihood classes), and ΔR² is
unchanged. A Wald test of a large block in a small sample rejects too
often, so read it with care when many predictors enter at once.
table_regression(list(m1, m2, m3), nested = TRUE, vcov = "HC3")
#> Hierarchical linear regression: wellbeing_score
#>
#> Model 1 Model 2
#> ──────────────────── ──────────────────────────────
#> Variable │ B SE p B SE p
#> ─────────────────┼──────────────────────────────────────────────────────
#> (Intercept) │ 64.70 1.60 <.001 65.00 1.62 <.001
#> age │ 0.05 0.03 .114 0.05 0.03 .106
#> sex: │
#> Female (ref.) │ – – – – – –
#> Male │ 3.89 0.91 <.001 3.88 0.91 <.001
#> smoking: │
#> No (ref.) │ – – –
#> Yes │ -1.68 1.12 .135
#> bmi │
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1163 1163
#> R² │ 0.02 0.02
#> Adj. R² │ 0.02 0.02
#> ΔR² │ – +0.00
#> Wald F-change │ – 2.24 (1, 1159)
#> p (change) │ – .135
#>
#> Model 3
#> ──────────────────────────────
#> Variable │ B SE p
#> ─────────────────┼────────────────────────────────
#> (Intercept) │ 80.57 3.27 <.001
#> age │ 0.07 0.03 .018
#> sex: │
#> Female (ref.) │ – – –
#> Male │ 4.21 0.91 <.001
#> smoking: │
#> No (ref.) │ – – –
#> Yes │ -1.71 1.11 .122
#> bmi │ -0.65 0.12 <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1163
#> R² │ 0.04
#> Adj. R² │ 0.04
#> ΔR² │ +0.02
#> Wald F-change │ 29.34 (1, 1158)
#> p (change) │ <.001
#>
#> Note. Linear regression models. Std. errors: heteroskedasticity-robust (HC3).Comparing models by their fit statistics
show_columns = "none" keeps the fit statistics and drops
every coefficient. Each model has one column, and the rows follow
show_fit_stats. With nested = TRUE the change
rows come last, so the table reads as the summary of the hierarchy.
table_regression(
list(m1, m2, m3),
show_columns = "none",
nested = TRUE,
show_fit_stats = c("nobs", "r2", "adj_r2", "aic", "bic",
"r2_change", "f_change", "p_change")
)
#> Hierarchical linear regression: wellbeing_score
#>
#> Variable │ Model 1 Model 2 Model 3
#> ────────────┼───────────────────────────────────────────────
#> n │ 1163 1163 1163
#> R² │ 0.02 0.02 0.04
#> Adj. R² │ 0.02 0.02 0.04
#> AIC │ 9683.9 9683.6 9657.7
#> BIC │ 9704.2 9708.9 9688.1
#> ΔR² │ – +0.00 +0.02
#> F-change │ – 2.28 (1, 1159) 28.13 (1, 1158)
#> p (change) │ – .132 <.001
#>
#> Note. Linear regression models. AIC = Akaike information criterion; BIC = Bayesian information criterion.The title still names the outcome, and the note keeps what describes
the models. Arguments that act on coefficients only, such as
exponentiate or keep, are ignored with a
warning. Without nested = TRUE, a table that compares an
information criterion across models fitted to different samples or
outcomes warns that the difference is not interpretable.
broom::glance() returns one row per model with its
model-level statistics, nobs, r.squared,
adj.r.squared, AIC and BIC among
them. The change rows compare two models, so they are in the table
only.
Robust variance
The default vcov = "classical" reports the OLS standard
error, valid under homoskedastic, independent errors. Two robust
alternatives are first-class arguments; bootstrap and jackknife variants
are also available via vcov = "bootstrap" /
"jackknife".
Heteroskedasticity-consistent (HC)
When error variance plausibly depends on the predictors (a ubiquitous
concern in cross-sectional social-science data), set
vcov = "HC*" for sandwich-style standard errors via
sandwich::vcovHC(). The valid types are HC0
through HC5; HC3 is the small-sample-friendly
default (Long and Ervin 2000):
fit <- lm(wellbeing_score ~ age + sex + smoking, data = sochealth)
table_regression(fit, vcov = "HC3")
#> Linear regression: wellbeing_score
#>
#> Variable │ B SE 95% CI p
#> ─────────────────┼──────────────────────────────────────
#> (Intercept) │ 65.20 1.61 [62.05, 68.35] <.001
#> age │ 0.05 0.03 [-0.01, 0.11] .127
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 3.86 0.91 [ 2.07, 5.64] <.001
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ -1.72 1.11 [-3.91, 0.47] .123
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² │ 0.02
#> Adj. R² │ 0.02
#>
#> Note. Linear regression. Std. errors: heteroskedasticity-robust (HC3). CI = confidence interval.The footer reads “Std. errors: heteroskedasticity-robust
(HC3)”; the column header for the confidence interval automatically
tracks ci_level.
Cluster-robust (CR)
For clustered observations (repeated measures on the same person,
students within schools, observations within regions),
vcov = "CR*" requests cluster-robust variance via
clubSandwich::vcovCR(), with the cluster identifier
supplied through cluster. Three forms are accepted:
-
Formula —
cluster = ~region. The variables are looked up inmodel.frame(fit)first, then in the model’s originaldataargument. Recommended: independent of the dataset’s name, composable for multi-way clustering (cluster = ~region:year). -
String —
cluster = "region". Single column name resolved the same way. Convenient but cannot express interactions. -
Vector —
cluster = df$region. Atomic vector of lengthnobs(fit). Use this when the cluster key is derived on the fly (cluster = interaction(df$region, df$year)) or pulled from a different dataset with matching row order.
Bare unquoted names (cluster = region) are
not accepted — they would require non-standard
evaluation that breaks under programmatic use (function wrapping, loops,
dynamic column choice).
fit <- lm(wellbeing_score ~ age + sex + smoking, data = sochealth_cc)
table_regression(
fit,
vcov = "CR2",
cluster = ~region,
show_columns = c("b", "se", "ci", "p", "ame", "ame_p")
)
#> Linear regression: wellbeing_score
#>
#> Variable │ B SE 95% CI p AME p
#> ─────────────────┼────────────────────────────────────────────────────
#> (Intercept) │ 65.00 1.74 [60.49, 69.51] <.001
#> age │ 0.05 0.04 [-0.05, 0.15] .247 0.05 .247
#> sex: │
#> Female (ref.) │ – – – – – –
#> Male │ 3.88 0.85 [ 1.68, 6.07] .006 3.88 .006
#> smoking: │
#> No (ref.) │ – – – – – –
#> Yes │ -1.68 1.55 [-5.72, 2.37] .331 -1.68 .331
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1163
#> R² │ 0.02
#> Adj. R² │ 0.02
#>
#> Note. Linear regression. Std. errors: cluster-robust (CR2), clusters by region. AME inference: t-test with Satterthwaite df. AME of a continuous predictor: average slope (dY/dX). CI = confidence interval; AME = average marginal effect.CR2 (the Bell-McCaffrey adjustment) is the recommended
default under few clusters (Pustejovsky and Tipton 2018; Imbens and
Kolesár 2016). Coefficient inference uses
clubSandwich::coef_test() with Satterthwaite-corrected
degrees of freedom. When AME columns are requested, the same
Satterthwaite framework is applied to the AME contrast via
clubSandwich::linear_contrast(), so B and AME
share the same t-distribution with the same df — visible in the
footer.
Multiple-comparison adjustment
p_adjust applies a family-wise or false-discovery-rate
adjustment to the p-values via stats::p.adjust(). The
family is the model’s full coefficient set (intercept and reference rows
excluded); B and AME are adjusted as
independent families within the same call. Available methods are the
same as in stats::p.adjust(): "none"
(default), "holm", "hochberg",
"hommel", "bonferroni", "BH" /
"fdr", "BY".
fit <- lm(wellbeing_score ~ age + sex + smoking + education,
data = sh)
table_regression(fit, p_adjust = "bonferroni")
#> Linear regression: wellbeing_score
#>
#> Variable │ B SE 95% CI p
#> ──────────────────────────┼──────────────────────────────────────
#> (Intercept) │ 53.96 1.66 [50.70, 57.22] <.001
#> age │ 0.03 0.03 [-0.03, 0.08] 1.000
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 3.57 0.81 [ 1.99, 5.15] <.001
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ 0.68 0.99 [-1.27, 2.63] 1.000
#> education: │
#> Lower secondary (ref.) │ – – – –
#> Upper secondary │ 11.89 1.05 [ 9.82, 13.96] <.001
#> Tertiary │ 19.71 1.12 [17.51, 21.91] <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² │ 0.22
#> Adj. R² │ 0.22
#>
#> Note. Linear regression. Std. errors: classical (OLS). P-values adjusted for 5 tests (Bonferroni method). CI = confidence interval.The footer documents the chosen method and the family size; the SE
column is unchanged. p_adjust is orthogonal to
vcov: combining a robust SE with a family-wise correction
is fully supported, e.g. vcov = "HC3", p_adjust = "holm".
The adjustment is applied before keep /
drop filtering, so the family stays the model’s full
coefficient set and the displayed adjusted p-values reflect the right
denominator regardless of which subset is shown.
A methodological note. Adjusting the p-values of every coefficient of
a single regression model is not the standard convention in
social-science or clinical reporting (Rothman 1990; Gelman, Hill, and
Yajima 2012; Greenland 2017; Harrell 2015 §5.4; APA Manual 7 §6.46).
Each coefficient tests a scientifically distinct hypothesis on a
distinct predictor, which is not the situation that family-wise
procedures were designed for. The default "none" reflects
this consensus.
Adjustment is nonetheless legitimate in three contexts:
-
Mass screening with many candidate predictors and
no prior hypothesis — typically
"BH"(Benjamini–Hochberg, false discovery rate). -
Pre-registered multi-endpoint confirmatory designs
— typically
"holm"(Holm’s step-down, strong family-wise error rate control). - When a reviewer, an editor, or a statistical analysis plan explicitly requests it — apply the requested method and document it in the footer.
The argument follows the same transparency over rejection
rule as standardized: the method is available, the choice
among methods is the analyst’s, and the footer makes that choice visible
to the reader.
Filtering displayed coefficients
keep and drop accept regular expressions
matched against coefficient names (as returned by
stats::coef()). They are mutually exclusive — use
keep to retain the focal predictors, or drop
to hide a few control variables. Multiple patterns combine with logical
OR. The intercept is not filtered: show_intercept decides
whether it shows.
fit <- lm(wellbeing_score ~ age + sex + smoking + bmi + education,
data = sh)
table_regression(fit, keep = c("^smoking", "^bmi$"))
#> Linear regression: wellbeing_score
#>
#> Variable │ B SE 95% CI p
#> ─────────────┼──────────────────────────────────────
#> (Intercept) │ 51.41 3.52 [44.50, 58.32] <.001
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ 0.79 1.00 [-1.17, 2.75] .428
#> bmi │ 0.10 0.12 [-0.14, 0.33] .418
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1163
#> R² │ 0.23
#> Adj. R² │ 0.22
#>
#> Note. Linear regression. Std. errors: classical (OLS). CI = confidence interval.
table_regression(fit, drop = "^education")
#> Linear regression: wellbeing_score
#>
#> Variable │ B SE 95% CI p
#> ─────────────────┼──────────────────────────────────────
#> (Intercept) │ 51.41 3.52 [44.50, 58.32] <.001
#> age │ 0.02 0.03 [-0.03, 0.08] .407
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 3.50 0.81 [ 1.91, 5.10] <.001
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ 0.79 1.00 [-1.17, 2.75] .428
#> bmi │ 0.10 0.12 [-0.14, 0.33] .418
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1163
#> R² │ 0.23
#> Adj. R² │ 0.22
#>
#> Note. Linear regression. Std. errors: classical (OLS). CI = confidence interval.Together with p_adjust, this is the article-ready
workflow: adjust on the full model, display only the rows the reader
cares about —
table_regression(fit, p_adjust = "BH", keep = "^treatment").
Significance stars
Stars are off by default. APA 7 §6.46 explicitly discourages
asterisks-only reporting in favor of exact p-values, and ASA’s post-2019
guidance (Wasserstein, Schirm, and Lazar 2019) reinforces the same
point. Set stars = TRUE for the APA preset
(*** p < .001, ** p < .01, * p < .05) or pass a
named numeric vector for custom thresholds:
fit <- lm(wellbeing_score ~ age + sex + smoking, data = sochealth)
table_regression(fit, stars = TRUE)
#> Linear regression: wellbeing_score
#>
#> Variable │ B SE 95% CI p
#> ─────────────────┼─────────────────────────────────────────
#> (Intercept) │ 65.20*** 1.66 [61.95, 68.45] <.001
#> age │ 0.05 0.03 [-0.01, 0.11] .130
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 3.86*** 0.91 [ 2.08, 5.63] <.001
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ -1.72 1.11 [-3.89, 0.45] .121
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² │ 0.02
#> Adj. R² │ 0.02
#>
#> Note. Linear regression. Std. errors: classical (OLS). CI = confidence interval.
#> *** p < .001, ** p < .01, * p < .05.Stars suffix the B column (or β when
standardization is requested); the threshold mapping is auto-documented
in the footer. The p column itself remains unstarred so the
numeric value stays readable.
Display options
labels names the variables in four forms:
FALSE (the default) keeps the names the model gives them,
and TRUE shows the label attribute of the outcome, the
numeric terms, and the factor headers. A named character vector keyed by
formula term labels ("sex") or coefficient names
("sexMale") sets the wording you publish, and
var_labels(age = "Age in years") renames a few variables
while the others keep their label attribute.
intercept_position switches the intercept to the bottom
(Stata convention). reference_style = "annotation" lifts
the reference level into the factor header
("sex: [ref: Female]") and drops the explicit reference
row, for compact tables. decimal_mark = "," switches to
European decimal notation and automatically changes the CI separator to
";" to avoid ambiguity:
fit <- lm(wellbeing_score ~ age + sex + education, data = sh)
table_regression(
fit,
labels = c(
age = "Age (years)",
sex = "Sex",
education = "Education"
),
reference_style = "annotation",
intercept_position = "last",
decimal_mark = ","
)
#> Linear regression: wellbeing_score
#>
#> Variable │ B SE 95% CI p
#> ───────────────────────────────────┼───────────────────────────────────────
#> Age (years) │ 0,03 0,03 [-0,03; 0,08] 0,343
#> Sex: [ref: Female] │
#> Male │ 3,65 0,80 [ 2,09; 5,22] <0,001
#> Education: [ref: Lower secondary] │
#> Upper secondary │ 11,85 1,04 [ 9,81; 13,89] <0,001
#> Tertiary │ 19,52 1,10 [17,36; 21,67] <0,001
#> (Intercept) │ 54,18 1,63 [50,99; 57,37] <0,001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1200
#> R² │ 0,22
#> Adj. R² │ 0,22
#>
#> Note. Linear regression. Std. errors: classical (OLS). CI = confidence interval.Generalized linear models (glm)
table_regression() accepts any glm() fit.
Its Wald test is the one summary.glm() reports:
t on the residual degrees of freedom when the dispersion is
estimated (gaussian, Gamma, inverse.gaussian and the quasi- families),
z when the family fixes it (binomial, Poisson). The
confidence interval and the average marginal effects use the same
distribution. The table title becomes family-aware (“Logistic
regression”, “Poisson regression”, “Probit regression”, …) and the
default fit-statistics block swaps in nobs, the R² rows of
the family and aic instead of R² and
Adj. R². A logit on a binary outcome shows Tjur’s and
McFadden’s R² (pseudo_r2_tjur,
pseudo_r2_mcfadden), the other binomial, Poisson and
negative binomial models McFadden’s alone, and the gaussian, Gamma,
inverse Gaussian and Tweedie families the deviance R²
(r2_deviance).
We illustrate with a logistic regression of smoking on
sex, age, and education from
sochealth, reusing sh — the copy of
sochealth with education as a plain
(unordered) factor, introduced in Per-coefficient effect sizes
— so contrasts are dummy-coded, one row per level. The model mixes a
binary factor (sex), a numeric predictor
(age), and a 3-level factor (education) —
enough to exercise every glm-specific feature below.
fit <- glm(smoking ~ sex + age + education, data = sh,
family = binomial)
table_regression(fit)
#> Logistic regression: smoking
#>
#> Variable │ B SE 95% CI p
#> ──────────────────────────┼──────────────────────────────────────
#> (Intercept) │ -1.11 0.29 [-1.67, -0.55] <.001
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ -0.04 0.14 [-0.32, 0.25] .800
#> age │ 0.01 0.00 [-0.00, 0.02] .214
#> education: │
#> Lower secondary (ref.) │ – – – –
#> Upper secondary │ -0.48 0.17 [-0.82, -0.14] .005
#> Tertiary │ -0.91 0.20 [-1.29, -0.52] <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Tjur) │ 0.02
#> R² (McFadden) │ 0.02
#> AIC │ 1200.9
#>
#> Note. Logistic regression. Std. errors: classical (Fisher information). CI = confidence interval; AIC = Akaike information criterion.Read the fit block with care. Neither row is a share of explained
variance, and neither compares with an OLS R². Tjur’s (2009) R² is the
mean predicted probability of the smokers minus that of the non-smokers.
McFadden’s (1974) is the proportional reduction of the log-likelihood
from the intercept-only model. Both read 0.02 here. Other measures give
other numbers on the same fit: Nagelkerke’s R², which
show_fit_stats adds with
"pseudo_r2_nagelkerke", is 0.03. When reporting one, name
it, as the row labels do: a generic “pseudo-R²” is ambiguous across
software (Long & Freese 2014, p. 126). A threshold quoted for one
measure does not carry over to another (Hemmert et al. 2018). Long and
Freese (2014, p. 120) find these measures of limited use in their own
research, with the exception of the information criteria, hence AIC in
the default block. The reference page, ?table_regression,
gives the measure of each family and why.
Response-scale display: exponentiate = TRUE
Set exponentiate = TRUE to switch the B
column to the ratio scale and relabel its header according to the family
and link: OR for binomial(logit),
IRR for poisson(log), HR for
binomial(cloglog), RR for
binomial(log), MR for Gamma(log),
and the generic exp(B) for other log-link families (a
genuine ratio of means). Links whose exponential is not a ratio
— probit, cauchit, inverse (the Gamma() default), among
others — are refused with a clear error rather than silently mislabeled;
report response-scale effects for those models via the AME column
instead. The standard error follows the delta-method approximation
SE_OR = OR × SE_log-odds (the Stata logit, or
convention). The test statistic and the p-value stay on the link scale,
where B = 0 and OR = 1 are the same hypothesis
and the Wald approximation is most accurate, so they match the
unexponentiated table verbatim. (Wald tests are not invariant to
reparameterization — a z rebuilt from SE_OR would give a
different statistic and p-value, which is why the link-scale test is
kept; Long & Freese 2014 §3.2.2.):
table_regression(fit, exponentiate = TRUE)
#> Logistic regression: smoking
#>
#> Variable │ OR SE 95% CI p
#> ──────────────────────────┼────────────────────────────────────
#> (Intercept) │ 0.33 0.09 [0.19, 0.58] <.001
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 0.96 0.14 [0.73, 1.28] .800
#> age │ 1.01 0.00 [1.00, 1.02] .214
#> education: │
#> Lower secondary (ref.) │ – – – –
#> Upper secondary │ 0.62 0.11 [0.44, 0.87] .005
#> Tertiary │ 0.40 0.08 [0.27, 0.59] <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Tjur) │ 0.02
#> R² (McFadden) │ 0.02
#> AIC │ 1200.9
#>
#> Note. Logistic regression. Coefficients exponentiated and displayed as OR; SE on the OR scale (delta method); CI bounds exponentiated (asymmetric). The (Intercept) row estimates the baseline odds. Std. errors: classical (Fisher information). CI = confidence interval; AIC = Akaike information criterion.Read exponentiated rows with the factor-change template (Long and
Freese 2014 §6.1.1): for Tertiary, OR = 0.40, 95% CI [0.27,
0.59] — 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, a 60% decrease in the odds (percentage change = 100
× (OR − 1)). An odds ratio is a ratio of odds, not of
probabilities: with a smoking prevalence around 21% it sits farther from
1 than the corresponding risk ratio, so reading it as “a 60% lower
probability of smoking” overstates the effect. The Average marginal
effects section below reports the probability-scale counterpart for
the same contrast (−0.15, that is, 15 percentage points).
The (Intercept) row is not a ratio, and the note says
what it is: the baseline odds, the odds of smoking when every factor
sits at its reference level and every numeric predictor at zero. Here
that is a woman with lower secondary education aged 0, younger than
anyone in the sample (25 to 75), so the 0.33 anchors the model and
describes no respondent. Centering age on its mean would
make it the odds of such a woman at the mean age.
Event counts per level
Clinical and epidemiological tables conventionally carry the raw
events / N next to each odds ratio (STROBE checklist
item 15; the NEJM house style). The opt-in "n_events"
column token renders exactly that — per factor level, with the
full-sample count on numeric predictors and the intercept row:
table_regression(fit, show_columns = c("n_events", "b", "ci", "p"))
#> Logistic regression: smoking
#>
#> Variable │ Events/N B 95% CI p
#> ──────────────────────────┼─────────────────────────────────────────────
#> (Intercept) │ 249/1175 -1.11 [-1.67, -0.55] <.001
#> sex: │
#> Female (ref.) │ 131/606 – – –
#> Male │ 118/569 -0.04 [-0.32, 0.25] .800
#> age │ 249/1175 0.01 [-0.00, 0.02] .214
#> education: │
#> Lower secondary (ref.) │ 78/257 – – –
#> Upper secondary │ 112/527 -0.48 [-0.82, -0.14] .005
#> Tertiary │ 59/391 -0.91 [-1.29, -0.52] <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Tjur) │ 0.02
#> R² (McFadden) │ 0.02
#> AIC │ 1200.9
#>
#> Note. Logistic regression. Std. errors: classical (Fisher information). CI = confidence interval; AIC = Akaike information criterion.The token is off by default and deliberately narrow: it is
defined for binomial models with an ungrouped 0/1 response and for
right-censored Cox models (events per factor level; Cox tables also
carry N events as a default fit statistic — see Survival regression
tables). Requesting it anywhere else is a hard, classed error
rather than a silently empty column.
Term-level partial chi-square
partial_chi2 is the glm analog of
partial_f2: for each model term, the partial
likelihood-ratio chi-square from the same Type-II nested comparison —
both models exclude any higher-order term containing the tested one, so
the test respects marginality and is contrast-invariant (Fox and
Weisberg 2019; equal to drop1(test = "LRT") for additive
models; on the LRT itself see Long & Freese 2014 §3.2.2, §3.2.4).
Rendered as value (df) so factor terms (k − 1
df) and numeric terms (1 df) read at a glance:
fit2 <- glm(smoking ~ age + education, data = sh, family = binomial)
table_regression(fit2, show_columns = c("b", "partial_chi2", "p"))
#> Logistic regression: smoking
#>
#> Variable │ B χ² p
#> ──────────────────────────┼───────────────────────────
#> (Intercept) │ -1.13 <.001
#> age │ 0.01 1.55 (1) .213
#> education: │
#> Lower secondary (ref.) │ – –
#> Upper secondary │ -0.48 21.70 (2) .005
#> Tertiary │ -0.91 21.70 (2) <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Tjur) │ 0.02
#> R² (McFadden) │ 0.02
#> AIC │ 1198.9
#>
#> Note. Logistic regression. Std. errors: classical (Fisher information). χ² = partial likelihood-ratio chi-squared; AIC = Akaike information criterion.The p column and the χ² column answer
different questions. Each per-level p-value compares one education level
with the reference, so it changes when the reference level changes; the
joint test — 21.70 on 2 df, p < .001 — asks whether
education as a whole improves the model, and is invariant
to the reference choice. When a factor’s levels straddle significance,
the joint test is the one to report. Wald and likelihood-ratio tests of
the same joint hypothesis are asymptotically equivalent but differ in
finite samples; many statisticians prefer the likelihood-ratio form when
both are available (Long & Freese 2014 §3.2.2), which is what
partial_chi2 computes.
The term-level test of a glm is the F that
anova(reduced, full, test = "F") reports when the
dispersion is estimated, the likelihood-ratio chi-square otherwise. A
gaussian, Gamma, inverse Gaussian, Tweedie or quasi-family
glm therefore asks for its twin, partial_f,
the column an lm takes too, and refuses
partial_chi2. Its value is the one
car::Anova(type = 2, test.statistic = "F", error.estimate = "dispersion")
reports (car’s default Pearson estimate differs by the IRLS convergence
error only), and its cell gives both of its degrees of freedom,
3.45 (2, 28) for an F on 2 and 28, as Stata and SAS report
them. A gaussian glm and its lm() twin print
the same F. A table that mixes both kinds of model asks for both tokens,
and each model fills its own column.
The same token extends to mixed-effects models (lmer,
glmer, glmmTMB, nlme::lme), where
the term-level test is the Wald chi-square built from the same
marginality-respecting Type-II hypothesis, matching
car::Anova(type = 2). This is a deliberate departure from
SAS PROC MIXED and lmerTest, which report Type-III tests by
default; the two conventions coincide for models without
interactions.
Standardized coefficients: refit and
pseudo
For glm, standardized = "refit" z-scores
numeric predictors only and refits the model — the response
stays on its observed scale because the link function is fixed. This is
the “x-standardization” convention (Long and Freese 2014 §4.7.2). The
other algebraic methods ("posthoc", "basic",
"smart") apply X-only scaling using the same algebra as in
the lm case.
standardized = "pseudo" (glm only) is the
latent-scale variant of Menard’s (2004, 2011) fully standardized
coefficient: Y* is the latent variable on the link scale,
with SD(Y*) = sqrt(var(linear-predictor) + var_link) and
var_link = π²/3 for logit, 1 for probit, π²/6 for cloglog.
Numeric predictors scale by SD(X) / SD(Y*); factor dummies
stay on their 0/1 scale and divide by SD(Y*) only — the
same dummy convention as the other methods, which is why the factor rows
in the table below are not B × SD(X) / SD(Y*). Defined for
binomial families; non-binomial returns NA with a
spicy_caveat:
table_regression(fit, standardized = "pseudo")
#> Logistic regression: smoking
#>
#> Variable │ B β SE 95% CI p
#> ──────────────────────────┼─────────────────────────────────────────────
#> (Intercept) │ -1.11 – 0.29 [-1.67, -0.55] <.001
#> sex: │
#> Female (ref.) │ – – – – –
#> Male │ -0.04 -0.02 0.14 [-0.32, 0.25] .800
#> age │ 0.01 0.05 0.00 [-0.00, 0.02] .214
#> education: │
#> Lower secondary (ref.) │ – – – – –
#> Upper secondary │ -0.48 -0.26 0.17 [-0.82, -0.14] .005
#> Tertiary │ -0.91 -0.49 0.20 [-1.29, -0.52] <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Tjur) │ 0.02
#> R² (McFadden) │ 0.02
#> AIC │ 1200.9
#>
#> Note. Logistic regression. Std. errors: classical (Fisher information). β = standardized coefficient ("pseudo": latent-scale SD(X)/SD(Y*) for numeric predictors, 1/SD(Y*) for factor dummies); CI = confidence interval; AIC = Akaike information criterion.Average marginal effects: probability units, not log-odds
For glm, AME is averaged over the observed sample — by
default the partial derivative for numeric predictors, the discrete
change for factor levels — on the response scale, not the link
scale. For logistic regression this means the displayed AME is
in probability units, not log-odds: an AME of
-0.15 for education = Tertiary reads “on
average, holding sex and age constant, a respondent with tertiary
education is 15 percentage points less likely to be a current smoker
than one with lower secondary education”. This is almost always the
quantity a substantive reader wants from a logistic regression, and it
is the reason AME is increasingly recommended over odds ratios for
communicating logistic effects (Mood 2010; Long and Freese 2014
§6.2.3).
table_regression(fit, show_columns = c("b", "p", "ame", "ame_ci", "ame_p"))
#> Logistic regression: smoking
#>
#> Variable │ B p AME 95% CI p
#> ──────────────────────────┼──────────────────────────────────────────────
#> (Intercept) │ -1.11 <.001
#> sex: │
#> Female (ref.) │ – – – – –
#> Male │ -0.04 .800 -0.01 [-0.05, 0.04] .800
#> age │ 0.01 .214 0.00 [-0.00, 0.00] .214
#> education: │
#> Lower secondary (ref.) │ – – – – –
#> Upper secondary │ -0.48 .005 -0.09 [-0.16, -0.02] .007
#> Tertiary │ -0.91 <.001 -0.15 [-0.22, -0.09] <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Tjur) │ 0.02
#> R² (McFadden) │ 0.02
#> AIC │ 1200.9
#>
#> Note. Logistic regression. Std. errors: classical (Fisher information). AME of a continuous predictor: average slope (dY/dX). AME = average marginal effect; CI = confidence interval; AIC = Akaike information criterion.Point estimate and inference are delegated to
marginaleffects::avg_slopes(), or to
avg_comparisons() under the ame_type of the
next section. Under cluster-robust variance, the Satterthwaite-df
handling described in the Robust variance section applies to
glm AME as well.
Slope or discrete change: ame_type
For a continuous predictor, ame_type sets what the AME
columns report:
-
"slope", the default: the average slope dY/dX, the derivative of the prediction at each observation’s own value, averaged over the sample; -
"unit": the average change in the prediction when the predictor rises by one unit from each observation’s own value, from x to x + 1; -
"sd": the same change for a rise of one standard deviation, from x to x + SD, with SD the unweighted standard deviation of the predictor in the estimation sample.
Factors, logical and 0/1 predictors stay contrasts with their reference level under all three, and the note names the quantity of the continuous rows. In a linear model without transformed terms the change for one unit is the slope, and the change for one SD is the slope times the SD. Under a non-identity link the prediction curves between x and x + 1. The slope then approximates the change for one unit only as far as the curve is straight over that step, and the gap grows as it bends (Long and Freese 2014, §6.2, p. 241).
The probability of rating one’s health good or very good, modeled on
monthly household income in thousands of francs, shows the gap. Each
call below asks for one quantity, and table_stack() sets
the three one under the other:
sh <- sh |>
dplyr::mutate(
good_health = self_rated_health %in% c("Good", "Very good"),
income_k = income / 1000
)
fit_inc <- glm(good_health ~ income_k + age + sex, data = sh,
family = binomial)
ame_income <- function(type) {
table_regression(fit_inc, show_columns = c("ame", "ame_ci"),
ame_type = type, keep = "income_k",
show_intercept = FALSE, show_fit_stats = FALSE,
digits = 3)
}
table_stack(
Slope = ame_income("slope"),
Unit = ame_income("unit"),
SD = ame_income("sd")
)
#> Logistic regression: good_health
#>
#> Variable │ AME 95% CI
#> ──────────┼───────────────────────
#> Slope: │
#> income_k │ 0.080 [0.062, 0.097]
#> ╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> Unit: │
#> income_k │ 0.073 [0.058, 0.088]
#> ╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> SD: │
#> income_k │ 0.098 [0.079, 0.118]
#>
#> Note. Logistic regression models. Std. errors: classical (Fisher information). Slope: AME of a continuous predictor: average slope (dY/dX). Unit: AME of a continuous predictor: average change in the prediction for a one-unit increase. SD: AME of a continuous predictor: average change in the prediction for a one-standard-deviation increase. AME = average marginal effect; CI = confidence interval.The slope says the probability rises at a rate of 0.080 per thousand
francs. One more thousand francs a month raises it by 0.073 on average:
71% of the respondents already rate their health good or very good, the
logistic curve flattens as the probability approaches 1, and the step
gains less than the tangent at its start promises. A sentence such as
“one more thousand francs a month goes with a probability 7.3 points
higher” describes "unit", not the slope. A rise of one
standard deviation, 1,395 francs, raises the probability by 0.098, less
than the 0.111 that the slope times the SD would give. Report
"unit" when the text speaks of one more unit,
"sd" when the natural unit is too small or too large to
read (Long and Freese 2014, p. 494), and the slope where the field reads
marginal changes, as economics does (p. 241).
A discrete change predicts at values the data may not hold. Incomes in the sample run from 1,000 to 7,388 francs, and the step carries some respondents past the highest of them:
c(unit = sum(sh$income_k + 1 > max(sh$income_k)),
sd = sum(sh$income_k + sd(sh$income_k) > max(sh$income_k)))
#> unit sd
#> 41 84For these 41 respondents, 84 under "sd", the changed
income lies beyond every income observed, and their part of the average
is an extrapolation of the logistic curve. Gelman and Pardoe (2007,
p. 32) keep predictive comparisons within the support of the data for
that reason. The share is small here. On a bounded scale it need not be:
on a rating from 1 to 5, everyone at 5 is moved to an impossible 6. The
slope, a derivative at the observed value, does not leave the data.
The "sd" step starts from each observation’s own value:
the uncentered change of Long and Freese (2014, p. 168), and the
discrete change of Mize, Doan, and Long (2019, §3.2), for whom “one or a
standard deviation are common values” of the step. It is not the
"sd" shortcut of
marginaleffects::avg_comparisons(), which moves every
observation from the mean minus half a standard deviation to the mean
plus half, whatever its own value: 0.114 on this model, against
0.098.
Profile-likelihood CIs: ci_method = "profile"
The default is ci_method = "wald" (symmetric
estimate ± q × SE, with q the quantile of the
coefficient’s test, t or z as summary.glm() reports it).
For glm you can also request
ci_method = "profile", which calls
MASS::confint.glm() to compute the profile-likelihood CI —
asymmetric, obtained by inverting the likelihood-ratio statistic rather
than the quadratic Wald approximation, and the preferred choice under
sparse data or near-boundary estimates (Venables & Ripley
MASS §7.2). The estimate, SE, statistic, and p-value all remain
Wald; "profile" only refines the CI bounds:
table_regression(fit, ci_method = "profile")
#> Logistic regression: smoking
#>
#> Variable │ B SE 95% CI p
#> ──────────────────────────┼──────────────────────────────────────
#> (Intercept) │ -1.11 0.29 [-1.68, -0.55] <.001
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ -0.04 0.14 [-0.32, 0.25] .800
#> age │ 0.01 0.00 [-0.00, 0.02] .214
#> education: │
#> Lower secondary (ref.) │ – – – –
#> Upper secondary │ -0.48 0.17 [-0.82, -0.14] .005
#> Tertiary │ -0.91 0.20 [-1.29, -0.52] <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Tjur) │ 0.02
#> R² (McFadden) │ 0.02
#> AIC │ 1200.9
#>
#> Note. Logistic regression. Std. errors: classical (Fisher information). 95% CIs: profile likelihood. CI = confidence interval; AIC = Akaike information criterion.One consequence worth naming when combining with
exponentiate = TRUE: the CI bounds are the exponentiated
profile bounds (asymmetric), while the SE column is the
delta-method SE on the ratio scale — the interval can no longer be
reconstructed as estimate ± q × SE. Both choices are
disclosed in the footer, which is why the table names them rather than
leaving the mismatch for a careful reader to discover.
Hierarchical glm (LRT)
nested = TRUE on a list of nested glm fits appends
Δχ² and p (change) rows to the fit-statistics
block — default change tokens c("lrt_change", "p_change"),
as in Hierarchical / nested regression above — the
APA-conventional layout for hierarchical logistic regression (Hosmer
& Lemeshow §3.5; Long & Freese 2014 §3.2.4). The chi-square
statistic comes from anova(test = "LRT"). The
variance-explained change tokens (r2_change,
f_change) are not defined for glm — requesting them in
show_fit_stats raises a spicy_invalid_input
error that points to the pseudo-R² family and to lrt_change
+ p_change instead:
m1 <- glm(smoking ~ sex, data = sh, family = binomial)
m2 <- glm(smoking ~ sex + age, data = sh, family = binomial)
m3 <- glm(smoking ~ sex + age + education, data = sh, family = binomial)
table_regression(list(m1, m2, m3), nested = TRUE)
#> Hierarchical logistic regression: smoking
#>
#> Model 1 Model 2
#> ──────────────────── ────────────────────────
#> Variable │ B SE p B SE p
#> ──────────────────────────┼────────────────────────────────────────────────
#> (Intercept) │ -1.29 0.10 <.001 -1.54 0.26 <.001
#> sex: │
#> Female (ref.) │ – – – – – –
#> Male │ -0.05 0.14 .713 -0.05 0.14 .722
#> age │ 0.01 0.00 .294
#> education: │
#> Lower secondary (ref.) │
#> Upper secondary │
#> Tertiary │
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175 1175
#> R² (Tjur) │ 0.00 0.00
#> R² (McFadden) │ 0.00 0.00
#> AIC │ 1217.6 1218.5
#> Δχ² │ – 1.10 (1)
#> p (change) │ – .294
#>
#> Model 3
#> ────────────────────────
#> Variable │ B SE p
#> ──────────────────────────┼──────────────────────────
#> (Intercept) │ -1.11 0.29 <.001
#> sex: │
#> Female (ref.) │ – – –
#> Male │ -0.04 0.14 .800
#> age │ 0.01 0.00 .214
#> education: │
#> Lower secondary (ref.) │ – – –
#> Upper secondary │ -0.48 0.17 .005
#> Tertiary │ -0.91 0.20 <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² (Tjur) │ 0.02
#> R² (McFadden) │ 0.02
#> AIC │ 1200.9
#> Δχ² │ 21.64 (2)
#> p (change) │ <.001
#>
#> Note. Logistic regression models. Std. errors: classical (Fisher information). AIC = Akaike information criterion.Gaussian glm caveat
A glm() with family = gaussian and
link = "identity" is mathematically equivalent to
lm() but lacks the variance- explained effect-size family
(partial_f2 / η² / ω²) and the Satterthwaite-corrected AME
path. Following the transparency over rejection rule, spicy
accepts the fit and emits a spicy_caveat suggesting a refit
with lm().
Mixed-effects models
table_regression() supports the mixed-effects fits of
lme4::lmer() and lme4::glmer(),
glmmTMB::glmmTMB(), and nlme::lme(), and, for
an ordinal outcome, of ordinal::clmm() and glmmTMB’s
ordinal family (see Ordinal regression tables). The
fixed-effects section of the table follows the same conventions as
lm / glm. Below it, a subordinate
Random effects block reports the variance components as
rows — one per random standard deviation (σ), correlation (ρ), and the
residual — each with its estimate, SE, and CI in the shared coefficient
columns. The number of groups — one
N (<grouping factor>) row per grouping factor, here
N (Subject) — and the Nakagawa marginal / conditional R²
(Nakagawa & Schielzeth 2013) join the fit-statistic rows, so they
align per model in a comparison table:
library(lme4)
fit <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
table_regression(fit)
#> Linear mixed-effects regression: Reaction
#>
#> Variable │ B SE 95% CI p
#> ─────────────────────────────────┼────────────────────────────────────────
#> (Intercept) │ 251.41 6.82 [238.03, 264.78] <.001
#> Days │ 10.47 1.55 [ 7.44, 13.50] <.001
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> Random effects: │
#> σ Subject (Intercept) │ 24.74 5.84 [ 6.79, 34.32] –
#> σ Subject Days │ 5.92 1.25 [ 2.47, 8.00] –
#> ρ Subject ((Intercept), Days) │ 0.07 0.33 [ -0.57, 0.70] –
#> σ (Residual) │ 25.59 1.51 [ 22.44, 28.39] –
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 180
#> N (Subject) │ 18
#> R² (marginal) │ 0.28
#> R² (conditional) │ 0.80
#> AIC │ 1755.6
#> BIC │ 1774.8
#>
#> Note. Linear mixed-effects regression. Std. errors: Wald (model-based). P-values: Wald z-test, large-sample approximation. Random effects (REML): LR test vs linear regression, χ̄²(3) = 150.04, p < .001. CI = confidence interval; AIC = Akaike information criterion; BIC = Bayesian information criterion.No ICC row appears here: a random-slope fit has no
single intra-class correlation, so the row is dropped automatically (the
intercept-only case shows it).
Two deliberate conventions:
- The variance-component rows carry no p-value. The
null hypothesis σ = 0 sits on the boundary of the parameter space, where
a Wald test is invalid (Self & Liang 1987); no reporting guideline
requests a per-component p, and the R ecosystem (
lme4,nlme,parameters,broom.mixed), Stata, and MLwiN all omit it. The correct significance signal for the random part is the footer’s likelihood-ratio test against the no-random-effects model, whose p-value applies a boundary correction. With a single variance component the halvingp = 0.5 × P(χ²₁ > LR)is the exact Self & Liang (1987) 50:50 mixture — thechibar2(01)line Stata’smixedprints. With more random parameters, as here (two variances and a correlation), the exact reference is a multi-component mixture (Stram & Lee 1994); spicy reports the pragmatic halvedχ²_q, less conservative than the plainχ²_qthat Stata prints for joint tests with a caveat that it is conservative. - The footer also carries the estimator label —
(REML)for restricted maximum likelihood (thelmer/lmedefault) or(ML)for full ML (glmer,glmmTMBdefault, andlmer(..., REML = FALSE)). The label is informational: it tells the reader which likelihood the variance estimates are anchored to, without implying that REML and ML estimates are interchangeable for inference.
The display and inference controls are arguments, not conventions:
show_re = FALSE hides the block,
re_scale = "variance" switches the rows from σ to σ²,
re_columns restricts which quantities the variance rows
show, re_ci = "profile" replaces the Wald CIs on variance
components with boundary-respecting profile-likelihood CIs (Wald is the
default because profiling is far slower), and
re_test = "lrt" or "rlrt" adds an opt-in
per-component test (“do we need this random slope?”) with the boundary
correction applied. The footer names the fixed-effect inference method
per engine — Satterthwaite t when lmerTest is loaded,
Wald-z for glmer / glmmTMB and plain
lmer, containment-df t for nlme::lme. The AME
tokens of show_columns and ame_type (see
Average marginal effects above) work identically for
mixed-effects fits — always on the response scale, with Wald-z
inference.
nested = TRUE works on mixed-effects fits with the
change tokens listed in Hierarchical / nested regression; note
that its plain-χ² reference is appropriate for steps that add fixed
effects, while random-structure steps belong to re_test,
whose reference respects the boundary.
For the model-building sequence, the REML/ML estimator choice, the
Nakagawa R² and ICC mechanics, per-class inference paths, singular fits,
testing random components, and the glmmTMB /
nlme::lme specifics — including a worked glmer
AME example — see Mixed-effects regression
tables.
Output formats
The default print method draws an ASCII table to the console. The
same content is available as a plain data.frame with its
numbers (output = "data.frame": one row per printed row,
the interval in two bounds, the key columns variable,
label, level and .row_role in
front), a long broom-style data.frame (output = "long"), or
as a rich-format table for gt, tinytable,
flextable, excel, word, or
clipboard:
out <- table_regression(
lm(wellbeing_score ~ age + sex + smoking, data = sochealth),
output = "data.frame"
)
str(out)
#> 'data.frame': 11 obs. of 9 variables:
#> $ variable : chr "(Intercept)" "age" "sex" "sex" ...
#> $ label : chr "(Intercept)" "age" "sex" "sex" ...
#> $ level : chr NA NA NA "Female" ...
#> $ .row_role: chr "coef" "coef" "factor_header" "reference" ...
#> $ B : num 65.2009 0.0465 NA NA 3.8558 ...
#> $ SE : num 1.6567 0.0307 NA NA 0.9053 ...
#> $ 95% CI LL: num 61.9504 -0.0137 NA NA 2.0796 ...
#> $ 95% CI UL: num 68.451 0.107 NA NA 5.632 ...
#> $ p : num 1.59e-216 1.30e-01 NA NA 2.22e-05 ...
#> - attr(*, "note")= chr "Note. Linear regression. Std. errors: classical (OLS). CI = confidence interval."
out <- table_regression(
lm(wellbeing_score ~ age + sex + smoking, data = sochealth),
output = "long"
)
out[, c("model_id", "term", "estimate_type", "estimate",
"std.error", "p.value")]
#> # A tibble: 6 × 6
#> model_id term estimate_type estimate std.error p.value
#> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 M1 (Intercept) B 65.2 1.66 1.59e-216
#> 2 M1 age B 0.0465 0.0307 1.30e- 1
#> 3 M1 sexFemale B NA NA NA
#> 4 M1 sexMale B 3.86 0.905 2.22e- 5
#> 5 M1 smokingNo B NA NA NA
#> 6 M1 smokingYes B -1.72 1.11 1.21e- 1
pkgdown_dark_gt(
table_regression(
lm(wellbeing_score ~ age + sex + smoking, data = sochealth),
output = "gt"
)
)| Linear regression: wellbeing_score | |||||
|
Variable
|
B
|
SE
|
95% CI
|
p
|
|
|---|---|---|---|---|---|
| LL | UL | ||||
| (Intercept) | 65.20 | 1.66 | 61.95 | 68.45 | <.001 |
| age | 0.05 | 0.03 | −0.01 | 0.11 | .130 |
| sex: | |||||
| Female (ref.) | – | – | – | – | – |
| Male | 3.86 | 0.91 | 2.08 | 5.63 | <.001 |
| smoking: | |||||
| No (ref.) | – | – | – | – | – |
| Yes | −1.72 | 1.11 | −3.89 | 0.45 | .121 |
| n | 1175 | ||||
| R² | 0.02 | ||||
| Adj. R² | 0.02 | ||||
broom integration
broom::tidy() returns a long tibble with one row per
(model_id, term, estimate_type) and broom-canonical column
names (estimate, std.error,
conf.low, conf.high, statistic,
p.value). broom::glance() returns one row per
(model_id, outcome) with model-level statistics;
df.residual is preserved as a numeric so cluster-robust
Satterthwaite df flows through verbatim:
fit <- lm(wellbeing_score ~ age + sex + smoking, data = sochealth)
out <- table_regression(fit, standardized = "refit")
broom::tidy(out)
#> # A tibble: 8 × 17
#> model_id outcome outcome_level term estimate_type estimate std.error conf.low
#> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 M1 wellbe… NA (Int… B 65.2 1.66 62.0
#> 2 M1 wellbe… NA (Int… beta -0.0961 0.0431 -0.181
#> 3 M1 wellbe… NA age B 0.0465 0.0307 -0.0137
#> 4 M1 wellbe… NA age beta 0.0439 0.0290 -0.0130
#> 5 M1 wellbe… NA sexM… B 3.86 0.905 2.08
#> 6 M1 wellbe… NA sexM… beta 0.247 0.0579 0.133
#> 7 M1 wellbe… NA smok… B -1.72 1.11 -3.89
#> 8 M1 wellbe… NA smok… beta -0.110 0.0708 -0.249
#> # ℹ 9 more variables: conf.high <dbl>, statistic <dbl>, df <dbl>,
#> # df.residual <dbl>, p.value <dbl>, test_type <chr>, is_intercept <lgl>,
#> # factor_term <chr>, factor_level <chr>
broom::glance(out)
#> # A tibble: 1 × 15
#> model_id outcome nobs weighted_nobs r.squared adj.r.squared omega2 sigma
#> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 M1 wellbeing_s… 1175 NA 0.0190 0.0165 0.0165 15.5
#> # ℹ 7 more variables: rmse <dbl>, f2 <dbl>, AIC <dbl>, AICc <dbl>, BIC <dbl>,
#> # deviance <dbl>, df.residual <dbl>The long format is the right entry point when the table is one step
in a larger pipeline — saving to disk for the manuscript appendix,
faceting by subgroup, or feeding a downstream post-estimation analysis.
The tidy() output keeps the estimate_type
column so the same data frame can hold rows for B, β, AME, and
per-coefficient effect-size estimates without ambiguity.
See also
- Supported models — the map of every supported class: estimands, variance estimators, standardized-coefficient scope, and fit statistics per family.
- Categorical predictors for categorical predictors across all model families: dummy coding and reference levels, joint tests, ordinal predictors (scores vs dummies), and the categorization trap.
Model-family articles build on the mechanics shown here:
- Mixed-effects regression tables for mixed-effects (multilevel) models: random effects as table rows, ICC, the boundary-corrected chi-bar-squared test, and opt-in per-term LRT / RLRT.
- GEE regression tables for population-averaged (GEE) models: the fit’s own sandwich inference, working-correlation choice via QIC, and jackknife variants for few clusters.
- Count and two-part regression tables for count and two-part models: Poisson and negative-binomial rate ratios, offsets, and zero-inflated / hurdle components as labeled blocks.
- Ordinal regression tables for ordinal (proportional-odds) models: shared slopes, ordered thresholds as rows, and per-category marginal effects.
- Multinomial regression tables for multinomial logistic models: outcome categories as columns and changing the reference category.
- Survival regression tables for survival models: Cox hazard ratios with events and concordance, accelerated failure time models, and survival estimands.
- Bayesian regression tables for Bayesian fits: posterior medians and credible intervals, the probability of direction, and random effects from the draws.
Descriptive and reporting companions:
- Categorical summary tables for the APA Table 1 categorical descriptors (factors, labelled variables, chi-squared or Fisher’s exact tests chosen per variable, association measures).
- Continuous summary tables for the APA Table 1 / Table 2 continuous descriptors and unadjusted group-comparison tests.
- Model-based continuous summary tables for the one-predictor-by-many-outcomes linear-model counterpart (estimated marginal means, contrast or slope, four effect-size families with noncentral CIs).
- Summary tables for reporting for an end-to-end reporting workflow that combines the spicy table functions along the APA Table 1 / 2 / 3 sequence.
References
Aiken, L. S., and West, S. G. (1991). Multiple Regression: Testing and Interpreting Interactions. Sage.
American Psychological Association (2020). Publication Manual of the American Psychological Association (7th ed.). Section 6.46.
Cohen, J., Cohen, P., West, S. G., and Aiken, L. S. (2003). Applied Multiple Regression / Correlation Analysis for the Behavioral Sciences (3rd ed.). Lawrence Erlbaum.
Fox, J. (2016). Applied Regression Analysis and Generalized Linear Models (3rd ed.). Sage.
Fox, J., and Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). Sage.
Gelman, A. (2008). Scaling regression inputs by dividing by two standard deviations. Statistics in Medicine, 27(15), 2865–2873.
Gelman, A., and Pardoe, I. (2007). Average predictive comparisons for models with nonlinearity, interactions, and variance components. Sociological Methodology, 37(1), 23–51.
Gelman, A., Hill, J., and Yajima, M. (2012). Why we (usually) don’t have to worry about multiple comparisons. Journal of Research on Educational Effectiveness, 5(2), 189–211.
Greenland, S. (2017). For and against methodologies: Some perspectives on recent causal and statistical inference debates. European Journal of Epidemiology, 32(1), 3–20.
Harrell, F. E. (2015). Regression Modeling Strategies (2nd ed.). Springer. Section 5.4 on multiplicity.
Hemmert, G. A. J., Schons, L. M., Wieseke, J., and Schimmelpfennig, H. (2018). Log-likelihood-based pseudo-R² in logistic regression: Deriving sample-sensitive benchmarks. Sociological Methods & Research, 47(3), 507–531.
Hosmer, D. W., Lemeshow, S., and Sturdivant, R. X. (2013). Applied Logistic Regression (3rd ed.). Wiley.
Imbens, G. W., and Kolesár, M. (2016). Robust standard errors in small samples: Some practical advice. Review of Economics and Statistics, 98(4), 701–712.
Long, J. S., and Ervin, L. H. (2000). Using heteroscedasticity consistent standard errors in the linear regression model. The American Statistician, 54(3), 217–224.
Long, J. S., and Freese, J. (2014). Regression Models for Categorical Dependent Variables Using Stata (3rd ed.). Stata Press.
McFadden, D. (1974). Conditional logit analysis of qualitative choice behavior. In P. Zarembka (Ed.), Frontiers in Econometrics (pp. 105–142). Academic Press.
Menard, S. (2004). Six approaches to calculating standardized logistic regression coefficients. The American Statistician, 58(3), 218–223.
Menard, S. (2011). Standards for standardized logistic regression coefficients. Social Forces, 89(4), 1409–1428.
Mize, T. D., Doan, L., and Long, J. S. (2019). A general framework for comparing predictions and marginal effects across models. Sociological Methodology, 49(1), 152–189.
Mood, C. (2010). Logistic regression: Why we cannot do what we think we can do, and what we can do about it. European Sociological Review, 26(1), 67–82.
Nagelkerke, N. J. D. (1991). A note on a general definition of the coefficient of determination. Biometrika, 78(3), 691–692.
Nakagawa, S., and Schielzeth, H. (2013). A general and simple method for obtaining R² from generalized linear mixed-effects models. Methods in Ecology and Evolution, 4(2), 133–142.
Olejnik, S., and Algina, J. (2003). Generalized eta and omega squared statistics: Measures of effect size for some common research designs. Psychological Methods, 8(4), 434–447.
Pinheiro, J. C., and Bates, D. M. (2000). Mixed-Effects Models in S and S-PLUS. Springer.
Pustejovsky, J. E., and Tipton, E. (2018). Small-sample methods for cluster-robust variance estimation and hypothesis testing in fixed effects models. Journal of Business & Economic Statistics, 36(4), 672–683.
Rothman, K. J. (1990). No adjustments are needed for multiple comparisons. Epidemiology, 1(1), 43–46.
Self, S. G., and Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under non-standard conditions. Journal of the American Statistical Association, 82(398), 605–610.
Smithson, M. (2003). Confidence Intervals. Sage.
Steiger, J. H. (2004). Beyond the F test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis. Psychological Methods, 9(2), 164–182.
Stram, D. O., and Lee, J. W. (1994). Variance components testing in the longitudinal mixed effects model. Biometrics, 50(4), 1171–1177.
Tjur, T. (2009). Coefficients of determination in logistic regression models — A new proposal: The coefficient of discrimination. The American Statistician, 63(4), 366–372.
Venables, W. N., and Ripley, B. D. (2002). Modern Applied Statistics with S (4th ed.). Springer. Section 7.2 on profile likelihood.
Wasserstein, R. L., Schirm, A. L., and Lazar, N. A. (2019). Moving to a world beyond “p < 0.05”. The American Statistician, 73(sup1), 1–19.