---
title: "Diagnostics and Model Selection"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Diagnostics and Model Selection}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

Before the final table goes into a manuscript, check the model. These helpers
keep diagnostics close to the regression workflow. They are screening aids:
interpret them with the study design, clinical or subject-matter judgement, and
the diagnostics from the fitted model.

```{r diag-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"))
  )

exposures <- c("age", "lwt", "race", "smoke", "ht", "ui", "ptl_cat")
```

## Convergence Screening

Use `check_convergence()` before interpreting model estimates, especially for
log-binomial and small or sparse binary-outcome models. A non-converged model is
a fitting warning, not a finding.

```{r diag-convergence, message=FALSE, warning=FALSE}
check_convergence(
  data = birthwt_data,
  exposures = exposures,
  outcome = low,
  approach = logit,
  multivariate = TRUE
)
```

For risk-ratio workflows, this same check helps users decide whether a
log-binomial model fitted cleanly or whether a robust Poisson approach may be a
more practical sensitivity analysis.

```{r diag-convergence-logbin, message=FALSE, warning=FALSE}
check_convergence(
  data = birthwt_data,
  exposures = c("smoke", "ht", "ui", "ptl_cat"),
  outcome = low,
  approach = logbinomial,
  multivariate = TRUE
)
```

## Collinearity Screening

`check_collinearity()` reports VIF-style diagnostics for multivariable models.
High VIF values are prompts to inspect coding, overlap between predictors, and
the scientific purpose of the model.

```{r diag-collinearity, message=FALSE, warning=FALSE}
birthwt_multi <- multi_reg(
  data = birthwt_data,
  outcome = low,
  exposures = exposures,
  approach = logit
)

check_collinearity(birthwt_multi, format = gt)
```

Adjusted-mode `multi_reg()` objects contain one model per exposure. The
collinearity output keeps that list structure so each model can be inspected
separately.

```{r diag-collinearity-adjusted, message=FALSE, warning=FALSE}
birthwt_adjusted <- multi_reg(
  data = birthwt_data,
  outcome = low,
  exposures = c("smoke", "ht", "ui", "ptl_cat"),
  adjust_for = c("age", "lwt", "race"),
  approach = logit
)

check_collinearity(birthwt_adjusted, format = tibble)
```

## Model Fit Plots

`plot_model_fit()` turns fitted models into quick diagnostic plots. It accepts
raw `lm()` and `glm()` objects, and it also works with models saved inside
`uni_reg()` and `multi_reg()` results.

For logistic regression, the calibration plot compares predicted probabilities
with observed event proportions. Points close to the diagonal line suggest that
the model predictions are reasonably aligned with the observed data. Calibration
is usually most useful for multivariable models because predicted probabilities
vary across many people.

```{r diag-fit-logistic, message=FALSE, warning=FALSE}
plot_model_fit(
  birthwt_multi,
  type = calibration,
  bins = 6
)
```

When a `uni_reg()` object contains several models, use `model_name` to choose
the exposure you want to inspect. For a simple binary exposure, calibration may
only show two points because the model has only two fitted probabilities; in
that situation, residual and influence plots are usually more useful.

```{r diag-fit-uni, message=FALSE, warning=FALSE}
birthwt_uni <- uni_reg(
  data = birthwt_data,
  outcome = low,
  exposures = c("age", "lwt", "smoke"),
  approach = logit
)

plot_model_fit(
  birthwt_uni,
  model_name = smoke,
  type = residual
)
```

For logistic models, residual plots often form two visible bands. That is a
normal consequence of a 0/1 outcome and should be interpreted as a screening
plot rather than a linear-model residual plot.

For linear regression, `type = all` shows the classic residual, Q-Q,
scale-location, and Cook's distance views.

```{r diag-fit-linear, message=FALSE, warning=FALSE}
fit_lm <- lm(bwt ~ age + lwt, data = birthwt_data)
plot_model_fit(fit_lm)
```

## Proportional Hazards Screening

For Cox models, use `check_ph()` before treating hazard ratios as final. It
reports Schoenfeld residual tests from `survival::cox.zph()`, including a global
test. Small p-values suggest possible non-proportional hazards and should be
reviewed with plots, follow-up pattern, and clinical judgement.

```{r diag-ph-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")),
    celltype = factor(
      celltype,
      levels = c("squamous", "smallcell", "adeno", "large"),
      labels = c("Squamous", "Small cell", "Adenocarcinoma", "Large cell")
    )
  )

cox_fit <- cox_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = c(trt, celltype, prior),
  adjust_for = c(age, karno)
)
```

```{r diag-ph-table, message=FALSE, warning=FALSE}
check_ph(cox_fit, format = gt)
```

Use `format = tibble` when you want to inspect or filter the diagnostic results.

```{r diag-ph-tibble, message=FALSE, warning=FALSE}
check_ph(cox_fit, transform = rank, format = tibble)
```

## Stepwise Model Selection

`compare_models()` is for prespecified candidate models that have already been
fitted with gtregression. It answers a different question from stepwise
selection: "How do these planned models compare?" The inputs should be
`multi_reg()`, `cox_reg()`, or `surv_reg()` outputs, not raw `lm()`, `glm()`,
`coxph()`, or `survreg()` objects. This keeps the workflow consistent with the
publication-ready tables created by the package.

```{r diag-compare-logit, message=FALSE, warning=FALSE}
logit_m0 <- multi_reg(
  data = birthwt_data,
  outcome = low,
  exposures = smoke,
  approach = logit
)

logit_m1 <- multi_reg(
  data = birthwt_data,
  outcome = low,
  exposures = c(smoke, age, lwt),
  approach = logit
)

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

logit_model_comparison <- compare_models(
  logit_m0,
  logit_m1,
  logit_m2,
  model_names = c(
    "Smoking only",
    "Add age and weight",
    "Full clinical model"
  ),
  primary_exposure = smoke,
  format = gt
)

logit_model_comparison$table
```

The table reports N, number of parameters, AIC, BIC, log-likelihood, and
likelihood-ratio comparisons when `nested = TRUE`. Lower AIC or BIC identifies
better relative fit among the compared models. When `primary_exposure` is
supplied, the table also tracks that effect estimate and the percentage change
across models.

`compare_models()` automatically checks whether the candidate models appear to
use the same analysis sample. It uses retained row identifiers when the fitted
model stores them; otherwise it compares N and event counts. If the models use
different complete-case samples, the table still displays AIC, BIC,
log-likelihood, and likelihood-ratio statistics for transparency, but the footer
warns that these values should not be interpreted as formal model-selection
criteria across different datasets. In that situation, use the primary exposure
estimate, percentage change, confidence intervals, and clinical or
epidemiological reasoning to judge robustness.

### Build Candidates in the App

The function API intentionally accepts fitted gtregression objects. In the
gtregression app, **Advanced > Compare models** provides a guided layer before
that API: choose the shared outcome and approach, name two to six candidates,
and select each candidate's exposures, adjustment variables, and optional
interaction. The app fits `multi_reg()`, `cox_reg()`, or `surv_reg()` objects and
then calls `compare_models()` with those objects.

For the birth weight models above, enter *Smoking only*, *Add age and weight*,
and *Full clinical model* as the candidate names, place `smoke` in the exposure
control, and add the corresponding variables under adjustment. Select `smoke`
as the primary exposure to reproduce the estimate-change columns. The generated
code records every fitted model and is suitable for direct use in an R script.

Candidate order is deliberate. Put the simplest or primary model first because
percentage change is calculated from that model and likelihood-ratio statistics
compare sequential pairs. If the fitted samples differ, or a sequential pair is
not nested, the rendered table displays only the warning relevant to that
comparison. See the **gtregression App** article for the complete ordinary,
Cox, and parametric survival walkthrough.

For Cox and parametric survival models, fit the candidate models with `cox_reg()`
or `surv_reg()` first. `compare_models()` then keeps survival-specific columns
such as events and Cox concordance.

```{r diag-compare-cox, message=FALSE, warning=FALSE}
cox_m0 <- cox_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = trt
)

cox_m1 <- cox_reg(
  data = lung_data,
  time = time,
  event = status,
  exposures = trt,
  adjust_for = c(age, karno)
)

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

cox_model_comparison <- compare_models(
  list(
    "Treatment only" = cox_m0,
    "Add age and performance" = cox_m1,
    "Full clinical model" = cox_m2
  ),
  primary_exposure = trt,
  format = gt
)

cox_model_comparison$table
```

`select_models()` compares candidate models step by step. It is useful for
exploration, teaching, and sensitivity checks. It should not replace a planned
model based on study design or a causal framework.

```{r diag-select, message=FALSE, warning=FALSE}
selected <- select_models(
  data = birthwt_data,
  outcome = low,
  exposures = exposures,
  approach = logit,
  direction = forward
)

selected
```

The selected direction is recorded in the formatted table footer. Backward and
both-direction searches are available using the same interface.

```{r diag-select-other-directions, message=FALSE, warning=FALSE}
select_models(
  data = birthwt_data,
  outcome = low,
  exposures = exposures,
  approach = logit,
  direction = backward,
  format = tibble
)$results_table

select_models(
  data = birthwt_data,
  outcome = low,
  exposures = exposures,
  approach = logit,
  direction = both,
  format = tibble
)$results_table
```

## What To Inspect

- `check_convergence()`: convergence status and maximum fitted probabilities.
  Use `format = gt` or `format = flextable` for viewing tables.
- `check_collinearity()`: VIF and interpretation. Nested model outputs keep
  their list structure when formatted.
- `plot_model_fit()`: residual, calibration, observed-versus-predicted, and
  influence plots for `lm`/`glm` models and stored `uni_reg()` / `multi_reg()`
  fitted models.
- `check_ph()`: Schoenfeld residual proportional hazards tests for Cox models,
  including term-level and global tests.
- `compare_models()`: AIC, BIC, log-likelihood, likelihood-ratio tests, sample
  size, events for survival models, and optional primary-exposure tracking for
  gtregression candidate models.
- `select_models()`: `$results_table`, `$best_model`, `$all_models`, and
  `$direction`; `$table` is added when `format = gt` or `format = flextable`.