Regression Tables

Univariable Models

uni_reg() fits one model per exposure and returns a table ready for reports. Variable labels set with attr(x, "label") or labelled::var_label() are used automatically in the displayed table, while raw column names remain in $table_body.

birthwt_uni <- uni_reg(
  data = birthwt_data,
  outcome = "low",
  exposures = birthwt_exposures,
  approach = "logit",
  theme = clinical
)

birthwt_uni$table

Multivariable Models

With adjust_for = NULL, multi_reg() fits all supplied exposures in one multivariable model. This is the usual fully adjusted model when every exposure listed should appear in the same formula.

birthwt_full <- multi_reg(
  data = birthwt_data,
  outcome = low,
  exposures = c("age", "lwt", "race", "smoke", "ht", "ui"),
  approach = logit,
  theme = clinical
)

birthwt_full$table

Exposure-Specific Adjusted Models

Use adjust_for when you want one adjusted model per exposure, each adjusted for the same core covariate set. This is useful for screening several clinically important exposures while keeping the adjustment strategy explicit.

birthwt_multi <- multi_reg(
  data = birthwt_data,
  outcome = low,
  exposures = c("smoke", "ht", "ui", "ptl_cat", "ftv_cat"),
  adjust_for = c("age", "lwt", "race"),
  approach = logit,
  theme = striped
)

birthwt_multi$table

The adjustment variables are recorded in a compact table footnote, so the result is ready for manuscript-style reporting without making the table unnecessarily tall.

Optional Model Statistics

Publication tables should stay readable. When you need model-fit information, set model_stats = TRUE and inspect the returned object’s $model_stats element. This keeps AIC, BIC, log-likelihood, deviance, pseudo R-squared, and linear-model R-squared values available without adding clutter to the main table.

birthwt_uni_stats <- uni_reg(
  data = birthwt_data,
  outcome = low,
  exposures = birthwt_exposures,
  approach = logit,
  model_stats = TRUE
)

birthwt_uni_stats$model_stats

##     model      AIC      BIC    logLik deviance null_deviance  pseudo_r2
## 1     age 235.9120 242.3955 -115.9560 231.9120       234.672 0.01176126
## 2     lwt 232.6907 239.1742 -114.3453 228.6907       234.672 0.02548803
## 3    race 235.6616 245.3869 -114.8308 229.6616       234.672 0.02135051
## 4   smoke 233.8046 240.2881 -114.9023 229.8046       234.672 0.02074128
## 5      ht 234.6499 241.1334 -115.3249 230.6499       234.672 0.01713938
## 6      ui 233.5959 240.0794 -114.7979 229.5959       234.672 0.02163060
## 7 ptl_cat 225.8978 232.3812 -110.9489 221.8978       234.672 0.05443445
## 8 ftv_cat 238.0851 247.8103 -116.0425 232.0851       234.672 0.01102346
##   r_squared adj_r_squared   n
## 1        NA            NA 189
## 2        NA            NA 189
## 3        NA            NA 189
## 4        NA            NA 189
## 5        NA            NA 189
## 6        NA            NA 189
## 7        NA            NA 189
## 8        NA            NA 189

For adjusted mode, multi_reg() returns one row of statistics per adjusted exposure-specific model.

birthwt_multi_stats <- multi_reg(
  data = birthwt_data,
  outcome = low,
  exposures = c("smoke", "ht", "ui", "ptl_cat"),
  adjust_for = c("age", "lwt", "race"),
  approach = logit,
  model_stats = TRUE
)

birthwt_multi_stats$model_stats

##     model      AIC      BIC    logLik deviance null_deviance  pseudo_r2
## 1   smoke 226.5772 246.0277 -107.2886 214.5772       234.672 0.08562914
## 2      ht 227.7487 247.1991 -107.8743 215.7487       234.672 0.08063739
## 3      ui 231.0194 250.4699 -109.5097 219.0194       234.672 0.06669985
## 4 ptl_cat 222.4421 241.8926 -105.2211 210.4421       234.672 0.10324998
##   r_squared adj_r_squared   n
## 1        NA            NA 189
## 2        NA            NA 189
## 3        NA            NA 189
## 4        NA            NA 189

Clustered Data and Multilevel Models

Use multilevel_reg() when observations are grouped and responses from the same group may be correlated. For example, patients may be clustered within clinics. The function fits one random intercept for the grouping variable and reports publication-ready fixed effects. The grouping variable defines the correlation structure; it does not receive its own effect estimate.

set.seed(71)
clinic_data <- data.frame(
  clinic = factor(rep(seq_len(20), each = 12)),
  treatment = factor(rep(c("Control", "Intervention"), 120)),
  age = rnorm(240, 58, 9)
)
clinic_effect <- rnorm(20, 0, 3)
clinic_data$score <- 50 +
  2 * (clinic_data$treatment == "Intervention") +
  0.2 * clinic_data$age +
  clinic_effect[clinic_data$clinic] +
  rnorm(240)

clustered_model <- multilevel_reg(
  clinic_data,
  outcome = score,
  exposures = c(treatment, age),
  group = clinic,
  approach = linear,
  model_stats = TRUE
)

clustered_model
clustered_model$model_stats

For a continuous outcome, the table reports adjusted beta coefficients and the ICC is calculated on the observed scale. With approach = "logit", the table reports adjusted odds ratios and the ICC uses the usual latent logistic scale. The fixed-effect p-values use a Wald normal approximation. A singular-fit warning indicates that the estimated between-group variance is near zero; it is a prompt to reconsider whether the random intercept is supported by the data, not a result to hide.

The focused API deliberately supports one random intercept. More complex random slopes, crossed effects, or nested structures require explicit model specification and specialist interpretation outside this workflow.

Ordered Outcomes

Use ordinal_reg() when the outcome has at least three categories with a meaningful order. Create the outcome as an ordered factor and place its levels from lower to higher. The package deliberately does not infer this order from alphabetical labels.

recovery_data$recovery <- ordered(
  recovery_data$recovery,
  levels = c("Poor", "Fair", "Good", "Excellent")
)

recovery_model <- ordinal_reg(
  recovery_data,
  outcome = recovery,
  exposures = c(treatment, age),
  adjust_for = sex,
  model_stats = TRUE,
  assumption_test = TRUE
)

recovery_model
recovery_model$model_stats
recovery_model$assumption_tests

The table reports common odds ratios. An odds ratio above one indicates greater odds of being in a higher outcome category, assuming proportional odds. The outcome order is repeated in the table footnote so the direction remains clear after export. $assumption_tests provides a term-level check: a small p-value suggests that a single common odds ratio may not adequately describe that predictor across every outcome threshold.

Unordered Outcomes

Use multinomial_reg() when an outcome has at least three categories without a natural order. The first factor level is the outcome reference category. Set it deliberately with factor() or change it reproducibly with relevel().

discharge_data$destination <- factor(
  discharge_data$destination,
  levels = c("Home", "Rehabilitation", "Hospital")
)

destination_model <- multinomial_reg(
  discharge_data,
  outcome = destination,
  exposures = c(treatment, age),
  adjust_for = sex,
  model_stats = TRUE
)

destination_model
destination_model$model_stats
names(destination_model$per_outcome)

The main publication table presents each comparison separately, such as Rehabilitation versus Home and Hospital versus Home. Odds ratios from different outcome comparisons should not be collapsed into one estimate. For that reason, plots use one named comparison at a time:

rehabilitation_model <- destination_model$per_outcome[["Rehabilitation"]]
plot_reg(rehabilitation_model)
forest_reg(forest_df(rehabilitation_model))

Other Effect Measures

Switch the approach to change the estimand.

birthwt_rr <- uni_reg(
  data = birthwt_data,
  outcome = low,
  exposures = c("smoke", "ht", "ui", "ptl_cat"),
  approach = logbinomial
)

birthwt_rr$table

Use approach = firth when a binary-outcome logistic model has sparse cells, very wide intervals, or separation concerns. The output remains an odds-ratio table, but the model is fitted with Firth penalized logistic regression.

The built-in data_endometrial dataset is a useful teaching example because neovascularization is completely absent among low-grade cases, a pattern that can make ordinary logistic regression unstable.

data("data_endometrial", package = "gtregression")

endometrial_data <- data_endometrial |>
  mutate(
    HG = factor(HG, levels = c(0, 1),
                labels = c("Low grade", "High grade")),
    NV = factor(NV, levels = c(0, 1), labels = c("Absent", "Present"))
  )

endometrial_firth <- multi_reg(
  data = endometrial_data,
  outcome = HG,
  exposures = c(NV, PI, EH),
  approach = firth
)

## logistf::logistf(formula = formula, data = data)
## 
## Model fitted by Penalized ML
## Coefficients:
##                    coef   se(coef) lower 0.95  upper 0.95      Chisq
## (Intercept)  3.77455951 1.43900672  1.0825371  7.20928050  8.1980136
## NVPresent    2.92927330 1.46497415  0.6097244  7.85463171  6.7984572
## PI          -0.03475175 0.03789237 -0.1244587  0.04045547  0.7468285
## EH          -2.60416387 0.75362838 -4.3651832 -1.23272106 17.7593175
##                        p method
## (Intercept) 4.193628e-03      2
## NVPresent   9.123668e-03      2
## PI          3.874822e-01      2
## EH          2.506867e-05      2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=43.65582 on 3 df, p=1.78586e-09, n=79
## Wald test = 21.66965 on 3 df, p = 7.641345e-05

endometrial_firth$table

Cox Survival Models

cox_reg() uses direct time and event arguments and returns hazard ratios. Without adjust_for, the table shows crude HRs. With adjust_for, the table shows adjusted HRs.

data("data_lungcancer", package = "gtregression")

lung_data <- data_lungcancer |>
  mutate(
    trt = factor(trt, levels = c(1, 2),
                 labels = c("Standard treatment", "Test treatment")),
    prior = factor(prior, levels = c(0, 10), labels = c("No", "Yes"))
  )

attr(lung_data$trt, "label") <- "Treatment group"
attr(lung_data$celltype, "label") <- "Cancer cell type"
attr(lung_data$karno, "label") <- "Karnofsky performance score"
attr(lung_data$age, "label") <- "Age"
attr(lung_data$prior, "label") <- "Prior therapy"

lung_hr <- cox_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = c("trt", "celltype", "karno", "age"),
  theme = clinical
)

lung_hr$table

lung_adj_hr <- cox_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = c(trt, celltype, prior),
  adjust_for = c(age, karno),
  model_stats = TRUE,
  theme = striped
)

lung_adj_hr$table

lung_adj_hr$model_stats

##      model      AIC      BIC    logLik concordance events   n
## 1      trt 973.7560 982.3121 -483.8780   0.7119491    128 137
## 2 celltype 961.0882 975.3484 -475.5441   0.7349500    128 137
## 3    prior 974.7459 983.3019 -484.3729   0.7134257    128 137

Use baseline_strata when a categorical variable should have a separate baseline hazard but should not receive its own hazard ratio. This fits one Cox model with common predictor effects. It is different from stratifier, which fits separate models and produces stratum-specific estimates.

lung_baseline_strata <- cox_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = c(trt, age, karno),
  baseline_strata = celltype,
  multivariable = TRUE
)

lung_baseline_strata$table

Use interaction = exposure*modifier for planned effect modification. In the default exposure-by-exposure workflow, provide one exposure.

cox_interaction <- cox_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = trt,
  adjust_for = c(age, karno),
  interaction = trt*prior
)

cox_interaction$table

Parametric Survival Models

surv_reg() uses the same time, event, exposures, and adjust_for grammar as cox_reg(), but fits parametric survival models with survival::survreg(). The table reports time ratios rather than hazard ratios. A time ratio above 1 suggests longer survival time; below 1 suggests shorter survival time, conditional on the chosen distribution.

lung_time_ratio <- surv_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = c("trt", "celltype", "karno", "age"),
  distribution = weibull,
  theme = clinical
)

lung_time_ratio$table

lung_adj_time_ratio <- surv_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = c(trt, celltype, prior),
  adjust_for = c(age, karno),
  distribution = lognormal,
  model_stats = TRUE,
  theme = striped
)

lung_adj_time_ratio$table

lung_adj_time_ratio$model_stats

##      model distribution      AIC      BIC    logLik    scale events   n
## 1      trt    lognormal 1451.264 1465.864 -720.6318 1.110644    128 137
## 2 celltype    lognormal 1444.379 1464.819 -715.1894 1.064351    128 137
## 3    prior    lognormal 1451.768 1466.368 -720.8840 1.112741    128 137

surv_interaction <- surv_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = trt,
  adjust_for = c(age, karno),
  interaction = trt*prior,
  distribution = weibull
)

surv_interaction$table

Continuous Outcomes

Linear regression outputs beta coefficients and keeps diagnostics under $reg_check.

birthwt_linear <- multi_reg(
  data = birthwt_data,
  outcome = bwt,
  exposures = c("age", "lwt", "race", "smoke", "ht", "ui"),
  approach = linear
)

birthwt_linear$table

What To Inspect

• $table: publication-ready table.
• $table_body: numeric estimates behind the display.
• $models: fitted model objects.
• $model_summaries: model-level summaries.
• $model_stats: optional model-fit statistics when model_stats = TRUE.
• $variable_labels: labels used in publication output.
• $reg_check: diagnostics for linear models.