---
title: "Diagnostics, classification and uncertainty"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Diagnostics, classification and uncertainty}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: ../inst/REFERENCES.bib
csl: apa.csl
link-citations: true
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = FALSE, comment = "",
                      fig.width = 7, fig.height = 4.5, dpi = 96,
                      dev.args = list(bg = "transparent"))
# Console colour carries no meaning on a rendered page. pkgdown turns it on for
# its own build, and the escape sequences then reach the reader as literal text,
# so colour is switched off here for a plain vignette render and a site build
# alike. The fixed width keeps printed output inside the documentation column.
options(cli.num_colors = 1, cli.hyperlink = FALSE, crayon.enabled = FALSE,
        width = 80)

# Figures on the package website sit on a warm off-white page in light mode and
# are inverted by pkgdown in dark mode, so an opaque background would read as a
# pale slab one way and a black plate the other. Two things paint one. The
# device canvas is made transparent by `dev.args` above, and theme_depictr()
# then inherits theme_minimal()'s white plot.background, which is drawn over
# that canvas, so it is cleared as each figure is printed. This is deliberately
# a vignette-level choice: theme_depictr() keeps its opaque background, which is
# what a figure saved for a paper wants.
transparent_bg <- ggplot2::theme(
  plot.background  = ggplot2::element_rect(fill = NA, colour = NA),
  panel.background = ggplot2::element_rect(fill = NA, colour = NA)
)
knit_print.ggplot <- function(x, ...) knitr::normal_print(x + transparent_bg)
knit_print.patchwork <- function(x, ...) knitr::normal_print(x & transparent_bg)

library(depictr)
has_simr <- requireNamespace("simr", quietly = TRUE)
```

## Regression diagnostics

`residual_diagnostics_plot()` assembles the classic panel. `influence_plot()`
and `qq_plot()` zoom in on particular checks. We fit a linear model to the
crop-yield trial.

```{r}
fit <- lm(yield ~ rainfall + fertiliser + soil_ph, data = crop_yield)
```

```{r, fig.height = 6}
residual_diagnostics_plot(fit, title = "Crop-yield model")
```

```{r}
influence_plot(fit)
```

The bubble area is Cook's distance, the standard measure of an observation's
influence on the fitted coefficients [@cook1977].

`vif_plot()` checks for multicollinearity among the predictors. The crop-yield
predictors are close to orthogonal, so to see the diagnostic do its job we add a
soil-moisture measure that is largely driven by rainfall.

```{r, fig.height = 3.2}
set.seed(1)
collinear <- crop_yield
collinear$soil_moisture <- 0.05 * collinear$rainfall +
  rnorm(nrow(collinear), sd = 2)
vif_plot(lm(yield ~ rainfall + soil_moisture + fertiliser, data = collinear))
```

Rainfall and soil moisture each carry a variance inflation factor near 6, above
the rule-of-thumb line at 5, while fertiliser stays at 1.

## Diagnostics for a generalised linear model

For a binary `glm` the raw residuals take only a few values, so a plain
residual-versus-fitted plot is hard to read. `binned_residual_plot()` instead
splits the data into equal-count bins of fitted values and plots the mean
residual per bin against a +/- 2 standard-error band [@gelman2007]. Most points
should sit inside the band, scattered around zero.

We model the *rare* adverse-event outcome from the clinical trial (about a 10%
base rate) on the baseline biomarker, age and treatment arm.

```{r}
gfit <- glm(adverse_event ~ biomarker + age + arm,
            data = clinical_trial, family = binomial)
```

```{r}
binned_residual_plot(gfit, title = "Binned residuals: adverse-event model")
```

`residual_diagnostics_plot()` recognises a `glm` and switches to
GLM-appropriate panels (binned residuals and a Q-Q plot of randomised quantile
residuals, which are standard-normal under a correct model regardless of the
response distribution):

```{r, fig.height = 6}
residual_diagnostics_plot(gfit, title = "Adverse-event model")
```

## Classification on an imbalanced outcome

The adverse-event outcome is imbalanced, the case the precision-recall, gains
and lift charts are built for. depictr's classification plots each read a
binomial `glm` directly, a pair of vectors or, to compare
several models, a named list of models.

The ROC curve reports the AUC. Passing `youden = TRUE` marks the Youden's J
operating point (the threshold maximising sensitivity + specificity - 1):

```{r, fig.width = 5, fig.height = 5}
roc_curve_plot(gfit, youden = TRUE)
```

When the positive class is rare the precision-recall curve is more informative
than the ROC curve, because it ignores the many true negatives. The baseline is
the positive prevalence, and `f1 = TRUE` marks the maximum-F1 operating point:

```{r, fig.width = 5, fig.height = 5}
pr_curve_plot(gfit, f1 = TRUE)
```

### Comparing two models in one figure

A *named list* overlays one colour-coded curve per model, with a per-curve AUC
(or average precision) in the legend. Here the full model is compared with a
biomarker-only model.

```{r}
reduced <- glm(adverse_event ~ biomarker, data = clinical_trial,
               family = binomial)
models <- list(Full = gfit, `Biomarker only` = reduced)
```

Because an ROC curve hugs the top-left, the bottom-right corner is always free,
so `legend_inside = TRUE` tucks the per-model legend there and reclaims the
right-hand margin:

```{r, fig.width = 5, fig.height = 5}
roc_curve_plot(models, youden = TRUE, legend_inside = TRUE)
```

```{r, fig.width = 5, fig.height = 5}
pr_curve_plot(models, f1 = TRUE)
```

For ranking and targeting tasks, the cumulative gains and lift charts show how
many positive cases are captured as more of the score-ordered population is
targeted, again overlaying both models:

```{r, fig.width = 5, fig.height = 5}
gain_plot(models, legend_inside = TRUE)
```

```{r, fig.height = 3.6}
lift_plot(models, legend_inside = TRUE)
```

### Choosing an operating point

`threshold_plot()` sweeps the decision threshold across the full range of scores
and plots each metric as it trades off, marking the Youden and maximum-F1
optimal thresholds with dashed lines. It is the natural companion to the ROC and
PR curves for actually *choosing* a cut-off.

```{r}
tp <- threshold_plot(gfit, title = "Metrics across the decision threshold")
tp
attr(tp, "thresholds")   # the Youden and max-F1 thresholds
```

### Calibration

`calibration_plot()` checks whether predicted probabilities match observed
frequencies. Each bin's observed rate carries a Wilson binomial confidence
interval, so bins backed by few observations (common in the sparse upper tail
when the outcome is rare) are not over-interpreted.

```{r, fig.width = 5, fig.height = 5}
calibration_plot(gfit, bins = 6)
```

A confusion-matrix heatmap completes the picture. Passing `threshold = "youden"`
reuses the same operating point the ROC curve marks, so the two agree:

```{r, fig.width = 5, fig.height = 4.5}
confusion_matrix_plot(gfit, threshold = "youden", normalise = "row")
```

## Uncertainty

`posterior_plot()` summarises draws (posterior, bootstrap, simulation) as a
distribution per parameter. Here are the posterior draws from a Bayesian fit of
the lexical-decision model, shown as point-and-interval forests.

```{r}
draws <- readRDS(
  system.file("extdata", "lexdec_draws.rds", package = "depictr")
)
posterior_plot(draws[c("conditionunrelated", "modalityauditory",
                       "word_frequency")],
               labels = c(conditionunrelated = "condition",
                          modalityauditory = "modality",
                          word_frequency = "word frequency"),
               style = "interval", title = "Posterior estimates (ms)")
```

## Power curves

`power_curve_plot()` reads a `simr::powerCurve()` object or a tidy data frame, so
a slow power simulation does not have to be re-run to redraw it. The package
ships a `powerCurve` object produced by a `simr` analysis of the lexical-decision
design, which simulated power for the word-frequency slope as the number of
participants grows. We read it straight from disk.

```{r, eval = has_simr}
pc <- readRDS(
  system.file("extdata", "powercurve_lexdec.rds", package = "depictr")
)
power_curve_plot(pc, x_lab = "Number of participants",
                 title = "Power for the word-frequency effect")
```

```{r, eval = !has_simr, echo = !has_simr}
# Summarising a powerCurve object needs the 'simr' package. Without it we read
# the same five points from the summary stored alongside the object, as a tidy
# data frame, which power_curve_plot() also accepts.
pc_df <- readRDS(system.file("extdata", "powercurve_lexdec_summary.rds",
                             package = "depictr"))
```

```{r, eval = !has_simr, echo = !has_simr, fig.height = 4.5}
power_curve_plot(pc_df, x_lab = "Number of participants",
                 title = "Power for the word-frequency effect")
```

The analysis ran 100 simulations at each of 12, 24, 36, 48 and 60 participants.
Power for the word-frequency slope is 88% with 12 participants and 99% with 24,
and it stays there for the larger samples. The confidence band is widest at the
smallest sample size, where its lower limit falls on the 80% target line, so 24
participants is the first design the simulation clears comfortably.

## Composing and saving

Combine any of these with `arrange_plots()` and save with `save_plot()`:

```{r, fig.height = 3.2}
panel <- arrange_plots(
  qq_plot(fit), influence_plot(fit),
  ncol = 2, title = "Diagnostics", tag_levels = "A"
)
panel
```

`save_plot()` writes it out at a print-ready 300 dpi by default, creating any
missing directories along the way.

```{r, eval = FALSE}
save_plot("figures/diagnostics.png", panel, width = 7, height = 3.2, dpi = 300)
```

## References