Regression coefficient summary table

Publication-ready coefficient table from one or more fitted lm / glm models. Supports standardized coefficients ($\beta$), average marginal effects (AME), partial effect sizes ($f^2$ / $\eta^2$ / $\omega^2$ and partial F for lm; partial $\chi^2$ or partial F for glm), pseudo-$R^2$ (glm), and a full vocabulary of variance estimators (classical / HC* / cluster-robust with Satterthwaite-corrected df / bootstrap / jackknife). glm covers binomial / poisson / Gamma / inverse.gaussian / quasi families with any link.

Usage

table_regression(
  models,
  vcov = "classical",
  cluster = NULL,
  ci_level = 0.95,
  ci_method = c("wald", "profile", "boot_percentile", "hdi"),
  boot_n = 1000L,
  tau = NULL,
  at_time = NULL,
  standardized = c("none", "refit", "posthoc", "basic", "smart", "pseudo"),
  exponentiate = FALSE,
  ame_type = c("slope", "unit", "sd"),
  p_adjust = "none",
  show_columns = NULL,
  keep = NULL,
  drop = NULL,
  show_intercept = TRUE,
  show_thresholds = TRUE,
  show_components = TRUE,
  intercept_position = c("first", "last"),
  factor_layout = c("grouped", "flat"),
  reference_style = c("row", "annotation", "footer", "none"),
  reference_label = NULL,
  show_fit_stats = NULL,
  fit_stats_layout = c("first_col", "merged"),
  layout = c("columns", "stacked"),
  show_re = TRUE,
  re_scale = c("sd", "variance"),
  re_columns = c("est", "se", "ci"),
  re_test = c("none", "lrt", "rlrt"),
  re_ci = c("wald", "profile"),
  model_labels = NULL,
  outcome_labels = NULL,
  stars = FALSE,
  nested = FALSE,
  digits = 2L,
  p_digits = 3L,
  effect_size_digits = 2L,
  fit_digits = 2L,
  ic_digits = 1L,
  decimal_mark = ".",
  align = c("decimal", "center", "right"),
  padding = 0L,
  labels = getOption("spicy.var_labels", FALSE),
  note_labels = NULL,
  title = NULL,
  note = NULL,
  note_add = NULL,
  output = c("default", "data.frame", "long", "tinytable", "gt", "flextable", "excel",
    "clipboard", "word"),
  excel_path = NULL,
  excel_sheet = NULL,
  clipboard_delim = "\t",
  word_path = NULL,
  word_template = NULL,
  style = NULL
)

Arguments

models

A fitted model object, or a list of such fits (named or unnamed; classes may be mixed). Single fits are auto-promoted to a 1-element list. A broad set of model classes is supported – linear / generalized linear (lm, glm, MASS::glm.nb), mixed-effects (lmer, lme, glmmTMB), survival (coxph, survreg), ordinal (polr, clm, and the cumulative-link mixed models ordinal::clmm and glmmTMB with family = glmmTMB::ordinal()), mgcv::gam/bam, betareg, mlogit, survey::svyglm, rms (ols/lrm/cph/Glm), and Bayesian (rstanarm/brms), among others. An unsupported class raises spicy_unsupported. Raw data + formula is not accepted – fit-only API.

vcov

Variance-covariance estimator: "classical", "HC0"-"HC5", "CR0"-"CR3", "bootstrap", or "jackknife". A scalar is recycled to all models; a list (one string per model) allows mixed estimators. Default "classical". The resampling estimators (lm / glm, including MASS::glm.nb) refit each replicate on resampled rows of the fixed evaluated design; for glm.nb the dispersion parameter theta is held at its full-sample estimate (Stata's nbreg, vce(bootstrap) re-estimates it per replicate – the two conventions differ slightly). Resampling that fails on nearly every replicate raises spicy_resampling_failed rather than silently reporting a different estimator. Quantile regression (quantreg::rq) uses its own estimator family instead: "classical" resolves to the heteroskedasticity-robust "nid" sandwich (Hendricks-Koenker, Hall-Sheather bandwidth – quantreg's own large-sample default and the Stata vce(robust) analog), with "iid" (Koenker-Bassett homoskedastic, Stata vce(iid) parity), "ker" (Powell kernel) and "rank" (rank-score inversion: genuine CIs, no SE / t / p – those cells render as dashes) as opt-ins, and "bootstrap" running quantreg's native boot.rq (bsmethod = "xy"; with cluster =, the Hagemann 2017 wild gradient cluster bootstrap – the one cluster-robust route for rq; CR* and HC* are refused for rq, as is "jackknife", whose leave-one-out form is inconsistent for quantiles). With nested = TRUE, the change test of each pair uses the same estimator (a Wald test; see Hierarchical (nested) model comparison). See Inference and standard errors.

cluster

Cluster identifier for cluster-robust variance (used when vcov is "CR0"-"CR3" or a cluster-bootstrap / cluster-jackknife). Three accepted forms (see How to specify cluster in the details):

Formula: ~region, ~region:year (recommended).
String column name: "region".
Atomic vector of length nobs(fit): df$region, interaction(df$region, df$year), ... (for keys derived on the fly).

For multi-model use, pass a list of one form per model (mix-and-match allowed). A bare unquoted name (cluster = region) is not column selection: it is evaluated as an ordinary variable (the vector form when region exists in the calling environment, a migration error pointing at ~region / "region" otherwise). Default NULL (no clustering).

ci_level

Confidence level for all reported CIs (B, $\beta$, AME, partial effect sizes). Default 0.95.

ci_method

CI construction. "wald" (default) uses estimate +/- z x SE (t x SE for lm). "profile" (glm and ordinal MASS::polr / ordinal::clm) uses the profile-likelihood CI from stats::confint() – MASS::confint.glm() for glm, confint.polr() / confint.clm() for ordinal fits – asymmetric, exact for likelihood-based inference (Venables & Ripley MASS Section 7.2). Only the CI bounds change; estimate, SE, statistic and p-value remain Wald. For ordinal fits the profile covers the predictor coefficients (the cut-point thresholds stay Wald), and a robust vcov takes precedence (its Wald-robust CIs are used instead, with a consolidated warning). "profile" with lm raises spicy_invalid_input. "boot_percentile" (requires vcov = "bootstrap" for every model) replaces the coefficient CI bounds with equal-tailed percentile intervals of the bootstrap replicates (the boot::boot.ci() type = "perc" convention; Davison & Hinkley 1997, ch. 5), reusing the same resamples as the bootstrap SEs. Only the CI bounds change; estimate, SE, statistic and p-value remain Wald from the bootstrap covariance – the Stata convention (normal-based table CIs by default, percentile on request via estat bootstrap). AME and standardized-beta CIs are not covered. With exponentiate = TRUE the percentile bounds are exponentiated (percentile intervals are transformation-respecting). "hdi" (Bayesian stanreg / brmsfit fits only) replaces the default equal-tailed credible interval with the highest-density interval – the shortest interval containing ci_level of the posterior draws (Kruschke 2015), relabeling the column header 95% HDI. Unlike the equal-tailed interval the HDI is not transformation-invariant, so under exponentiate = TRUE it is recomputed on the exponentiated draws rather than transformed. Requesting "hdi" for a frequentist fit, or "profile" / "boot_percentile" for a Bayesian fit, raises spicy_invalid_input.

boot_n

Number of bootstrap replicates when vcov = "bootstrap", and of the bootstrap of the RMST and risk-difference columns: a positive whole number, as in table_continuous_lm() and table_multiresponse(). Default 1000L, the default of table_continuous_lm(), which computes the same bootstrap variance. That variance, of a lm or glm fit, needs at least 10 valid replicates: a smaller boot_n is an error. Typed in a table that runs no bootstrap, it is ignored with a warning (class spicy_ignored_arg).

tau

RMST horizon: the "rmst" column family reports the restricted-mean-survival-time difference over [0, tau] (coxph and survreg fits). A positive time on the outcome's scale, or "minmax" for the smallest per-group maximum follow-up (the resolved value is disclosed in the table note). No default: the horizon defines the estimand. Stratified fits (strata()) are supported: standardization keeps each subject's own stratum baseline, disclosed in the table note.

at_time

Landmark time for the "risk_diff" column family: the difference in cumulative incidence at at_time (coxph and survreg fits). A positive time on the outcome's scale; no default.

standardized

Standardization method for the "beta" column. One of "none" (default), "refit", "posthoc", "basic", "smart", "pseudo". "pseudo" is glm only (Menard 2011 fully-standardized); using it with lm() raises spicy_invalid_input. Supported classes: lm, glm (incl. MASS::glm.nb), the mixed engines (lmer / glmer / glmmTMB / nlme::lme), and fixed-effects stan_glm-style stanreg fits for the algebraic flavors "posthoc" / "basic" / "smart" – exact affine rescales of the posterior draws (median, MAD SD and credible bounds all scale by the same positive SD ratio; the beta stays on the link scale under exponentiate = TRUE). Standard-formula fixed-effects brmsfit models are supported too (their design matrix is recovered through insight); the scale factors are engine-invariant, so the same model fit with rstanarm or brms gets identical betas up to sampling noise. Bayesian "refit" (would re-run the sampler) and "pseudo" (per-draw latent variance; planned) are refused. Also refused: multilevel fits (stan_glmer / stan_lmer and brm() with group-level terms; a Bayesian beta would need an explicit sd(Y) decomposition that the frequentist mixed engines resolve by refitting on z-scored data), non-GLM stanreg subclasses (stan_polr, stan_betareg), and brms formulas with distributional, multivariate or special terms (mo(), s(), ...). In all refused Bayesian cases, standardize predictors before fitting instead. Any other class raises spicy_unsupported_standardized rather than rendering an empty beta column. See the Standardized coefficients section.

exponentiate

Logical. When TRUE, B, the CI bounds, and the SE (delta method: SE_OR = OR x SE_log-odds) are exp()-transformed for links where the result is a ratio: logit (OR), log (IRR / RR / MR per family, generic exp(B) for other log-link families – a genuine ratio of means), and binomial / ordinal cloglog (HR; grouped-time proportional hazards, Prentice & Gloeckler 1978). For the ordinal (cumulative) cloglog model the displayed HR is exp(-B), not exp(B): the cumulative parametrization cloglog P(Y <= j) = zeta_j - xB places the hazard of the grouped event on -B, so exp(B) would be the reciprocal of the hazard ratio (the CI endpoints are negated and swapped accordingly, and the table note discloses the convention; a cloglog clm with nominal = terms refuses exponentiate because threshold-side coefficients have no covariate-HR reading). The statistic and p-value stay on the link scale (invariant under monotone transformation). Identity-link fits are left untouched (a spicy_ignored_arg warning fires when no model in the table exponentiates), so mixed lm + logit tables keep working. Any other link (probit, cauchit, inverse – the Gamma() default –, 1/mu^2, sqrt, ordinal loglog, ...) raises spicy_invalid_input: the exponential of such a coefficient has no ratio interpretation, and reporting it would mislabel the estimate. Report response-scale effects for those models via the AME column instead. Note that the displayed CI is the link-scale CI with exponentiated endpoints (Wald or profile per ci_method), so it is asymmetric around the ratio and cannot be reconstructed as estimate ± q × SE; the delta-method SE follows the Stata convention ([R] logistic, Methods and formulas: se(OR) = OR × se(b)) and the log-scale CI endpoints carry the uncertainty more faithfully than this SE. The footer states the SE scale and the CI asymmetry whenever an SE column is displayed. It also names what a displayed intercept estimates once exponentiated, which is no ratio: the baseline odds (logit), risk (binomial log), rate (count families, log) or mean (other log-link families) at the reference profile. Default FALSE.

ame_type

The quantity the AME columns report for a continuous predictor. "slope" (default): the average slope dY/dX, the derivative of the prediction averaged over the observations (Stata's margins, dydx()). "unit": the average change in the prediction when the predictor rises by one unit from each observation's own value. "sd": the same for a rise of one standard deviation, the unweighted standard deviation of the predictor in the estimation sample, again from each observation's own value. These are the uncentered unit and standard-deviation changes of Long and Freese (2014), averaged over the observations as in Mize, Doan and Long (2019). The step is taken on the variable as it is in the data, before any transformation the formula applies (log(x), poly(x, 2)). For the observations with the highest values, it can carry the predictor beyond the values observed: the prediction there is an extrapolation of the model. Gelman and Pardoe (2007) consider changes only in the region of support of the data. "sd" is not the "sd" shortcut of marginaleffects::avg_comparisons(), which moves every observation between the same two values, the mean minus and the mean plus half a standard deviation. For an unweighted fit and marginaleffects' default forward contrast, it equals avg_comparisons(fit, variables = list(x = sd(x))), the standard deviation taken over the rows the fit used. The three agree when the prediction is linear in the predictor (an lm without transformed terms: "unit" equals the slope, "sd" the slope times the standard deviation). Under a non-identity link they differ, the more so as the coefficient grows: with a log link and no interaction, the slope is b * mean(mu) and the change for one unit (exp(b) - 1) * mean(mu). "unit" is the quantity the sentence "one more unit of X goes with ... more of Y" describes; the slope only approximates it. Categorical, logical and 0/1 predictors are contrasts with their reference level under every value. Weighted fits average with the same weights under every value, per-category models (ordinal, multinomial) report the quantity for each category probability, and Bayesian fits summarize its draws. The table note names the quantity. Ignored, with a warning, when no AME column is shown.

p_adjust

Multiple-comparison adjustment method applied to the family of estimated coefficient p-values within each model (intercept and reference rows excluded). One of "none" (default), "holm", "hochberg", "hommel", "bonferroni", "BH" / "fdr", or "BY". Delegated to stats::p.adjust(); applied per-model and per estimate_type (B and AME p-values are adjusted independently within their own families). Active adjustments are documented in the footer (method + family size). See the Multiple-comparison adjustment section for when this is and is not appropriate.

show_columns

Character vector of tokens selecting the per-coefficient columns and their display order. Accepts atomic tokens ("b", "se", "ci", "t", "p", "beta", "n", "n_events", "pd" (probability of direction, Bayesian fits only), "rhat" / "ess_bulk" / "ess_tail" / "mcse" (per-coefficient sampler diagnostics and the Monte Carlo standard error of the displayed posterior median, all-Bayesian tables only), "ame", "ame_se", "ame_ci", "ame_p", "rmst" + "rmst_se" / "rmst_ci" / "rmst_p", "risk_diff" + its _se / _ci / _p companions, "partial_f2" + "partial_f2_ci", "partial_eta2" + "partial_eta2_ci", "partial_omega2" + "partial_omega2_ci", "partial_chi2", "partial_f" (lm, and glm whose dispersion is estimated), "r2" / "adj_r2" (per-fit variance explained, lm screens)) and group tokens ("all_b", "all_b_compact", "all_b_full", "all_beta", "all_ame", "all_ame_compact", "all_f2", "all_eta2", "all_omega2"). See Vocabulary tokens in the details for the full enumeration. Default NULL selects a context-aware layout: "all_b" (single model) or "all_b_compact" (multi-model, and the single-multinomial outcome-as-columns layout, which has the same width pressure – restore CIs with atomic tokens, e.g. c("b", "se", "ci", "p")). The "p" token is always the B / beta p-value; for the AME-specific p-value use "ame_p". "none" shows the fit statistics alone, one column per model (see Fit statistics alone in the details).

keep

Character vector of regexes. Only coefficient rows whose term name (as in stats::coef() – e.g. "wt", "cyl6", "factor(cyl)8"; the terms of an rms fit in the spelling of lm(), "cyl6" for rms's "cyl=6" and "wt:cyl6" for "wt * cyl=6") matches at least one pattern are kept. The intercept and the subordinate blocks (thresholds, random effects) are exempt from keep and drop: show_intercept, show_thresholds and show_re govern them. Mutually exclusive with drop. Filtering is a display choice; p_adjust runs against the full coefficient family before filtering. Default NULL (no filter).

drop

Character vector of regexes. Coefficient rows matching any pattern are removed, with the exemptions of keep. Mutually exclusive with keep. Default NULL.

show_intercept

Whether to display the intercept row. Default TRUE (APA convention). Hide via FALSE.

show_thresholds

For ordinal cumulative-link models (MASS::polr, ordinal::clm, ordinal::clmm, survey::svyolr, glmmTMB with the ordinal family, rms::lrm with more than two outcome levels), whether to display the estimated category thresholds (cut-points) as a subordinate "Thresholds" block of rows below the predictors, carrying B / SE / CI / p like the predictor rows. Default TRUE. FALSE collapses them to a compact one-line footer note instead. Thresholds are reported on the log-odds (B) scale and are never exponentiated (under exponentiate = TRUE their rows stay on the log-odds scale). They follow the polr convention, link(P(Y <= j)) = threshold - linear predictor: an rms::lrm fit, which models P(Y >= j), shows its intercepts with the sign reversed, named "3|4" (displayed "3 | 4") rather than "y>=4". Has no effect on non-ordinal models, and the rows are shown only when a coefficient column ("b"/"beta") is in show_columns.

show_components

For models with secondary components, whether to display them as labeled subordinate blocks of rows below the primary (count / conditional / location) coefficients. Default TRUE:

pscl::zeroinfl and glmmTMB(ziformula = ): a Zero-inflation block – the model for the probability of a structural (excess) zero.
pscl::hurdle: a Zero hurdle, Pr(y > 0) block – the model for the probability of a nonzero count (note the opposite direction vs zero-inflation; the footer names each block's meaning).
A glmmTMB hurdle model (a zero-truncated family, glmmTMB::ziGamma(), glmmTMB::lognormal() or glmmTMB::beta_family(), with ziformula): a Zero hurdle, Pr(y = 0) block – glmmTMB models the probability of a zero, the opposite event of the pscl::hurdle block.
glmmTMB(dispformula = ): a Dispersion block (only when dispersion was actually modeled; log scale, never exponentiated).
ordinal::clm(scale = ~): a Scale effects block – the covariate effects on the log standard deviation of the latent response. Never exponentiated: their exponential is a ratio of latent standard deviations, not an odds ratio.

Component rows carry full Wald inference (B / SE / z / p / CI), join the p_adjust family, and take significance stars. Under exponentiate = TRUE a component is exponentiated only when its link makes the result an odds ratio (the logit zero components); probit / cauchit / cloglog zero links and count-type hurdle zero parts stay on the link scale, disclosed in the footer. FALSE omits the blocks (the title still names the model type).

intercept_position

Where to place the intercept when shown. "first" (default, APA) or "last" (Stata-style, intercept just above the fit-stats footer). Ignored when show_intercept = FALSE (with spicy_ignored_arg warning).

factor_layout

Layout of factor predictors. Applies to any categorical predictor – factor, ordered, character, or logical (R coerces the latter two to factors at fit time). Two options:

"grouped" (default): the variable name on its own header row ending with : (e.g., education:); each level follows as an indented sub-row with the bare level name. APA convention.
"flat": each non-reference dummy is one row with the <variable><level> form (e.g., educationUpper); no header, no indent. Econometrics convention.

reference_style

Rendering of factor reference levels. Four modes, distinguishing WHERE the reference information is exposed (in a row, inline, in the footer, or nowhere):

"row" (default): explicit row Female (ref.) with en-dashes in all stat columns (NEJM / BMJ clinical convention). reference_label controls the suffix.
"annotation": the row is dropped and the reference is shown inline. Under factor_layout = "grouped" the factor header reads education: [ref: Lower]; under factor_layout = "flat" the marker [vs Lower] is attached to the first non-reference dummy of each factor (subsequent dummies inherit the same reference).
"footer": the row is dropped and a sentence Reference categories: education = Lower; sex = Female. is added to the footer note. SAS PROC LOGISTIC / SPSS "Categorical Variables Codings" convention. Best for publication-grade dense multi-factor tables.
"none": the row is dropped and no reference information is displayed anywhere. The user is responsible for stating the reference convention elsewhere (article text, table caption). Under factor_layout = "flat", an informational message is emitted to flag the silent omission.

Ordered factors with AME: under R's default contr.poly, ordered factors have B coefficients named .L / .Q / .C (orthogonal polynomial trends) which have no per-level reference semantics. When "ame" is in show_columns, however, the AME block is per-level contrasts against levels()[1]. A synthetic reference row anchored on levels()[1] is therefore emitted so the reader sees the AME baseline explicitly, with the same reference_style handling as plain treatment-coded factors. The [vs <ref>] annotation in "annotation" mode is attached to the first AME row, not to the polynomial-trend rows.

reference_label

Suffix shown after the reference level in reference_style = "row" mode. NULL (default) prints "(ref.)", translated under options(spicy.language = ); a string you supply is printed as typed. Ignored by the other three modes, whose fixed wording follows the language: [ref: ...], [vs ...], and the footer line Reference categories: .... The suffix follows the level after a space, unless it opens on ,, ; or :, which attach to the level directly: a level such as T0 (2018) reads T0 (2018) (ref.); reference_label = ", ref." gives T0 (2018), ref.; reference_style = "annotation" moves the reference out of the row.

show_fit_stats

Character vector of tokens for the model-level rows below the coefficients; row order follows token order. NULL (default) resolves class-aware:

lm, estimatr::lm_robust(): c("nobs", "r2", "adj_r2"). estimatr::iv_robust() gets "nobs" alone, and an IV fixest::feols() neither R-squared nor AIC – the 2SLS R-squared is not the classical one. Asked for, it is the R-squared of ivreg(), on the structural residuals.
glm, ordinal polr / clm: c("nobs", "pseudo_r2_mcfadden", "pseudo_r2_nagelkerke", "aic") (McFadden = Stata ologit default, Nagelkerke = SPSS PLUM).
mixed lm + glm: the union of the two (the renderer en-dashes per cell the stat not defined for a given model class).

Under nested = TRUE the default is extended with the class-appropriate change-stat tokens (e.g. "r2_change", "f_change" for lm). See Vocabulary tokens (show_fit_stats subsection) and Hierarchical (nested) model comparison in the details for the full vocabulary.

fit_stats_layout

Layout of the fit-stat values (n, R^2, AIC, ...) within each model's column group. Two options:

"first_col" (default): the value is placed in the FIRST numeric sub-column of each model (typically B); the model's remaining sub-columns (SE, LL, UL, p, ...) are left empty for that row. The APA Manual 7 Table 7.13 layout.
"merged": the model's numeric sub-columns are merged into a single wide cell containing the fit-stat value, centered under the model spanner. Stata esttab layout / Econometrica and AER journal convention. Resolves the mixed-precision look of "first_col" (an integer n row sharing the B column with two-decimal coefficients).

Cell merging is supported by excel, flextable, and word (via flextable). gt, tinytable, clipboard, and default (console) always render in "first_col" mode regardless of this setting:

gt lacks a native row-spanning cell-merge API (tab_spanner covers columns, not row-cell ranges).
tinytable's style_tt(colspan = N) emits HTML colspan only on header rows, not on body cells.
clipboard ships TSV plaintext.
default ships fixed-width ASCII.

Decimal alignment of every numeric column is preserved in both modes: the B column decimal-aligns its coefficient values plus any fit-stat value(s) in "first_col" mode (native primitives handle the mixed-precision case), and trivially decimal-aligns in "merged" mode (the fit-stat values move out of the B column into the merged cell).

layout

How several models are laid out. "columns" (default) sets them side by side, one group of columns per model. "stacked" sets them one under the other, one panel per model: each panel is the table of that model alone (its default columns, the confidence interval included, its blocks and fit statistics), opened by a row with its label. See Multi-model semantics and table_stack(), which stacks tables built separately.

show_re

Logical. TRUE (default) renders the random-effects variance components of a mixed-effects fit (lmer, glmer, glmmTMB, lme, clmm) as a subordinate "Random effects" block of table rows below the fixed effects: one row per standard deviation / correlation per grouping factor, plus the residual, each with its estimate, SE, and CI in the shared coefficient columns. The group sizes (N (groups)) and the ICC render as fit-statistic rows; the footer reports the estimation method (REML / ML) and the likelihood-ratio test of the whole random part against the no-random-effects model, with the boundary-corrected chi-bar-squared p-value (Self & Liang 1987; Stram & Lee 1994). Variance-component rows deliberately carry no per-row p-value: a Wald test of a variance is invalid at the boundary of the parameter space, and no reporting guideline requests one (see the Mixed-effects models section of the Publication-ready regression tables article). FALSE suppresses the block. No effect on fits without random effects (lm, glm, coxph, ...).

re_scale

One of "sd" (default) or "variance". Controls the display scale of the random-effects rows:

"sd": report the random-effect standard deviation $\sigma$ (Gelman 2005, Technometrics: "directly interpretable as the size of the variation across groups"). Standard error and CI converted via the Delta method: $SE(\sigma) = SE(\sigma^2) / (2\sigma)$; $CI(\sigma) = \sqrt{CI(\sigma^2)}$.
"variance": report $\sigma^2$ (the canonical internal scale; SE and CI come straight from the Hessian / nlme::intervals() / glmmTMB::confint() without rescaling).

Correlation rows ($\rho$) are unitless and pass through either way.

re_columns

Character vector. Subset of c("est", "se", "ci") controlling which cells of the random-effects rows are displayed ("est" is mandatory); deselected SE / CI cells render as an en-dash on those rows only. Useful for slimming output (re_columns = "est") or for journals that want only standard errors (re_columns = c("est", "se")). Display-only: the underlying data (broom::tidy(), as_structured()) always carries the full SE + CI.

Note. Under the default re_ci = "wald" the SE and CI are Wald (est ± z * SE, clamped at 0 for variances). Wald can be optimistic near the variance boundary (Self & Liang 1987 chi-bar-squared); request boundary-respecting profile-likelihood intervals with re_ci = "profile" when robustness is critical. See the Mixed-effects models section of the Publication-ready regression tables article.

For lmer / glmer fits these SEs come from merDeriv, whose cost grows superlinearly with the number of observations (about a minute at n of roughly 2,700). Above options("spicy.re_se_max_n") (default 1000) they are skipped: the rows keep their estimates, the SE / CI cells render as en-dashes, a table note states the omission, and a spicy_caveat warning points here. Raise the cap (e.g. options(spicy.re_se_max_n = Inf)) to force the computation, or test the random terms with re_test = "lrt".

re_test

One of "none" (default), "lrt", or "rlrt". Opt-in per-term significance test for the random-effect variance components, filling the otherwise-empty p column of the Random effects rows. Never a Wald test (invalid at the boundary sigma = 0):

"lrt": likelihood-ratio test of each random term vs the model refitted without it (the term's variance plus its covariances with the other terms of its bar), referred to the boundary-corrected chi-bar-squared mixture 0.5 chi2(q-1) + 0.5 chi2(q) (Self & Liang 1987; Stram & Lee 1994). The reduction scheme matches lmerTest::ranova(); the mixture reference makes the p-value exact-asymptotic rather than conservative. A bar's intercept is tested only when it is the bar's single term. Supported: lmer, glmer, glmmTMB, and lme with a simple random = ~ terms | group structure.
"rlrt": exact restricted likelihood-ratio test with a simulated finite-sample null (RLRsim::exactRLRT(); Crainiceanu & Ruppert 2004). Only defined for a Gaussian lmer / lme fit with a single variance component.

The test statistic and df stay out of the displayed t/z column (they are chi-square-scale, not t/z) but are carried in broom::tidy() (test_type "chibar2" / "rlrt"). The whole-block LR test in the footer is unaffected. Correlation and residual rows are never tested (a correlation is tested jointly with its slope; the residual has no zero-variance null). Refits happen once per random term: expect a noticeable cost on large models. For structures outside these two routes (e.g. multiple variance components needing a finite-sample null), pbkrtest::PBmodcomp() offers a parametric-bootstrap LR test on the same nested pair of fits (Halekoh & Hojsgaard 2014) – run it directly and report its p-value alongside the table.

re_ci

One of "wald" (default) or "profile". Uncertainty route for the random-effect variance-component rows of lmer / glmer fits:

"wald": SE and symmetric CI from the observed information (merDeriv), subject to the options("spicy.re_se_max_n") size cap (see re_columns).
"profile": profile-likelihood CIs via confint(fit, method = "profile") – the route lme4 itself documents and defaults to. The intervals are asymmetric and respect the boundary at 0; no SE is shown (lme4's position: a symmetric SE misdescribes the skewed sampling distribution of a variance). Sidesteps the size cap entirely (roughly two seconds per variance parameter even at n ~ 7,000; glmer profiles cost more). The footer discloses the method, and likelihood intervals transform exactly, so re_scale = "variance" shows the profile CI of the variance itself.

glmmTMB and nlme::lme fits keep their engine-native CIs (TMB's sdreport; nlme's apVar) and refuse "profile".

model_labels

Per-model labels used as the column-group spanner above each model's sub-columns (console + gt / flextable / tinytable / Excel / Word renderers). NULL (default) resolves automatically; see Multi-model semantics for the full rule. A character vector of length length(models) overrides. Refused (error) for a single multinomial model: there the column groups are outcome categories, relabeled via outcome_labels.

outcome_labels

Optional Outcome body row override. NULL (default) hides the row entirely – under the multi-model spanner the DV is already visible above the data. A character vector of length length(models) forces an explicit Outcome row with those values (the spanner stays as "Model 1, ..." unless model_labels is also supplied). FALSE also suppresses the row. Single multinomial model (outcome-as-columns layout): repurposed as the category spanner override – a character vector with one label per non-reference outcome category (unique, in model order), e.g. c("Student vs Employed", ...); the reference category's AME-only group (when displayed) keeps its own name, and FALSE is a no-op (there is no Outcome row to suppress).

stars

Significance asterisks. FALSE (default, APA 7 Section 6.46) – no stars. TRUE – APA cutoffs c("*" = 0.05, "**" = 0.01, "***" = 0.001). A named numeric vector specifies custom thresholds, e.g. c("+" = 0.10, "*" = 0.05, "**" = 0.01, "***" = 0.001). With output = "excel" a marked estimate is written as text ("64.63***" is not a number); the cells with no marker stay numeric.

nested

Whether to inject pairwise change-statistic rows for adjacent models (M2 vs M1, M3 vs M2, ...). FALSE (default) – pure side-by-side display. TRUE – requires identical nobs and identical response variable across all models. Under a non-classical vcov, each pair is tested by a Wald test of the coefficients the current model adds, on that model's matrix, shown as "Wald F-change" or "Wald $\chi^2$ (change)" by the regime of its coefficient tests. See Hierarchical (nested) model comparison for the rule, the classes it covers and its small-sample caution.

digits

Decimal places for general numeric tokens (b, beta, se, ci, t, f_change, lrt_change, ame, ame_se, weighted_nobs; the Wald change statistics under a robust vcov follow f_change). Default 2L.

p_digits

Decimal places for p-values (p, ame_p, p_change). APA-strict: leading zero stripped, <.001 (or <.0001 etc. depending on p_digits) for small values. Default 3L; at least 1L, as in the descriptive tables. Every *_digits argument is a whole number: a fraction is an error.

effect_size_digits

Decimals for per-coefficient effect sizes (partial_f2, partial_eta2, partial_omega2) and the term-level tests (partial_chi2, partial_f). Default 2L.

fit_digits

Decimals for variance-explained / model-level effect-size fit stats (r2, adj_r2, r2_change, adj_r2_change, omega2, f2, f2_change, sigma, rmse). Default 2L.

ic_digits

Decimals for information criteria and deviance (AIC, AICc, BIC, deviance, and their _change form). Default 1L.

decimal_mark

The decimal mark, a single character: "." (the default) or "," in most styles; the Lancet style of spicy_style() uses the middle dot. When "," is used, the CI bracket separator switches to "; " automatically to avoid "0,18 [0,07, 0,30]" ambiguity. With output = "excel" a numeric cell is displayed by Excel with the viewer's locale separator, which a file cannot set, so a non-default mark makes the body go out as pre-formatted text (the default "." keeps every cell a real number).

align

Numeric column alignment. "decimal" (default) – pre-pad cells so decimal marks line up vertically (publication-style). For CI cells ([LL, UL]) the left bracket, the LL decimal point, the comma separator, the UL decimal point, and the right bracket are independently aligned across rows. "center" and "right" apply the same alignment to every numeric column.

padding

Non-negative integer giving the extra characters added to each data column's auto-computed width when the default print method renders the table. Default 0L (compact – fits more models in the same console / page width). Use 2L (Stata-like) or 4L for a more spacious look. Headers stay centered above the data region regardless of padding.

labels

How the variables are named in the table: the outcome in the title (and above its columns when the models have different outcomes), a numeric term in its row, and a factor in its header. Four forms:

labels = FALSE shows the names the model gives them (bmi, age, sex).
labels = TRUE shows each variable's label attribute (as set by haven, labelled, or label_from_names()) when it is present, not NA, and not empty, and the name otherwise: "Linear regression: Body mass index", "Age (years)", "Sex:". The attribute is read off the model frame, or off the data the model was fit on where the model frame lost it. Factor levels, interactions, and transformed terms (log(income), poly(age, 2)) are not variables and keep their names, as do the coefficient rows of factor_layout = "flat".
labels = c(age = "Age (years)", sexMale = "Men"), a named character vector, gives the variables and rows it names these labels. Every other variable follows options(spicy.var_labels): its name when the option is unset or FALSE, its label attribute when it is TRUE.
labels = var_labels(age = "Age in years") gives the variables and rows it names these labels, and every other variable its label attribute, whatever the option says. var_labels() alone is labels = TRUE.

The names are model identifiers, not display text: the outcome, the term labels of stats::terms() ("age", "sex", "age:sex", a factor's term naming its header, an interaction's term its row when it has one row), and the coefficient names ("(Intercept)", "sexMale"; those of an rms fit in the spelling of lm(), "sexMale" for rms's "sex=Male" and "age:sexMale" for "age * sex=Male"), each naming one row. They keep that form whatever the language of the table. Classes with no terms component (nls(), sampleSelection::selection()) are keyed on their coefficient (or parameter) names alone. With several models, a name is looked up in the terms of all of them and lands on the rows where it exists. A name that matches nothing in the table warns (class spicy_unknown_name), with the list of its terms, and is left out, so one vector can still serve every table of a project. The headers of the subordinate blocks (Random effects, Thresholds, Zero-inflation, ...) are set by the package and cannot be relabeled here.

The default reads options(spicy.var_labels), and a TRUE or FALSE typed in the call wins over it. Any other value, NULL included, is an error. A variable label documents the data file; a published table carries the wording you set here. note_labels gives the outcome a shorter label in the title.

note_labels

The label a variable takes where the table names it in a sentence, in the title and in the notes, in place of the one labels gives it: a named character vector, note_labels = c(sex = "sex at birth"), or a var_labels() object, which reads as one. labels can then name a row with a whole question while the note says "Missing values removed: sex at birth (2).". The rows, the column headers and the spanners keep the labels of labels, and so do the title and the notes for a variable note_labels does not name. A name that matches no column warns and is left out, as in labels, and NULL (the default) renames nothing. A note label that another variable of the table shows is an error (class spicy_label_collision): choose another one, or relabel the other variable.

title, note

Override or suppress the auto-built caption / methodological footer. Three modes per argument, those note has in every table function:

NULL (default): the package builds the standard caption ("Linear regression on <DV>" / "Hierarchical linear regression on <DV>" / ..., the outcome named by its label in note_labels, else in labels) and a methodological note, one paragraph in the order of every table of the package: what the model is (its family, design, scale of the coefficients, reference categories, ...), the cases a univariable screen fits on, the inference (VCV type, p-adjust method, ...), the glosses of the blocks, then the abbreviations the column headers and the rows of the fit statistics and of the random effects show ("AIC = Akaike information criterion"), defined as the style in force asks (spicy_style(abbreviations = ); under "jama" they open the note on a line of their own), and the sentences of note_add. Under stars, the legend of the stars follows on a line of its own.
FALSE: the corresponding banner row is omitted from every output engine, for note the legend of the stars included. Use when the surrounding manuscript provides its own caption / note.
character string (length 1): replaces the auto-built text. A title is printed verbatim. A note replaces the general note and is printed as typed (multi-line notes accepted via embedded "\n") after the lead-in of the style, which a string opening with a lead-in of its own does not repeat, the sentences of note_add after it; the legend of the stars stays on its line below, and note = "" keeps that legend alone.

Validation messages, the spanner row, and the in-body change- stat rows are not affected – they belong to the table structure, not to the banner.

note_add

Sentences to add to the note: a character vector of one or more strings, each printed as typed (no period is added) and joined to the note by a space. They close its paragraph, after the definitions of the abbreviations and before the legend of the stars, in every output and in attr(x, "note"), and the sentences the table writes stay. With a string in note, they follow it, so a wrapper can add a house sentence to every table while one table replaces its note. An NA or an empty string is an error, and so is note_add with note = FALSE, which leaves no note to add to.

output

Output type. "default" (a printable spicy_regression_table); "data.frame" (the displayed table with its numbers); "long" (one row per coefficient, broom's names); "tinytable" / "gt" / "flextable" (rich-format tables); "excel" (writes to excel_path); "clipboard" (copies to system clipboard); "word" (writes flextable to word_path).

excel_path

File path for output = "excel", ending in ".xlsx", in a directory that exists. Default NULL (required when output = "excel").

excel_sheet

Sheet name when writing to Excel. NULL (the default) uses "Regression".

clipboard_delim

Field delimiter for output = "clipboard". Default "\t" (tab-separated, pastes cleanly into Excel / Google Sheets / Word). The clipboard payload mirrors the Excel layout (title row, spanner row, header, body, footer note) but is plain text – horizontal rules, cell merging, decimal alignment and the monospace font cannot be encoded in delimited text and are therefore absent from the paste. A cell holding the delimiter itself (a label with a comma under clipboard_delim = ",", or any number under decimal_mark = ","), a double quote or a line break is quoted RFC 4180-style, so the grid survives whatever delimiter you choose.

Paste behavior by target:

Excel / Google Sheets: numerics are auto-detected and right-aligned; text cells stay left-aligned. (P-values such as .005 get re-parsed as 0.005 by Excel's auto-format – to preserve the APA leading-zero-dropped display, prefer output = "excel".)
Word: the paste is converted to a Word table; all cells start left-aligned. Apply a Table Style (Insert > Table > Design) for APA-style borders, and set right-alignment on numeric columns (Layout > Align Right). For a self-contained Word file with borders and alignment pre-applied, use output = "word" instead.

word_path

File path for output = "word", ending in ".docx", in a directory that exists. Default NULL (required when output = "word"). The Word table inherits the flextable styling (Calibri font, APA borders, decimal-aligned numerics) and adds Word-specific features: an auto-numbered caption ("Table 1: ...", "Table 2: ...") via Word's SEQ field so multiple table_regression() calls in one document number consecutively; a re-printed header row on each page break; row split prevention so a single coefficient row never wraps across two pages; and an APA-styled note line (*Note.* italic prefix per APA Manual 7 §7.14).

R Markdown / Quarto: for embedded use, prefer output = "flextable" (returns the flextable object that knits to docx/HTML/PDF natively). output = "word" writes a standalone .docx file, suited to scripted exports rather than chunk-level rendering.

word_template

Optional path to a custom .docx file used as the template for output = "word". The template's header, footer, page size, margins, and named styles ("Table Caption" in particular) are honored; the table is appended to the template body. Useful for institutional templates with pre-set headers ("APA Style", "Manuscript Submission Template", etc.). Default NULL (uses flextable's stock template).

Customizing the caption appearance: the table caption is tagged with the Word named style "Table Caption". The visual rendering (italic / bold / color / font) follows whatever that style is set to in the docx template. The stock Word template renders "Table Caption" in italic — the APA Manual 7 §7.10 condensed convention. For a different appearance (Nature-style bold non-italic, APA-strict 2-line bold-number / italic-title, etc.), edit the "Table Caption" style in a docx template and pass it via word_template = "your_template.docx". Style-based delegation keeps the rendered caption consistent with the surrounding document and lets editorial conventions (Nature, APA-strict, journal-specific) be applied without modifying the call site.

style

A journal style: a theme name ("jama", "nejm", "lancet", "annals", "apa", "aer"), a spicy_style() object, or NULL (the default). A style only changes DEFAULTS – any argument you pass explicitly wins over it. Set options(spicy.style = ) for document-wide scope. A theme covers numeric formatting, the word that opens the note and the abbreviations the note defines, not full editorial conformity; ?spicy_style lists the exact rules each one encodes and the official document they come from. An unknown name is an error.

Value

A spicy_regression_table object (a data.frame subclass with classes c("spicy_regression_table", "spicy_table", "data.frame")) when output = "default". The result carries rendering attributes (title, note, align, padding) and provenance attributes (outcome, model_ids) consumed by the print method and the broom methods. For other output values, returns the format-specific object: gt_tbl, flextable, tinytable; for "data.frame" the displayed table as a data.frame (see as.data.frame.spicy_regression_table()); for "long" a tbl_df (model_id, outcome, term, estimate_type, estimate, std.error, conf.low, conf.high, statistic, df, df.residual (the second df of an F test, NA elsewhere), p.value, test_type, nobs, r.squared, adj.r.squared, ...); the file path, invisibly, for "excel" and "word"; the table, invisibly, for "clipboard". The two data frames carry the note as printed in attr(x, "note"). With layout = "stacked", the stacked table of table_stack(), whose value is described there.

Vocabulary tokens

Two vector arguments – show_columns and show_fit_stats – accept named tokens that select what to display and in what order. All tokens are lowercase (snake_case for compound tokens). Group tokens ("all_b", "all_ame", ...) expand to a fixed vector of atomic tokens; see show_columns below.

`show_columns` – per-coefficient columns

Each token = one displayed column.

Coefficient family: "b", "beta" (standardized), "se", "ci", "t", "p".
Marginal effects: "ame", "ame_se", "ame_ci", "ame_p". "p" always refers to the B-coefficient p-value; for the AME-specific p-value use "ame_p".
Partial effect sizes – lm only: "partial_f2", "partial_eta2", "partial_omega2", each with a paired _ci companion ("partial_f2_ci", ...).
Term-level test – lm, glm and mixed-effects models: "partial_chi2" and its twin "partial_f", the term-level Type-II test. Each term is tested under the marginality-respecting hypothesis that excludes any higher-order interaction containing it (the car::Anova(type = 2) convention), so the value does not depend on the factor coding. For glm the term-level test is the F that anova(reduced, full, test = "F") reports when the dispersion is estimated, the likelihood-ratio chi-square otherwise: "partial_f", $F = (\Delta D / q) / \hat\phi$ on $(q, df)$ degrees of freedom, $q$ the term's and $df$ the residual degrees of freedom of the model, with $\hat\phi$ the dispersion summary() estimates, for gaussian, Gamma, inverse Gaussian, Tweedie, MASS::negative.binomial(theta) fitted by glm() and the quasi families (the value car::Anova(type = 2, test.statistic = "F", error.estimate = "dispersion") reports; car's default Pearson estimate differs by the IRLS convergence error only, and drop1(test = "F") divides by the deviance dispersion instead); "partial_chi2", the likelihood-ratio chi-square (Long & Freese 2014 Sections 3.2.2, 3.2.4), for binomial, Poisson and MASS::glm.nb. For lm, "partial_f" is the Type II F of car::Anova(type = 2), the F the variance-explained partials are built on, so a gaussian glm and its lm() twin print the same value. A token is refused when no model of the table carries it, and when a model carries the other one while that one is not requested (each test has its own column: ask for both). Both come from the refits and the model's own dispersion (the residual variance of an lm) whatever vcov. For lmer / glmer / glmmTMB / nlme::lme, "partial_chi2" is the Wald chi-square built from the same Type-II hypothesis – a deliberate departure from the Type-III default of SAS PROC MIXED and lmerTest (the two conventions coincide for models without interactions). Rendered as value (df) to disambiguate factor terms (k-1 df) from numeric terms (1 df), and the partial F as F (q, df) with both of its degrees of freedom, 3.45 (2, 28), as Stata and SAS report them: each cell carries the residual df of its own model, and the note defines the statistic only. The console and the rich outputs (tinytable, gt, flextable, Word, the clipboard) print the test in one column under the statistic's header, and Excel writes the statistic as a number with its df in the number format (0.00" (2, 28)"). output = "data.frame" and as_structured() keep the statistic and its df as numbers, in the statistic's column, then "<column>: df" and, for the F, "<column>: residual df". tidy() has them in df and df.residual, and inline() cites the statistic alone.
Variance explained per fit – lm only: "r2", "adj_r2". Like "n", these are per-fit columns: they are populated by table_regression_uv() screens, where each predictor block is its own model and the $R^2$ answers "how much of the outcome does this predictor explain on its own". A model whose $R^2$ is a single model-level number drops the column and reports it in the fit-statistics rows instead (the "r2" / "adj_r2" tokens of show_fit_stats, on by default for lm). Refused for classes with no least-squares variance partition (glm, coxph, ...), where only competing pseudo-$R^2$ measures exist: name the one you want in show_fit_stats.
Sample columns: "n" – per-row N, populated by table_regression_uv() screens (each predictor block is its own fit); models without per-row N data drop the column. "n_events" – outcome event counts as events/N, per factor level (reference row included) and model totals on continuous rows, computed on each model's own estimation sample (STROBE item 16's "data behind the association"; NEJM-style no. of events/total no.). Binary (binomial) outcomes and right-censored coxph fits – other models raise spicy_invalid_input.
Survival estimands – coxph and survreg: "rmst", "rmst_se", "rmst_ci", "rmst_p" (restricted-mean- survival-time difference over [0, tau]) and "risk_diff", "risk_diff_se", "risk_diff_ci", "risk_diff_p" (cumulative-incidence difference at at_time), by g-computation from the fitted model with bootstrap inference. The horizon arguments are required; see tau / at_time.

Group tokens (presets) expand to a fixed atomic vector before validation:

"all_b" -> c("b", "se", "ci", "p")
"all_b_compact" -> c("b", "se", "p")
"all_b_full" -> c("b", "se", "ci", "t", "p")
"all_beta" -> c("b", "beta", "se", "ci", "p")
"all_ame" -> c("ame", "ame_se", "ame_ci", "ame_p")
"all_ame_compact" -> c("ame", "ame_p")
"all_f2" / "all_eta2" / "all_omega2" -> partial_* + its _ci companion.

Mix groups and atomic tokens: show_columns = c("all_b", "ame", "ame_p"). Duplicates after expansion are deduplicated; the order of tokens controls the order of the displayed columns. If standardized != "none" and "beta" is not already requested, it is auto-injected after "b". Asking for "beta" while standardized = "none" raises spicy_invalid_input.

Default (show_columns = NULL) is context-aware: "all_b" for a single model (APA-7 Section 6.46 publication layout), "all_b_compact" for two or more models (CI dropped to fit the side-by-side layout; restore it explicitly when needed).

Fit statistics alone (show_columns = "none"): no coefficient row or column. Each model has one column, headed by its label (the outcome, "Model 1", a name of the list or model_labels, as for the column groups of the full table; "Model 1" for a single model), and the rows are the fit statistics of show_fit_stats, in its order, the change rows of nested = TRUE included. The title is the one of the full table; the note keeps what describes the models and their fit, and leaves out what describes coefficients. Models with different outcomes get an Outcome row when the column headers do not name them (outcome_labels = FALSE removes it). "none" goes alone, and show_fit_stats = FALSE is refused with it, as is layout = "stacked"; an empty show_columns (character(0), FALSE) is an error that points to "none". The arguments that act on coefficients only are ignored, with one spicy_ignored_arg warning naming those set to other than their default: ci_level, ci_method, tau, at_time, standardized, exponentiate, ame_type, p_adjust, keep, drop, show_intercept, show_thresholds, show_components, intercept_position, factor_layout, reference_style, reference_label, fit_stats_layout, show_re, re_scale, re_columns, re_test, re_ci, stars, effect_size_digits, and vcov, cluster and boot_n unless nested = TRUE (where a robust vcov gives the Wald change rows). An information criterion (AIC, BIC, deviance, ...) shown for models fitted to different numbers of observations or to different outcomes raises a spicy_caveat warning: its differences compare nothing. tidy() returns no row, glance() one row per model, and inline() cites a statistic of a model: inline(tbl, "aic", model = 2). table_regression_uv() refuses the token.

`show_fit_stats` – model-level rows below the coefficients

Counts: "nobs", "weighted_nobs", "n_events" (Cox models: number of events, reported alongside n per the field convention; blank for other classes). "nobs" and the BIC follow the engine's stats::nobs(), which for a weighted polr or clm fit is the sum of the weights.
Variance explained (lm only): "r2", "adj_r2", "omega2".
Pseudo-$R^2$ (glm and ordinal polr / clm): "pseudo_r2_mcfadden" (McFadden 1974), "pseudo_r2_nagelkerke" (Nagelkerke 1991), "pseudo_r2_tjur" (Tjur 2009; binomial only).
Residual scale: "sigma", "rmse". "sigma" is the lm $\hat{\sigma}$. For glm and the GLM-like fixest fits, it is the square root of the dispersion (so a gaussian glm reads as lm), and 1 when the family fixes it (binomial, Poisson). For survreg, it is the scale $\sigma$ of $\log T = X'\beta + \sigma W$ for the distributions fitted on $\log T$ (the standard deviation of $\log T$ for the log-normal), and of $T = X'\beta + \sigma W$ for those that model $T$ itself (gaussian, logistic, extreme value, t), blank when strata give one scale per stratum. For glmmTMB, it is the same square root for the gaussian, Gamma and Tweedie families, and the fixed 1 of Poisson and binomial. It is blank for the negative-binomial and beta fits below, whose dispersion has its own row, and for the other glmmTMB families, whose sigma() is a family-specific parameter. For mgcv::gam and bam, a glm family follows the glm rule; an extended family shows the square root of its scale only when the GAM estimates one (tw()), and is blank otherwise (an estimated nb(), betar(), scat(), ziP(), ocat()): mgcv fixes their scale at 1.
Survey designs (survey::svyglm only): "eff_p", the effective number of parameters of the design (the first element of survey's extractAIC, the sum of the Rao-Scott eigenvalues). Opt-in and distinct from "aic", which reports the design-based AIC of Lumley & Scott (2015). Blank for other classes.
Bayesian fits only: "r2_bayes" (posterior-median Bayesian $R^2$, Gelman et al. 2019 – in the all-Bayesian default), "elpd_loo" and "looic" (PSIS-LOO expected log predictive density and its deviance-scale twin, Vehtari et al. 2017; opt-in, a few seconds per model; the footer discloses the elpd standard error), and "waic" (Watanabe-Akaike; PSIS-LOO is generally preferred).
Negative-binomial dispersion. NB2 fits (MASS::glm.nb, lme4::glmer.nb, fixest::fenegbin, mgcv::gam with nb(), pscl::zeroinfl and pscl::hurdle with dist = "negbin", whose theta is that of the count part, and glmmTMB with nbinom2 or truncated_nbinom2): "theta" ($V = \mu + \mu^2/\theta$) and "alpha" ($= 1/\theta$, the Stata nbreg convention). NB1 fits (glmmTMB with nbinom1 or truncated_nbinom1): "delta" ($V = \mu(1 + \delta)$, the $\delta$ of Stata's nbreg, dispersion(constant)).
Beta precision: "phi". In beta regression (betareg::betareg, mgcv::gam with betar(), and glmmTMB with beta_family()), $Var(y) = \mu(1-\mu)/(1+\phi)$ (Ferrari & Cribari-Neto 2004; higher $\phi$ = less dispersion). For glmmTMB with betabinomial() or ordbeta(), $\phi$ is the precision parameter of the beta part. A constant precision is back-transformed from its link, so y ~ x | 1 reports the same $\phi$ as y ~ x.
For these dispersion rows ("theta", "alpha", "delta", "phi"), a column whose model lacks the parameter is blank, so a Poisson model can stand beside an NB2 model, or NB1 beside NB2. A parameter fixed in the call is not an estimate and has no row: MASS::negative.binomial(theta) in glm(), mgcv's negbin(theta), and nb(theta) or betar(theta) with theta > 0. Such a negative-binomial model keeps the glm rule for "sigma". A negative-binomial glmer() always shows its theta, which lme4 counts in the degrees of freedom of the model: a glmer.nb() fit, its update() refits, and a glmer() with negative.binomial(theta). For mgcv::gam and bam, theta and phi are estimated with the smoothing parameters (REML by default), so they differ from the maximum-likelihood values of glm.nb(), betareg or glmmTMB. method = "ML" gives those. A token is refused only when no model in the table has the parameter, or when a model of that family has dispersion covariates (a glmmTMB dispformula, or a betareg y ~ x | z), so the parameter is not a single number.
GEE (geepack::geeglm only): "qic" and "qicu" – the quasi-likelihood information criteria (Pan 2001; QIC for working-correlation choice, QICu for covariate choice) via geepack::QIC(), computed only when requested (QIC silently refits the independence working model); "scale" – the estimated scale (dispersion) parameter, shown only when the fit estimated it (blank under scale.fix = TRUE, matching geepack's own summary); "max_cluster_size" – the largest cluster in the estimation sample (on by default with "n_groups", which renders the N (<id>) cluster count). Cluster statistics count the clusters geepack itself used (geese$clusz, consecutive runs of id = values), so they always describe the displayed sandwich inference. Refused when no model in the table is a geeglm fit; the likelihood-based tokens ("aic", pseudo-$R^2$, ...) are refused in return for all-GEE tables – quasi-likelihood has no likelihood.
Absorbed fixed effects (fixest fits, and estimatr fits built with fixed_effects =): "fixed_effects" renders a Fixed effects: block at the top of the fit statistics – one Yes / No row per absorbed factor (the etable / esttab convention), blank cells for models without the concept in mixed tables – and is on by default for both. Varying-slope-only factors (Origin[[x]]) absorb no intercept: they read No when another model absorbs that factor, and contribute no row otherwise. "within_r2" is the FE-partialed within R-squared (a default for feols, except an IV one; opt-in for estimatr; GLM-family fixest fits report fixest's McFadden pr2 instead). Both tokens are refused when no model in the table absorbs fixed effects.
Effect size: "f2".
Information criteria: "aic", "aicc", "bic", "deviance" (lowercase like every other token; the rendered row labels stay AIC / AICc / BIC). The AIC and BIC of a fixest model are those of its lm(), glm() or MASS::glm.nb() twin with the fixed effects as dummies: glm()'s log-likelihood for a Gaussian, Gamma or inverse-Gaussian feglm(), with the residual variance, dispersion or theta counted, which fixest's AIC() leaves out. A Gaussian femlm() follows the same rule; Poisson, binomial and IV fits keep fixest's values. A lean fit (lean = TRUE) that lost what they need, its weights or the fitted values of a Gamma or inverse-Gaussian feglm(), shows no AIC, BIC or nested rows.
Change-stats for hierarchical comparison (active under nested = TRUE; see Hierarchical comparison below): "r2_change", "adj_r2_change", "f_change", "f2_change", "lrt_change", "aic_change", "aicc_change", "bic_change", "deviance_change", "p_change". The differences, current model minus previous ("r2_change", "adj_r2_change", "f2_change", and the AIC, AICc and BIC changes), print with an explicit sign (+0.07, -9.4) in every output. The test statistics and "p_change" are never signed. The test statistics carry their degrees of freedom, as the partial F cell does: 12.38 (1, 29) for an F, the number of added parameters and the residual df of the larger model (the denominator df for the Wald F of a robust vcov), 9.12 (1) for a chi-square. A df that is not a whole number (the HTZ denominator under "CR2", the df of a smooth term) prints with one decimal, 174.77 (2, 0.8). The rich outputs print the same cell, inline() cites the statistic alone, output = "data.frame" keeps it as a number, and Excel writes the number with the df in its number format. "deviance_change" is the drop in deviance (previous minus current), positive when the larger model fits better, and is not signed either.

Default (resolved when NULL) is class-aware: lm and lm_robust fits get c("nobs", "r2", "adj_r2"); glm and ordinal polr / clm fits get c("nobs", "pseudo_r2_mcfadden", "pseudo_r2_nagelkerke", "aic"); mixed lm + glm sets union both groups (the renderer per-row en-dashes the inappropriate cell); Cox fits get c("nobs", "n_events", "aic"); GEE fits get c("nobs", "n_groups", "max_cluster_size") (the cluster structure, Stata xtgee-header style). When nested = TRUE, the class-aware default is extended with change tokens (c("r2_change", "f_change", "p_change") for lm, c("lrt_change", "p_change") for glm, c("f_change", "p_change") for a glm whose dispersion is estimated). The order of tokens in show_fit_stats controls the order of the rows.

Multi-model semantics

Pass a single fit or a list() of fits. Multi-model layout draws a centered spanner label above each model's sub-columns:

list("Naive" = m1, "Adjusted" = m2) -> spanner labels "Naive" / "Adjusted". Partial naming (list("Naive" = m1, m2)) auto-fills missing slots as "Model <position>".
list(m1, m2) (unnamed) -> if all response variables differ, each outcome becomes the spanner label of its model and the redundant Outcome body row is suppressed. The outcome is named as in the title: its name (from formula(fit)[[2]]), or its label attribute under labels = TRUE. If DVs match, the labels default to "Model 1, 2, ...".
model_labels = c("A", "B") overrides everything.

A default label is worded in the language of the table (options(spicy.language)), over the columns and in the notes that cite the model; the column names of as.data.frame() and as_structured() keep the English one ("Model 1: B").

Duplicate explicit names in the list are rejected (spicy_invalid_input) – they would silently collide in the internal model_id key. So is a name that collides with the auto-filled label of another slot (list("Model 2" = m1, m2)): two models under one spanner cannot be told apart in the table, nor addressed by inline(). Rename the model, or pass model_labels.

Stacked layout

layout = "stacked" sets the models one under the other, one panel per model, through table_stack(). Each panel is exactly the table table_regression() gives that model alone, with the same arguments: its default columns ("all_b", the confidence interval included), its reference rows, blocks and fit statistics. A row with the label of the model opens each panel: model_labels, else the names of models, else the outcomes when they all differ, else "Model k". The panels share one column header when they have the same columns, and each opens with its own header otherwise (an odds ratio beside a coefficient, a multinomial model with its outcome categories in columns). vcov and cluster lists give each panel its own estimator; outcome_labels does not apply. Under nested = TRUE, the change statistics of each model sit in its own panel, labeled "vs" the model before it. The type of each model, its standard errors and the adjustment of its p-values are worded as in the note of the models side by side, each panel named by its label; the other sentences each open with the labels of the panels they concern. A warning or an error names a model as the side-by-side table does. tidy(), glance() and output = "long" return the rows of the models side by side, under a leading panel column, with one exception: under keep or drop, a multinomial model returns the rows its panel shows, selected by its terms as in its table alone ("age"), where side by side the patterns meet terms named with their outcome ("Student: age"). output = "data.frame" binds the data frames of the panels by rows. A keep, drop or show_intercept = FALSE that leaves a panel without a coefficient row is an error: side by side, such a model keeps its fit statistics beside the others.

Inference and standard errors

vcov selects the variance-covariance estimator:

"classical" – OLS (lm) / Fisher information (glm).
"HC0" to "HC5" – heteroskedasticity-consistent (via sandwich::vcovHC()).
"CR0" to "CR3" – cluster-robust with Satterthwaite-corrected df (via clubSandwich::vcovCR()). Requires cluster.
"bootstrap" – nonparametric or cluster bootstrap (boot_n replicates).
"jackknife" – leave-one-out / leave-one-cluster-out.

For multi-model use, both vcov and cluster accept a single value (recycled to all models) or a list (one per model). The same fit can appear several times with different estimators to compare standard errors side-by-side.

Inferential regimes for lm / glm (B and AME share the same regime):

classical, HC* -> t with df.residual (lm). For glm, the test summary() of the model reports: t on the residual degrees of freedom when the dispersion is estimated (gaussian, Gamma, inverse Gaussian, Tweedie, negative.binomial(theta) and the quasi families), z otherwise (binomial, Poisson, MASS::glm.nb).
bootstrap, jackknife -> z asymptotic.
CR0-CR3 -> t with Satterthwaite-corrected df (B via clubSandwich::coef_test(); AME via clubSandwich::linear_contrast(); Pustejovsky & Tipton 2018). Under non-linear terms (poly(), I(), log(), splines::ns()), AME falls back to z-asymptotic with a spicy_fallback warning.

For other classes the t-vs-z axis follows the estimator's native reference distribution (e.g. glmmTMB, survival, ordinal, betareg, mlogit use z; lm, lme, lmer, rms::ols use t).

Multinomial models: outcome categories as columns

A single nnet::multinom model renders the publication layout: predictors as rows, one column group per non-reference outcome category (spanner = category name), so a predictor's effect can be compared across equations along its row. The mandatory footer note Reference outcome: <level>. names the base category (Stata's "base outcome" line); outcome_labels relabels the spanners (e.g. "Student vs Employed"). Fit statistics print once, under the first group. When per-category AMEs are requested the reference category appears as a last, AME-only group (its coefficient cells are empty – the reference has no equation; its AMEs complete the zero-sum across categories). The default show_columns compacts to B / SE / p exactly like multi-model tables; request atomic tokens to restore CIs. Multi-model and nested = TRUE multinomial tables keep the one-row-per-(category, predictor) layout – categories-within-models would need two spanner levels. tidy() and output = "long" always return the long form ("<category>: <term>" rows), whatever the display; as_structured() and output = "data.frame" mirror the displayed table (one column set per category), as they do for every other layout.

Robust SE availability by model class

Not every estimator is defined for every class. A robust vcov the class cannot honor fails fast with spicy_unsupported_vcov – never a silent model-based result under a robust label:

lm, glm, MASS::glm.nb: all of classical, HC*, CR*, bootstrap, jackknife. lm additionally takes "CR1S" – the full Stata regress, vce(cluster) convention: the CR1S small-sample scaling (equal to sandwich::vcovCL() with type = "HC1") with t(G - 1) inference, G the number of clusters. Use it to reproduce Stata tables; "CR2" (Bell-McCaffrey, Satterthwaite df) remains the recommended modern choice. "CR1S" is refused for glm with an explanation: Stata's ML commands use a different convention (G/(G-1) scaling, z), so the label would not match Stata output there. For MASS::glm.nb, the sandwich holds theta at its estimate, so robust standard errors can differ by a few percent from software that includes the score of theta.
mlogit: classical + "CR1" only (cluster at the choice-situation level) – sandwich::vcovHC() mis-scales the sandwich for mlogit's per-choice-situation scores, so HC* is refused.
lmer, lme: classical + CR* only – HC* and the resamplers (which refit lm/glm) are not defined for these.
coxph, rms (ols/lrm/cph/Glm): classical + "CR0" only: the Lin-Wei sandwich of coxph(cluster = ) and rms::robcov() are one cluster estimator each, with no small-sample factor. "CR1" to "CR3" are refused, as are HC* and the resamplers.
survreg, mgcv::gam/bam, polr, clm, betareg, nnet::multinom, pscl::zeroinfl/hurdle: classical + "CR1" only: sandwich::vcovCL() computes one cluster estimator for these, the CR0 sandwich scaled by $G/(G-1)$, G the number of clusters (clubSandwich's CR1). "CR0", "CR2" and "CR3" are refused, as are HC* and the resamplers. clm with a scale / nominal (partial-PO) component is classical only. multinom needs sandwich >= 3.1-2 (which added its estfun() method); its cluster is one entry per observation. For the two-part pscl fits the cluster sandwich covers both components (count and zero) at once. A negative-binomial zeroinfl or hurdle is classical only: the sandwich of these fits leaves out the score of theta. When the count part is negative binomial, fit the model with glmmTMB instead (family nbinom2 with ziformula, or truncated_nbinom2 for a hurdle), whose sandwich includes it; glmmTMB's zero hurdle is binomial, as with zero.dist = "binomial". When only the zero hurdle is negative binomial, the default zero.dist = "binomial" keeps "CR1".
quantreg::rq: its own estimator family, not the sandwich vocabulary: classical (= "nid", quantreg's large-sample default), "iid", "ker", "rank" (intervals only) and a native "bootstrap", clustered via the wild gradient bootstrap. HC*, CR* and jackknife are refused, each with its own reason – see the vcov argument.
geepack::geeglm: no spicy-side estimator at all: GEE inference is robust by construction, so the fit's own sandwich (or jackknife) standard errors – chosen by geeglm's std.err = option, clustered on its id = – are the displayed inference. HC* / CR* and cluster are refused with a pointer to those fit options.
survey::svyglm: classical only, on principle: the design-based Taylor / replicate variance is already the robust variance for the declared design, and clustering belongs in the design itself (survey::svydesign(ids = )), not in the table call.
estimatr (lm_robust/iv_robust) and fixest: these fits carry the variance they were computed with – estimatr's se_type =, and whatever fixest's own vcov interface produced – and spicy never overwrites it. A vcov request is refused; set the estimator through the fitting package instead (se_type = for estimatr; vcov = at estimation or in summary() for fixest). With IID errors, fixest's default small-sample correction (fixest::ssc()) multiplies the standard errors of a feglm() fit by sqrt((n - 1) / (n - K)) over those of its glm() twin with the fixed effects as dummies (a feols() fit already matches lm()), and summary(fit, ssc = fixest::ssc(K.adj = FALSE)) gives the glm() values.
glmmTMB: classical, "CR0" and "CR1", from glmmTMB's own sandwich (glmmTMB >= 1.1.15), which covers the coefficients, the zero-inflation and dispersion blocks and the AME. "CR0" is the sandwich as glmmTMB computes it; "CR1" scales it by $G/(G-1)$, G the number of clusters, as for the sandwich::vcovCL() classes. Both use Wald z. Without cluster, the cluster is the top-level grouping factor of the random effects. The estimator needs independent clusters that contain the random effects – every group within one cluster – so crossed random effects need a cluster that contains them all. It keeps the standard errors valid when the random-effects distribution or the conditional variance is misspecified, provided the estimates themselves stay consistent (Liang & Zeger 1986; Verbeke & Lesaffre 1997); with few clusters it understates the variance, which the factor of "CR1" offsets only in part: this class has no bias-reduced correction (Cameron & Miller 2015). "CR2", "CR3" and HC* are refused (the sandwich has one score per cluster, with no residuals and no hat values), as are REML fits and the ordinal family.
Other classes (glmer, rstanarm/brms, ...): classical (model-based) only (clubSandwich has no working backend for glmer).

For the mixed-effects fits (lmer, lme, glmmTMB), the robust estimator covers the fixed effects: their coefficients, their joint tests ("partial_chi2"), the standardized coefficients and the AME. The random-effect variance components keep their model-based inference. glmmTMB's sandwich takes one likelihood score per cluster over the fixed and variance parameters, while the CR* estimators of lmer and lme come from clubSandwich over the fixed effects alone: on the same model the two agree as the clusters grow, but not exactly.

Cluster-robust backends differ by class but are each cross-validated to the field-standard oracle: lm/glm/lmer/lme use clubSandwich (CR2 = Bell-McCaffrey, with Satterthwaite df for lm/lme/lmer); coxph/cph use the Lin-Wei grouped-dfbeta sandwich (identical to coxph(..., cluster=)); survreg/gam/polr/clm/betareg/mlogit/multinom and the Poisson and geometric pscl two-part fits use sandwich::vcovCL(). For gam/bam, spicy forms the sandwich itself, from the GLM score and the penalized covariance of the fit, and lets vcovCL() do only the clustering arithmetic: the sandwich methods for these fits are wrong for a free dispersion, a non-canonical link and a non-Gaussian bam(). Extended families other than nb() (betar, ocat, scat, tw, ...) have another score and are refused. For all of these vcovCL() classes, the one estimator is CR1, requested as "CR1"; the other CR* tokens are refused rather than mapped. rms fits use rms::robcov() (which needs the fit's x = TRUE, y = TRUE). The Lin-Wei and robcov() sandwiches have no small-sample factor: their one estimator is requested as "CR0", and the other CR* tokens are refused rather than mapped. glmmTMB fits use glmmTMB::sandwich(), with one score per cluster (the gradient of the cluster's marginal log-likelihood) and no finite-sample factor: "CR0" requests it as computed, "CR1" scaled by $G/(G-1)$, and "CR2" and "CR3" are refused rather than mapped. cluster length is one entry per observation, except mlogit (one per choice situation) and censored coxph (one per subject).

How to specify `cluster`

Three accepted forms, in order of preference:

Formula – cluster = ~region (or cluster = ~region:year for the interaction of two variables). The variables are looked up in model.frame(fit) first, then in the original data argument captured by the fit. Recommended: independent of the dataset's name, composable for multi-way clustering, consistent with sandwich::vcovCL() / clubSandwich::vcovCR().
String – cluster = "region". A single column name resolved the same way as the formula. Convenient but cannot express interactions.
Vector – cluster = df$region. An atomic vector of length nobs(fit). Use this when the cluster key is derived on the fly (cluster = interaction(df$region, df$year), cluster = as.integer(format(df$date, "%Y"))), comes from a different dataset with matching row order, or is otherwise not a column of the model's data.

A bare unquoted name (cluster = region) is not column selection: it is evaluated as an ordinary R variable. If an object region exists in the calling environment, its value is used as the vector form above; if it does not (the typical "unquoted column name" intent), the call fails with a migration error pointing at ~region / "region". Bare-name column selection is deliberately unsupported – it would require non-standard evaluation magic that breaks under programmatic use (function wrapping, dynamic column choice, loops).

For multi-model use, mix forms freely: cluster = list(~region, "region", df$region).

Bayesian fits (rstanarm / brms)

stanreg and brmsfit models are summarized from their posterior draws: B is the posterior median, SE the posterior MAD SD (the scaled median absolute deviation – the pairing of Regression and Other Stories and rstanarm's own print; equal to the posterior SD for a normal posterior), and the interval an equal-tailed credible interval (header 95% CrI; ci_method = "hdi" opts into the highest-density interval, header 95% HDI). There is no p-value or t-statistic; group presets ("all_b", ...) expand without them, and the opt-in "pd" column reports the probability of direction (Makowski et al. 2019), with a footer definition. Under exponentiate = TRUE all quantities are computed from the exponentiated draws directly (the SE is the posterior MAD SD on the ratio scale, not a delta-method approximation). The AME columns are equally draws-native: avg_slopes() computes the marginal effects per posterior draw, and the table reports their posterior median, MAD SD and credible interval (equal-tailed or HDI per ci_method); "ame_p" has nothing to fill and follows the p-column policy (dropped from presets, refused as an atomic token, dashed in mixed tables). Standardized betas follow the same draws logic: "posthoc" / "basic" / "smart" are exact affine rescales of the link-scale summaries (the Gaussian family divides by SD(y); every other family standardizes predictors only, the frequentist glm convention); "refit" and "pseudo" are refused.

Fit statistics: "r2_bayes" (in the Bayesian default) plus the opt-in "elpd_loo" / "looic" / "waic", whose standard errors are disclosed in the footer; unreliable estimates (PSIS-LOO Pareto k above the sample-size-specific threshold $\min(1 - 1/\log_{10} S, 0.7)$ of Vehtari et al. 2024 – the same bound loo::loo() prints – and WAIC p_waic > 0.4) add a footer caveat instead of being silenced. Every table is backed by an automatic sampler-diagnostics guard (R-hat >= 1.01, ESS below 100 per chain – floored at 400 so fewer chains never weaken the bar –, divergent transitions, E-BFMI < 0.2, per Vehtari et al. 2021): problems add a footer line and raise a warning classed spicy_bayes_diagnostics (nested under spicy_caveat, so withCallingHandlers(spicy_bayes_diagnostics = ...) mutes the guard selectively); clean fits print nothing. Per-coefficient "rhat" / "ess_bulk" / "ess_tail" columns are available in all-Bayesian tables, as is "mcse" – the Monte Carlo standard error of the displayed posterior median, the criterion for how many digits a table can honestly show (a displayed digit is Monte-Carlo stable when twice the MCSE stays below it; Gelman, Vehtari, McElreath et al. 2026, sec. 11.6). Under exponentiate = TRUE the MCSE is recomputed on the exponentiated draws.

Refused on principle (classed error, never a silent fallback): likelihood-based fit statistics ("aic", pseudo-R², ...), p_adjust, robust / cluster vcov (model the clustering with group-level terms instead), ci_method = "profile" / "boot_percentile", and non-MCMC fits (algorithm = "meanfield" / "optimizing" – refit with algorithm = "sampling"). See the Bayesian regression tables article.

Hierarchical (nested) model comparison

nested = TRUE adds per-pair change statistics as in-table 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; the FIRST model column gets en-dashes (no previous model to compare to). Validation requires identical nobs and identical response variable across all models.

Structural refusals. Four hierarchies carry no valid change statistic whatever their nobs and response, and are refused with spicy_invalid_input rather than rendered with empty or meaningless cells:

A REML pair whose fixed effects differ (nlme::lme(method = "REML"), glmmTMB(REML = TRUE)). The restricted likelihood depends on the fixed-effects design matrix, so the two criteria are not on a common scale and the likelihood-ratio test is not valid (Pinheiro & Bates 2000 Section 2.4.2). Refit both by maximum likelihood (method = "ML", REML = FALSE). A REML pair that differs only in its random structure stays valid and is accepted.
Two fixest fits that absorbed different fixed effects. fixest counts every absorbed level as an estimated parameter, so a test across two | fe structures measures the change in the absorbed effects together with the change in the coefficients, and the two cannot be separated. Hold the | fe part fixed across the hierarchy.
Mixed-effects fits from different engines (lme4, glmmTMB, nlme::lme). The change statistics come from the engine's own two-model anova(), and no engine's method accepts a fit produced by another. Refit one of the two with the other's engine.
A MASS::rlm hierarchy. M-estimation minimizes a bounded loss on the scaled residuals: it partitions no sums of squares and it is not a likelihood, so there is neither a likelihood-ratio test nor a partial F to report. Every change token is refused – including on a single fit or a nested = FALSE pair, where the token is the explicit request. Compare the least-squares fits (lm()) if the hierarchy itself is the question.

In the first three, nested = FALSE renders the models side by side instead.

Default change tokens auto-injected when show_fit_stats is NULL:

All-lm (and nls): c("r2_change", "f_change", "p_change") – APA hierarchical regression standard.
All-glm: c("lrt_change", "p_change") – Hosmer & Lemeshow Section 3.5; Long & Freese 2014 Section 3.2.4.
All-glm whose dispersion is estimated (gaussian, Gamma, inverse Gaussian, Tweedie, quasi families), and the feglm fits of these families: c("f_change", "p_change"), the F test anova(test = "F") reports, $F = (\Delta D / q) / \hat\phi$ on $(q, df)$, with $\hat\phi$ the dispersion and $df$ the residual degrees of freedom of the larger model. "lrt_change" is refused for such a hierarchy.
Any other hierarchy that carries a likelihood (survreg, polr, clm, coxph, gls, betareg, multinom, fixest, ...): the same c("lrt_change", "p_change").
Mixed-effects (lmer / glmer / glmmTMB / nlme::lme): c("aic_change", "bic_change", "lrt_change", "p_change") – Pinheiro & Bates 2000 Section 2.4.1.
Quantile (rq): c("f_change", "p_change") – the Wald-type test quantreg::anova.rq() reports.

To customize, pass the change tokens directly to show_fit_stats. The variance-explained change tokens ("r2_change", "adj_r2_change", "f_change", "f2_change") raise spicy_invalid_input on any hierarchy whose nested comparison is a likelihood-ratio test – glm ("f_change" excepted when its dispersion is estimated), and equally every other likelihood class and the mixed-effects families: that comparison reports a chi-square, not a variance-explained change, and the message points at "lrt_change" + "p_change". A quantile hierarchy refuses them too, and "lrt_change" with them, pointing at "f_change" + "p_change" instead. lm and nls keep the least-squares tokens.

The change test follows vcov. Under the classical vcov the change test is the partial F, the likelihood-ratio test, the F test anova() reports for a pair whose dispersion is estimated, or quantreg::anova.rq()'s Wald test, as above. Under any other vcov, and under the classical token when the fit carries a robust variance of its own, it is the Wald test of the coefficients the current model adds, on the matrix its coefficient rows were computed from: $W = b' V^{-1} b$, with $b$ the $q$ added coefficients (aliased ones left out) and $V$ their block of the matrix. The test takes the regime of the current model's coefficient tests:

F, in a row labeled "Wald F-change": $W / q$ on $(q, df)$ degrees of freedom, the convention of lmtest::waldtest(test = "F") and car::linearHypothesis(). It applies to lm under "HC0"-"HC5" (and "CR1S", on $G - 1$ degrees of freedom), glm with an estimated dispersion under "HC0"-"HC5", rms::ols under "CR0" or after rms::robcov(), quantreg::rq under "iid" and "ker", and a fixest fit estimated with a non-iid vcov whose coefficient tests are t (on fixest's own degrees of freedom, as fixest::wald()): feols, and feglm with a family whose dispersion is estimated (gaussian, Gamma, inverse Gaussian and the quasi families). Under "CR0"-"CR3", lm, glm (MASS::glm.nb included), lme4::lmer and nlme::lme use the HTZ test of clubSandwich::Wald_test(), the small-sample correction of their coefficient tests.
Chi-square, in a row labeled "Wald $\chi^2$ (change)" that replaces the likelihood-ratio row: $W$ on $q$ degrees of freedom, the convention of lmtest::waldtest(test = "Chisq"). It applies to glm with a fixed dispersion and glm.nb under "HC0"-"HC5", lm and glm under "bootstrap" and "jackknife", rq under "bootstrap", glmmTMB under "CR0" and "CR1", under "CR0" to coxph and rms::lrm / cph / Glm, and under "CR1" to survreg, mgcv::gam and bam, polr, clm, betareg, multinom, mlogit, zeroinfl and hurdle (not negative binomial); also to survreg(robust = TRUE), a fixest fit with a non-iid vcov whose coefficient tests are z (fepois, femlm, fenegbin, and feglm with the binomial or Poisson family), and an rms fit of these classes passed through robcov().

$\Delta R^2$, the information criteria and the deviance change do not depend on the vcov and are unchanged, and a vcov list tests each pair with the current model's estimator. The previous model's coefficients must all appear, under the same names, in the current one, and a smooth term added to a gam under a robust vcov is refused (spicy_invalid_input in both cases). rq keeps anova.rq()'s test under "nid", its default, and under "rank", which yields no matrix. A Wald test of many constraints in a small sample is liberal: on 219 observations, a block of 19 predictors gave F(19, 178) = 6.91 classically and 9.06 as a Wald HC3 test. The HTZ correction under "CR0"-"CR3" mitigates it.

Standardized coefficients

standardized controls the method when "beta" is in show_columns:

"refit" – refit on z-scored data. For lm both X and Y are z-scored (Cohen et al. 2003 gold standard); for glm only numeric X (Long & Freese 2014 Section 4.7.2 "x-standardization").
"posthoc" – post-hoc scaling. lm: $\beta = B \times SD(X) / SD(Y)$; glm: X-only $\beta = B \times SD(X)$ (Y is undefined on the link scale).
"basic" – like "posthoc" but factor dummies are scaled by their column SD.
"smart" – Gelman (2008): continuous numeric inputs are scaled by 2 * SD(X) (a +/-1 SD swing then spans the same range as a binary's 0 to 1 step); binary inputs – numeric 0/1 and factor dummies – are left unscaled. (Before 0.13.0 the rule was applied inverted; see NEWS.)
"pseudo" – glm only. Menard (2004, 2011) fully-standardized $\beta = B \times SD(X) / SD(Y^*)$, with $Y^*$ the latent variable on the link scale and $SD(Y^*) = \sqrt{Var(\hat{\eta}) + Var_{link}}$ ($\pi^2/3$ logit, 1 probit, $\pi^2/6$ cloglog). Binomial families only; non-binomial returns NA with a spicy_caveat.
"none" (default) – no $\beta$ computed.

Interactions and transforms. Under "refit", an interaction's $\beta$ is the coefficient of the product of the z-scored components – the recommended treatment for such models (Cohen et al. 2003 Section 7.7; Aiken & West 1991; Friedrich 1982). Under "posthoc", "basic", and "smart", the product / transformed design column is treated as a single numeric column and scaled by its own SD (under "smart", by 2 * SD when the product column is continuous; a binary product column – e.g. a binary-by-binary interaction – stays unscaled like any binary input) – the convention of SPSS beta, Stata regress, beta, SAS PROC REG STB, and lm.beta::lm.beta(), identical to effectsize::standardize_parameters(method = "basic") on those columns. The two conventions differ whenever the components are correlated; a spicy_caveat warns and the footer names the convention actually used. "refit" declines formulas with inline transforms (log(x), poly(), factor() written in the formula): the model frame's evaluated columns cannot be re-evaluated on z-scored data, so spicy falls back to "posthoc" with a warning (effectsize's refit instead standardizes before the transform – a different estimand). Pre-build transformed columns in data for an exact refit.

Multiple-comparison adjustment

Adjusting the p-values of all coefficients of a single regression model is not the standard convention. Each coefficient tests a distinct hypothesis on a distinct predictor – not the situation multiple-testing procedures were designed for (Rothman 1990; Greenland 2017; APA Manual 7 Section 6.46; Harrell Regression Modeling Strategies Section 5.4; Gelman, Hill & Yajima 2012). Hence the default p_adjust = "none".

Adjustment is appropriate for: mass screening with no prior hypothesis (typically "BH" / FDR), pre-registered multi-endpoint confirmatory designs (typically "holm"), or when a journal / SAP explicitly requests it.

The adjustment runs before any keep / drop filtering, so the family is the model's full coefficient set (intercept and reference rows excluded), not the displayed subset – filtering is a display choice and must not change the inferential family.

Tied event times and the survival estimands

A Cox model has to decide what to do when several events share the same time, and survival::coxph() decides with ties = "efron" by default. The absolute estimands do not make a second decision:

The baseline hazard behind the RMST-difference and risk-difference columns follows the tie-handling convention of the fit, exactly as survival::survfit() and survival::basehaz() do. A ties = "breslow" fit gives a Breslow baseline; the default Efron fit gives an Efron baseline. One rule, so the hazard-ratio column and the dRMST column always come from the same likelihood. To change the convention, change the fit.

When the choice matters is not simply "when there are many ties". What governs it is the fraction of tied events within the risk set, together with the effect size and the spread of the covariates: a large risk set dilutes a tied block, a small one – fine strata, matched or nested case-control designs – does not. A file with hundreds of tied days can be insensitive to the convention while a small stratified one with a handful of ties is not.

And Efron is the better default, not the correct answer: the approximations are all biased at finite sample size and agree asymptotically, and Efron's virtue is a much smaller small-sample bias (Hertz-Picciotto & Rockhill 1997). Efron himself, quoting Peto, put it as "it probably doesn't make much difference" (Efron 1977).

Output formats and broom integration

output selects the return type:

"default" – a spicy_regression_table (data.frame subclass) printed via spicy_print_table().
"data.frame" – the displayed table with its numbers, as.data.frame() of the default object: one row per printed row, the key columns variable, label, level and .row_role, then the printed columns, numeric, each interval split into its two bounds.
"long" – the long-format tibble, one row per (model_id, term, estimate_type), with broom's column names.
"tinytable" / "gt" / "flextable" – rich-format HTML / Word / PDF tables (require the corresponding Suggests package). In a PDF, a "flextable" table wider than the page can be fitted with flextable::fit_to_width().
"excel" – writes to excel_path via openxlsx2::write_xlsx() and returns the path invisibly.
"word" – writes to word_path via flextable::save_as_docx() and returns the path invisibly.
"clipboard" – copies to the system clipboard via clipr::write_clip().

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 with the model-level statistics; df.residual is kept numeric so cluster-robust Satterthwaite df is preserved.

Global options

options(spicy.style = ) – the journal style applied by default to all nine table functions (table_categorical(), table_categorical_svy(), table_continuous(), table_continuous_svy(), table_continuous_lm(), table_outcome(), table_multiresponse(), table_regression() and table_regression_uv()). Takes a theme name ("jama", "nejm", "lancet", "annals", "apa", "aer") or a spicy_style() object. The style of a report is a property of the document, so it is set once in the setup chunk; the style argument overrides it per call, and any formatting argument you type overrides both. options(spicy.style = NULL) restores spicy's own defaults. See spicy_style() for the exact rules each theme encodes and the official document they come from.
options(spicy.note_style = ) – how the table note is rendered by the "tinytable" engine (all nine table functions). A note is subordinate to the table it documents, so it is set one size down, in black: 0.9em (\small in LaTeX), matching the "gt" engine and the 10pt-on-11pt of the Word engine. Two opt-ins:
- "none" – no intervention; the note is rendered exactly as the receiving document template styles it.
- any other string – extra arguments for the Typst text() call around the note, appended to the size (e.g. options(spicy.note_style = "fill: luma(89)") for a gray note). Typst only; the HTML note keeps the plain 0.9em.
options(spicy.note_style = NULL) restores the default.

Weights

No weights argument: weights are a property of the fit (extracted via stats::weights()). Pass them when fitting: lm(y ~ x, data = df, weights = w). All downstream computations (vcov, AME, standardization, weighted_nobs) extract them automatically.

Internationalization

Output follows options(spicy.language) (see spicy_labels()); the column names of as.data.frame() and as_structured() stay English in every language. Override user-facing strings via reference_label, model_labels, outcome_labels, and labels. The title and footer are post-processable via attr(result, "title") and attr(result, "note").

Classed conditions

Every error and warning emitted by table_regression() carries a classed condition for programmatic dispatch via tryCatch() or withCallingHandlers(). Errors inherit from spicy_error (root); warnings from spicy_warning. Specific leaves used by this function include spicy_invalid_input, spicy_invalid_data, spicy_unsupported, spicy_unsupported_vcov, spicy_unsupported_standardized, spicy_missing_pkg, spicy_missing_column, spicy_ignored_arg, spicy_caveat, spicy_fallback. See spicy for the full taxonomy.

References

APA Manual 7 (American Psychological Association, 2020), Tables 7.13-7.15.

Aiken, L.S. & West, S.G. (1991). Multiple regression: Testing and interpreting interactions.

Cameron, A.C. & Miller, D.L. (2015). A practitioner's guide to cluster-robust inference. Journal of Human Resources, 50(2), 317-372.

Cohen, J., Cohen, P., West, S.G., & Aiken, L.S. (2003). Applied multiple regression / correlation analysis for the behavioral sciences (3rd ed.). Lawrence Erlbaum.

Davison, A.C. & Hinkley, D.V. (1997). Bootstrap methods and their application. Cambridge University Press.

Efron, B. (1977). The efficiency of Cox's likelihood function for censored data. Journal of the American Statistical Association, 72(359), 557-565.

Friedrich, R.J. (1982). In defense of multiplicative terms in multiple regression equations. American Journal of Political Science, 26(4), 797-833.

Gelman, A. & Pardoe, I. (2007). Average predictive comparisons for models with nonlinearity, interactions, and variance components. Sociological Methodology, 37(1), 23-51. doi:10.1111/j.1467-9531.2007.00181.x

Hertz-Picciotto, I. & Rockhill, B. (1997). Validity and efficiency of approximation methods for tied survival times in Cox regression. Biometrics, 53(3), 1151-1156.

Liang, K.-Y. & Zeger, S.L. (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73(1), 13-22.

Long, J.S. & Freese, J. (2014). Regression models for categorical dependent variables using Stata (3rd ed.). Stata Press.

Mize, T.D., Doan, L., & Long, J.S. (2019). A general framework for comparing predictions and marginal effects across models. Sociological Methodology, 49(1), 152-189. doi:10.1177/0081175019852763

Pustejovsky, J.E. & 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.

Verbeke, G. & Lesaffre, E. (1997). The effect of misspecifying the random-effects distribution in linear mixed models for longitudinal data. Computational Statistics & Data Analysis, 23(4), 541-556.

Wasserstein, R.L., Schirm, A.L., & Lazar, N.A. (2019). Moving to a world beyond "p < 0.05". The American Statistician, 73(sup1), 1-19.

Examples

# ---- Single-model usage ------------------------------------------
fit <- lm(wellbeing_score ~ age + sex + smoking, data = sochealth)

# Default APA layout: B / SE / 95% CI / p plus the n / R^2 /
# Adj. R^2 fit-stats footer. Factor reference level is annotated
# with `(ref.)` and shows an en dash in the statistic columns.
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.

# \donttest{
# Standardized coefficients (beta) injected next to B. "refit"
# is the Cohen et al. (2003) refit-on-z-scores convention;
# "basic" reproduces the SPSS / Stata regress, beta definition.
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.
# }

# \donttest{
# Custom column set: B + AME + AME-specific p-value. Note that
# the `p` token always belongs to B, never to AME -- use the
# explicit `ame_p` token for AME inference.
if (requireNamespace("marginaleffects", quietly = TRUE)) {
  table_regression(
    fit,
    show_columns = c("b", "p", "ame", "ame_ci", "ame_p")
  )
}
#> Linear regression: wellbeing_score
#> 
#>  Variable        │    B       p     AME      95% CI        p   
#> ─────────────────┼─────────────────────────────────────────────
#>  (Intercept)     │   65.20  <.001                              
#>  age             │    0.05   .130   0.05  [-0.01, 0.11]   .130 
#>  sex:            │                                             
#>    Female (ref.) │     –     –       –          –         –    
#>    Male          │    3.86  <.001   3.86  [ 2.08, 5.63]  <.001 
#>  smoking:        │                                             
#>    No (ref.)     │     –     –       –          –         –    
#>    Yes           │   -1.72   .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.
# }

# \donttest{
# Group-token shortcut: "all_b" + "all_ame" expands to the full
# B / AME column families side by side.
if (requireNamespace("marginaleffects", quietly = TRUE)) {
  table_regression(fit, show_columns = c("all_b", "all_ame"))
}
#> Linear regression: wellbeing_score
#> 
#>  Variable        │    B      SE       95% CI        p     AME    SE  
#> ─────────────────┼───────────────────────────────────────────────────
#>  (Intercept)     │   65.20  1.66  [61.95, 68.45]  <.001              
#>  age             │    0.05  0.03  [-0.01,  0.11]   .130   0.05  0.03 
#>  sex:            │                                                   
#>    Female (ref.) │     –     –          –          –       –     –   
#>    Male          │    3.86  0.91  [ 2.08,  5.63]  <.001   3.86  0.91 
#>  smoking:        │                                                   
#>    No (ref.)     │     –     –          –          –       –     –   
#>    Yes           │   -1.72  1.11  [-3.89,  0.45]   .121  -1.72  1.11 
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  n               │ 1175                                              
#>  R²              │    0.02                                           
#>  Adj. R²         │    0.02                                           
#> 
#>  Variable        │ 95% CI (AME)   p (AME) 
#> ─────────────────┼────────────────────────
#>  (Intercept)     │                        
#>  age             │ [-0.01, 0.11]     .130 
#>  sex:            │                        
#>    Female (ref.) │       –           –    
#>    Male          │ [ 2.08, 5.63]    <.001 
#>  smoking:        │                        
#>    No (ref.)     │       –           –    
#>    Yes           │ [-3.89, 0.45]     .121 
#> 
#> Note. Linear regression. Std. errors: classical (OLS). AME of a continuous predictor: average slope (dY/dX). CI = confidence interval; AME = average marginal effect.
# }

# ---- Cluster-robust variance -------------------------------------
# CR2 (Bell-McCaffrey) with Satterthwaite-corrected df is the
# recommended default under few clusters. Three forms are accepted
# for `cluster`; the formula is preferred for composability with
# multi-way clustering and for programmatic robustness.
# \donttest{
if (requireNamespace("clubSandwich", quietly = TRUE)) {
  table_regression(fit, vcov = "CR2", cluster = ~region)
}
#> Linear regression: wellbeing_score
#> 
#>  Variable        │    B      SE       95% CI        p   
#> ─────────────────┼──────────────────────────────────────
#>  (Intercept)     │   65.20  1.76  [60.63, 69.78]  <.001 
#>  age             │    0.05  0.04  [-0.05,  0.15]   .285 
#>  sex:            │                                      
#>    Female (ref.) │     –     –          –          –    
#>    Male          │    3.86  0.85  [ 1.66,  6.05]   .007 
#>  smoking:        │                                      
#>    No (ref.)     │     –     –          –          –    
#>    Yes           │   -1.72  1.52  [-5.70,  2.27]   .313 
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  n               │ 1175                                 
#>  R²              │    0.02                              
#>  Adj. R²         │    0.02                              
#> 
#> Note. Linear regression. Std. errors: cluster-robust (CR2), clusters by region. CI = confidence interval.
if (requireNamespace("clubSandwich", quietly = TRUE)) {
  table_regression(fit, vcov = "CR2", cluster = "region")
}
#> Linear regression: wellbeing_score
#> 
#>  Variable        │    B      SE       95% CI        p   
#> ─────────────────┼──────────────────────────────────────
#>  (Intercept)     │   65.20  1.76  [60.63, 69.78]  <.001 
#>  age             │    0.05  0.04  [-0.05,  0.15]   .285 
#>  sex:            │                                      
#>    Female (ref.) │     –     –          –          –    
#>    Male          │    3.86  0.85  [ 1.66,  6.05]   .007 
#>  smoking:        │                                      
#>    No (ref.)     │     –     –          –          –    
#>    Yes           │   -1.72  1.52  [-5.70,  2.27]   .313 
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  n               │ 1175                                 
#>  R²              │    0.02                              
#>  Adj. R²         │    0.02                              
#> 
#> Note. Linear regression. Std. errors: cluster-robust (CR2), clusters by region. CI = confidence interval.
if (requireNamespace("clubSandwich", quietly = TRUE)) {
  table_regression(fit, vcov = "CR2", cluster = ~region:age_group)
}
#> Linear regression: wellbeing_score
#> 
#>  Variable        │    B      SE       95% CI        p   
#> ─────────────────┼──────────────────────────────────────
#>  (Intercept)     │   65.20  1.54  [61.91, 68.49]  <.001 
#>  age             │    0.05  0.03  [-0.01,  0.10]   .107 
#>  sex:            │                                      
#>    Female (ref.) │     –     –          –          –    
#>    Male          │    3.86  0.87  [ 2.05,  5.67]  <.001 
#>  smoking:        │                                      
#>    No (ref.)     │     –     –          –          –    
#>    Yes           │   -1.72  1.25  [-4.31,  0.88]   .183 
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  n               │ 1175                                 
#>  R²              │    0.02                              
#>  Adj. R²         │    0.02                              
#> 
#> Note. Linear regression. Std. errors: cluster-robust (CR2), clusters by region:age_group. CI = confidence interval.
# }

# \donttest{
# ---- Hierarchical (nested) regression ----------------------------
# Adds in-table change-statistic rows (Delta R^2 / F-change /
# p-change for lm; LRT / p-change for glm) below the fit-stats.
# Note: hierarchical comparison requires identical observations
# across all models -- prepare a complete-case subset first so
# R's listwise deletion does not produce different `nobs` per
# model (which the function rejects).
sochealth_cc <- na.omit(
  sochealth[, c("wellbeing_score", "age", "sex", "smoking")]
)
m1 <- lm(wellbeing_score ~ age,                  data = sochealth_cc)
m2 <- lm(wellbeing_score ~ age + sex,            data = sochealth_cc)
m3 <- lm(wellbeing_score ~ age + sex + smoking,  data = sochealth_cc)
table_regression(
  list("Step 1" = m1, "Step 2" = m2, "Step 3" = m3),
  nested = TRUE
)
#> Hierarchical linear regression: wellbeing_score
#> 
#>                           Step 1                     Step 2             
#>                    ────────────────────  ────────────────────────────── 
#>  Variable        │    B      SE     p            B           SE     p   
#> ─────────────────┼──────────────────────────────────────────────────────
#>  (Intercept)     │   66.85  1.59  <.001    64.90            1.65  <.001 
#>  age             │    0.04  0.03   .160     0.05            0.03   .143 
#>  sex:            │                                                      
#>    Female (ref.) │                           –               –     –    
#>    Male          │                          3.87            0.91  <.001 
#>  smoking:        │                                                      
#>    No (ref.)     │                                                      
#>    Yes           │                                                      
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  n               │ 1175                  1175                           
#>  R²              │    0.00                  0.02                        
#>  Adj. R²         │    0.00                  0.02                        
#>  ΔR²             │     –                   +0.02                        
#>  F-change        │     –                   18.26 (1, 1172)              
#>  p (change)      │     –                    <.001                       
#> 
#>                                Step 3             
#>                    ────────────────────────────── 
#>  Variable        │         B           SE     p   
#> ─────────────────┼────────────────────────────────
#>  (Intercept)     │   65.20            1.66  <.001 
#>  age             │    0.05            0.03   .130 
#>  sex:            │                                
#>    Female (ref.) │     –               –     –    
#>    Male          │    3.86            0.91  <.001 
#>  smoking:        │                                
#>    No (ref.)     │     –               –     –    
#>    Yes           │   -1.72            1.11   .121 
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  n               │ 1175                           
#>  R²              │    0.02                        
#>  Adj. R²         │    0.02                        
#>  ΔR²             │   +0.00                        
#>  F-change        │    2.41 (1, 1171)              
#>  p (change)      │     .121                       
#> 
#> Note. Linear regression models. Std. errors: classical (OLS).

# The same comparison by its fit statistics alone.
table_regression(
  list("Step 1" = m1, "Step 2" = m2, "Step 3" = m3),
  show_columns = "none",
  nested = TRUE
)
#> Hierarchical linear regression: wellbeing_score
#> 
#>  Variable   │ Step 1        Step 2             Step 3       
#> ────────────┼───────────────────────────────────────────────
#>  n          │ 1175     1175               1175              
#>  R²         │    0.00     0.02               0.02           
#>  Adj. R²    │    0.00     0.02               0.02           
#>  ΔR²        │     –      +0.02              +0.00           
#>  F-change   │     –      18.26 (1, 1172)     2.41 (1, 1171) 
#>  p (change) │     –       <.001               .121          
#> 
#> Note. Linear regression models.
# }

# \donttest{
# ---- Side-by-side variance comparison ----------------------------
# Same fit, three vcovs in one wide table. Useful for showing the
# sensitivity of inference to the variance assumption.
if (requireNamespace("clubSandwich", quietly = TRUE)) {
  table_regression(
    list("Classical" = fit, "HC3" = fit, "CR2" = fit),
    vcov    = list("classical", "HC3", "CR2"),
    cluster = list(NULL, NULL, ~region)
  )
}
#> Linear regression comparison: wellbeing_score
#> 
#>                         Classical                HC3                CR2      
#>                    ────────────────────  ────────────────────  ───────────── 
#>  Variable        │    B      SE     p       B      SE     p       B      SE  
#> ─────────────────┼───────────────────────────────────────────────────────────
#>  (Intercept)     │   65.20  1.66  <.001    65.20  1.61  <.001    65.20  1.76 
#>  age             │    0.05  0.03   .130     0.05  0.03   .127     0.05  0.04 
#>  sex:            │                                                           
#>    Female (ref.) │     –     –     –         –     –     –         –     –   
#>    Male          │    3.86  0.91  <.001     3.86  0.91  <.001     3.86  0.85 
#>  smoking:        │                                                           
#>    No (ref.)     │     –     –     –         –     –     –         –     –   
#>    Yes           │   -1.72  1.11   .121    -1.72  1.11   .123    -1.72  1.52 
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  n               │ 1175                  1175                  1175          
#>  R²              │    0.02                  0.02                  0.02       
#>  Adj. R²         │    0.02                  0.02                  0.02       
#> 
#>                     CR2  
#>                    ───── 
#>  Variable        │ p (B) 
#> ─────────────────┼───────
#>  (Intercept)     │ <.001 
#>  age             │  .285 
#>  sex:            │       
#>    Female (ref.) │  –    
#>    Male          │  .007 
#>  smoking:        │       
#>    No (ref.)     │  –    
#>    Yes           │  .313 
#> 
#> Note. Linear regression models. Std. errors: Classical: classical (OLS); HC3: heteroskedasticity-robust (HC3); CR2: cluster-robust (CR2), clusters by region.
# }

# ---- Tidy long format for downstream pipelines -------------------
if (requireNamespace("broom", quietly = TRUE)) {
  broom::tidy(table_regression(fit))
}
#> # A tibble: 4 × 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            age   B               0.0465    0.0307  -0.0137
#> 3 M1       wellbe… NA            sexM… B               3.86      0.905    2.08  
#> 4 M1       wellbe… NA            smok… B              -1.72      1.11    -3.89  
#> # ℹ 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>

# \donttest{
# ---- Mixed-effects models ----------------------------------------
# Linear mixed-effects (lme4). The footer adds a random-effects
# block with sigma + Wald SE / CI from `merDeriv`, the Nakagawa
# marginal / conditional R^2 fit-stats, and a per-class p-value
# annotation line.
if (requireNamespace("lme4", quietly = TRUE)) {
  fit <- lme4::lmer(Reaction ~ Days + (Days | Subject),
                     data = lme4::sleepstudy)
  table_regression(fit)

  # Switch to the variance scale (sigma^2 instead of sigma).
  table_regression(fit, re_scale = "variance")

  # Minimal random-effects display: estimates only, no SE / CI.
  table_regression(fit, re_columns = "est")

  # Suppress the random-effects block entirely.
  table_regression(fit, show_re = FALSE)
}
#> P-values of the `lme4::lmer()` model are Wald z-tests, a large-sample approximation.
#> ℹ For Satterthwaite t-tests, which suit models with few groups, refit with `lmerTest::lmer()` (install lmerTest if needed).
#> This message is displayed once per session.
#> Registered S3 method overwritten by 'merDeriv':
#>   method        from        
#>   bread.lmerMod clubSandwich
#> 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 
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  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. CI = confidence interval; AIC = Akaike information criterion; BIC = Bayesian information criterion.

# Hierarchical mixed-effects comparison (nested LRT).
if (requireNamespace("lme4", quietly = TRUE)) {
  m1 <- lme4::lmer(Reaction ~ 1     + (1 | Subject),
                    data = lme4::sleepstudy, REML = FALSE)
  m2 <- lme4::lmer(Reaction ~ Days  + (1 | Subject),
                    data = lme4::sleepstudy, REML = FALSE)
  m3 <- lme4::lmer(Reaction ~ Days  + (Days | Subject),
                    data = lme4::sleepstudy, REML = FALSE)
  table_regression(list(m1, m2, m3), nested = TRUE)
}
#> Hierarchical linear mixed-effects regression: Reaction
#> 
#>                                          Model 1              Model 2      
#>                                    ────────────────────  ───────────────── 
#>  Variable                        │    B      SE     p         B        SE  
#> ─────────────────────────────────┼─────────────────────────────────────────
#>  (Intercept)                     │  298.51  8.79  <.001   251.41      9.51 
#>  Days                            │                         10.47      0.80 
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  Random effects:                 │                                         
#>    σ Subject (Intercept)         │   34.59  6.72   –       36.01      6.45 
#>    σ Subject Days                │                                         
#>    ρ Subject ((Intercept), Days) │                                         
#>    σ (Residual)                  │   44.26  2.46   –       30.90      1.72 
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  n                               │  180                   180              
#>  N (Subject)                     │   18                    18              
#>  ICC                             │    0.38                  0.58           
#>  R² (marginal)                   │    0.00                  0.29           
#>  R² (conditional)                │    0.38                  0.70           
#>  AIC                             │ 1916.5                1802.1            
#>  BIC                             │ 1926.1                1814.9            
#>  ΔAIC                            │     –                 -114.5            
#>  ΔBIC                            │     –                 -111.3            
#>  Δχ²                             │     –                  116.46 (1)       
#>  p (change)                      │     –                    <.001          
#> 
#>                                    Model 2          Model 3          
#>                                    ───────  ──────────────────────── 
#>  Variable                        │   p (B)       B        SE     p   
#> ─────────────────────────────────┼───────────────────────────────────
#>  (Intercept)                     │   <.001   251.41      6.63  <.001 
#>  Days                            │   <.001    10.47      1.50  <.001 
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  Random effects:                 │                                   
#>    σ Subject (Intercept)         │    –       23.78      5.58   –    
#>    σ Subject Days                │             5.72      1.19   –    
#>    ρ Subject ((Intercept), Days) │             0.08      0.32   –    
#>    σ (Residual)                  │    –       25.59      1.51   –    
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#>  n                               │           180                     
#>  N (Subject)                     │            18                     
#>  ICC                             │              –                    
#>  R² (marginal)                   │             0.29                  
#>  R² (conditional)                │             0.79                  
#>  AIC                             │          1763.9                   
#>  BIC                             │          1783.1                   
#>  ΔAIC                            │           -38.1                   
#>  ΔBIC                            │           -31.8                   
#>  Δχ²                             │            42.14 (2)              
#>  p (change)                      │             <.001                 
#> 
#> Note. Linear mixed-effects regression models. Std. errors: Wald (model-based). P-values: Wald z-test, large-sample approximation. Model 1: Random effects (ML): LR test vs linear regression, χ̄²(1) = 50.51, p < .001. Model 2: Random effects (ML): LR test vs linear regression, χ̄²(1) = 106.21, p < .001. Model 3: Random effects (ML): LR test vs linear regression, χ̄²(3) = 148.35, p < .001. ICC = intraclass correlation coefficient; AIC = Akaike information criterion; BIC = Bayesian information criterion.
# }

if (FALSE) { # \dontrun{
# ---- Rich-format outputs (require optional Suggests packages) ----
table_regression(fit, output = "gt")
table_regression(fit, output = "flextable")
table_regression(fit, output = "tinytable")

# ---- File outputs ------------------------------------------------
table_regression(fit, output = "excel",
                 excel_path = tempfile(fileext = ".xlsx"))
table_regression(fit, output = "word",
                 word_path  = tempfile(fileext = ".docx"))

# ---- System clipboard (interactive use) --------------------------
table_regression(fit, output = "clipboard")
} # }