---
title: "Multivariate analysis and survival"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Multivariate analysis and survival}
  %\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)
```

Beyond regression, depictr covers three staples of applied data analysis:
principal component analysis, clustering (with quality diagnostics) and survival
curves.

## Principal component analysis

`pca_plot()` runs a PCA on the numeric columns of a data frame and draws a
biplot: the observations projected onto two components, with the variable
loadings as arrows. `scree_plot()` shows how much variance each component
explains.

```{r}
num <- c("rainfall", "fertiliser", "soil_ph", "yield")
pca_plot(crop_yield, cols = num, group = "treatment",
         title = "Crop-yield PCA")
```

```{r, fig.height = 4}
scree_plot(crop_yield, cols = num)
```

Both functions also accept a ready-made `prcomp()` object, so you can analyse
once and plot several views:

```{r}
pc <- prcomp(crop_yield[num], scale. = TRUE)
pca_plot(pc, components = c(1, 3))
```

## Clustering

`cluster_plot()` runs k-means and shows the clusters on the first two principal
components (so it works for any number of variables), with convex hulls and
labelled centroids. `dendrogram_plot()` draws a hierarchical-clustering tree and
can cut it into `k` groups.

```{r}
cluster_plot(crop_yield, cols = num, k = 3, seed = 1,
             title = "Crop-yield clusters")
```

```{r, fig.height = 4}
region_means <- aggregate(
  cbind(stress, sleep_hours, life_satisfaction, age, income) ~ region,
  data = wellbeing_survey, FUN = mean
)
rownames(region_means) <- region_means$region
dendrogram_plot(region_means[-1], k = 2, title = "Regions clustered")
```

## Choosing the number of clusters

Choosing `k` should not be guesswork. `k_diagnostic()` evaluates a
cluster-quality criterion across a range of `k` and suggests a value, using the
average silhouette width (the default, after Rousseeuw, 1987), the
within-sum-of-squares elbow or the gap statistic (Tibshirani, Walther & Hastie,
2001). See `?k_diagnostic` for the references.

```{r, fig.height = 3.5}
kd <- k_diagnostic(crop_yield, k_range = 2:6, cols = num,
                   method = "silhouette")
kd  # the diagnostic curve, with the suggested k marked
```

The suggested `k` and the underlying table are attached as attributes:

```{r}
attr(kd, "suggested")
knitr::kable(attr(kd, "k_table"), digits = 3)
```

The gap statistic compares the within-cluster dispersion against a null
reference, so unlike the other two criteria it can also support `k = 1`. Giving
it a range that starts there shows what it makes of these data. The reference
samples are drawn at random, so the chunk is seeded for a reproducible build.

```{r, fig.height = 3.5}
set.seed(1)
k_diagnostic(crop_yield, k_range = 1:6, cols = num, method = "gap")
```

The verdict is `k = 1`, which is worth taking seriously: the three-cluster
partition drawn above is a description imposed on the data rather than evidence
of groups within it.

`silhouette_plot()` then shows the quality of an actual clustering, one bar per
observation grouped by cluster, with the average silhouette width per cluster and
overall (the dashed line). Wide positive bars are well-placed observations, and a
bar below zero marks an observation that sits closer to a neighbouring cluster
than to its own.

```{r, fig.height = 5}
set.seed(1)
cl <- kmeans(scale(crop_yield[num]), centers = attr(kd, "suggested"),
             nstart = 10)$cluster
silhouette_plot(crop_yield, cl, cols = num,
                title = "Silhouette widths by cluster")
```

The overall average here is 0.24, and only six of the two hundred bars fall
below zero, none of them by much. Rousseeuw's rough guide treats an average at
or below 0.25 as no substantial structure, so this figure agrees with the gap
statistic: k-means returns clusters whether or not the data contains any, and
the silhouette is one way of asking which case you are in.

The same diagnostics apply to the wellbeing survey's numeric profile:

```{r, fig.height = 4}
wb_num <- c("age", "income", "stress", "sleep_hours",
            "exercise_days", "life_satisfaction")
wkd <- k_diagnostic(wellbeing_survey, k_range = 2:6, cols = wb_num,
                    method = "wss")
wcl <- kmeans(scale(na.omit(wellbeing_survey[wb_num])),
              centers = attr(wkd, "suggested"), nstart = 10)$cluster
silhouette_plot(na.omit(wellbeing_survey[wb_num]), wcl,
                title = sprintf("Wellbeing survey, k = %d",
                                attr(wkd, "suggested")))
```

## Survival curves

`survival_plot()` draws Kaplan-Meier curves [@kaplan1958] with stepwise
confidence limits from Greenwood's formula [@greenwood1926] and censoring marks.
The estimate is computed in base R, so no modelling package is needed. You can
pass follow-up times and an event indicator directly, a data frame or a
`survival::survfit` object.

The `clinical_trial` dataset has two arms with a real survival difference (the
treatment arm has a lower hazard), so the curves separate. Turning on the three
publication annotations gives a survminer-style figure: a number-at-risk
table beneath the curves, dashed guides to each arm's median survival and
the log-rank test of the difference as a subtitle. A survival curve is
monotone-decreasing, so its bottom-left corner is always empty, and
`legend_inside = TRUE` puts the group legend there.

```{r, fig.height = 6}
survival_plot(
  clinical_trial$time, clinical_trial$event, group = clinical_trial$arm,
  risk_table = TRUE, median_line = TRUE, logrank = TRUE, legend_inside = TRUE,
  x_lab = "Months", title = "Overall survival by arm"
)
```

The log-rank p-value is tiny and the median survival is clearly longer in the
treatment arm. Because its event times are longer, that arm is more often
right-censored at the 36-month study end, visible as the denser run of censoring
marks.

## References