---
title: "Regression Tables"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Regression Tables}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

Go from fitted models to publication-ready tables without hand-formatting
effect estimates. `gtregression` supports logistic, log-binomial, Poisson,
robust Poisson, negative binomial, Cox survival, parametric survival, and
linear regression.

```{r reg-setup, message=FALSE, warning=FALSE}
library(gtregression)
library(dplyr)

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

birthwt_data <- data_birthwt |>
  mutate(
    race = factor(race, levels = c(1, 2, 3),
                  labels = c("White", "Black", "Other")),
    smoke = factor(smoke, levels = c(0, 1), labels = c("No", "Yes")),
    ht = factor(ht, levels = c(0, 1), labels = c("No", "Yes")),
    ui = factor(ui, levels = c(0, 1), labels = c("No", "Yes")),
    low = factor(low, levels = c(0, 1), labels = c("Normal BW", "Low BW")),
    ptl_cat = factor(ifelse(ptl > 0, "Yes", "No"), levels = c("No", "Yes")),
    ftv_cat = factor(case_when(
      ftv == 0 ~ "None",
      ftv == 1 ~ "One",
      ftv >= 2 ~ "Two or more"
    ), levels = c("None", "One", "Two or more"))
  )

birthwt_exposures <- c(
  "age", "lwt", "race", "smoke", "ht", "ui", "ptl_cat", "ftv_cat"
)

attr(birthwt_data$age, "label") <- "Maternal age"
attr(birthwt_data$lwt, "label") <- "Maternal weight"
attr(birthwt_data$race, "label") <- "Maternal race"
attr(birthwt_data$smoke, "label") <- "Smoking during pregnancy"
attr(birthwt_data$ht, "label") <- "Hypertension"
attr(birthwt_data$ui, "label") <- "Uterine irritability"
attr(birthwt_data$ptl_cat, "label") <- "Previous preterm labour"
attr(birthwt_data$ftv_cat, "label") <- "First trimester visits"
```

## 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`.

```{r uni-logit, message=FALSE, warning=FALSE}
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.

```{r multi-full, message=FALSE, warning=FALSE}
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.

```{r multi-adjusted, message=FALSE, warning=FALSE}
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.

```{r model-stats, message=FALSE, warning=FALSE}
birthwt_uni_stats <- uni_reg(
  data = birthwt_data,
  outcome = low,
  exposures = birthwt_exposures,
  approach = logit,
  model_stats = TRUE
)

birthwt_uni_stats$model_stats
```

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

```{r adjusted-model-stats, message=FALSE, warning=FALSE}
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
```

## 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.

```{r multilevel-example, eval=FALSE}
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.

```{r ordinal-example, eval=FALSE}
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()`.

```{r multinomial-example, eval=FALSE}
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:

```{r multinomial-plot-example, eval=FALSE}
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.

```{r reg-risk-ratio, message=FALSE, warning=FALSE}
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.

```{r reg-firth, message=FALSE, warning=FALSE}
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
)

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.

```{r cox-setup, message=FALSE, warning=FALSE}
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"
```

```{r cox-crude, message=FALSE, warning=FALSE}
lung_hr <- cox_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = c("trt", "celltype", "karno", "age"),
  theme = clinical
)

lung_hr$table
```

```{r cox-adjusted, message=FALSE, warning=FALSE}
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
```

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.

```{r cox-baseline-strata, message=FALSE, warning=FALSE}
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.

```{r cox-interaction, message=FALSE, warning=FALSE}
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.

```{r surv-crude, message=FALSE, warning=FALSE}
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
```

```{r surv-adjusted, message=FALSE, warning=FALSE}
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
```

```{r surv-interaction, message=FALSE, warning=FALSE}
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`.

```{r reg-linear, message=FALSE, warning=FALSE}
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.