---
title: "Validating abbreviated symptom definitions"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Validating abbreviated symptom definitions}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width  = 7,
  fig.height = 4
)
```

A symptom definition is only useful if it holds up beyond the sample it was derived in. This vignette tests that, first by internal validation on a single dataset and then by external validation in an independent cohort.

## Why validation matters

The optimization fits the data it is given, so a definition that reproduces the full diagnosis almost perfectly in one sample may owe part of that fit to noise specific to that sample. Two questions separate a real result from an artefact of fitting. Does the same rule perform comparably on cases that were not used to derive it? And does it transport to a cohort with a different case mix? The first is answered by internal validation, holding out part of the sample or rotating through folds. The second is answered by external validation, deriving the rule in one cohort and applying it unchanged in another. Predictive values make the second question pressing: because they depend on how common PTSD is, a rule moved from a high-prevalence clinic to a low-prevalence community sample will behave differently even if its sensitivity and specificity are stable.

## Requirements for the input data

The input must be the 20 PCL-5 items in their standard order, scored 0 to 4, with no missing values, plus any identifier columns you name in `id_col`. The full contract is described in the [Getting started](getting-started.html) vignette.

## Internal validation

`holdout_validation()` partitions the sample into a training set and a test set, derives the top combinations on the training set, and reports their performance on the held-out test set. We optimize by balanced accuracy, the package default, which weighs performance in the diagnosed and the non-diagnosed group equally; the [Getting started](getting-started.html) vignette explains this choice. To keep the vignette fast we work with a 120-row subset of the bundled data.

```{r holdout, message = FALSE}
library(PTSDdiag)
library(dplyr)

data("simulated_ptsd")
ptsd <- rename_ptsd_columns(simulated_ptsd[1:120, ],
                            id_col = c("patient_id", "age", "sex"))

ho <- holdout_validation(
  ptsd,
  train_ratio = 0.7,
  n_symptoms  = 6,
  n_required  = 4,
  n_top       = 3,
  score_by    = "balanced_accuracy",
  seed        = 42
)
ho$without_clusters$best_combinations
ho$without_clusters$summary
```

`score_by = "balanced_accuracy"` maximizes the mean of sensitivity and specificity; `score_by = "sensitivity"` remains the conservative alternative when missing a true case is the costlier error, and `score_by = "accuracy"` minimizes total misclassification. The `seed` argument makes the split reproducible.

`cross_validation()` extends the same logic to k folds, deriving on k minus one folds and testing on the remaining one, then pooling the results. Combinations that recur across folds are reported in `combinations_summary`.

```{r cv, message = FALSE}
cv <- cross_validation(
  ptsd,
  k          = 2,
  n_symptoms = 6,
  n_required = 4,
  n_top      = 3,
  score_by   = "balanced_accuracy",
  seed       = 42
)
cv$without_clusters$summary_by_fold
cv$without_clusters$combinations_summary  # NULL if no combination repeats
```

## External validation in a second cohort

To test the generalizability of definitions based on fewer symptoms, one needs to evaluate the derived definitions in a second, independent dataset. The package ships a second simulated dataset, `simulated_ptsd_genpop`, whose PTSD prevalence is about 21%, well below the 94% of the included clinical sample. Deriving in the clinical sample and validating in the community sample therefore probes exactly the prevalence shift described above. We again use a 120-row subset of each for speed.

The rule we transport is the one the internal validation above already produced: the top combinations the holdout derivation selected on its training data. They are written to a small JSON file once, and read back before being applied. The export step is what makes the rule portable across sites and analysts, without a need for data sharing or manual transcription.

```{r external, message = FALSE}
data("simulated_ptsd_genpop")

# Export the holdout-derived combinations for reuse
tmp <- tempfile(fileext = ".json")
write_combinations(ho$without_clusters$best_combinations, tmp,
                   n_required = 4,
                   score_by = "balanced_accuracy",
                   description = "Six-symptom, four-required definition")

# A second analyst reads the file and applies it to the community sample.
# simulated_ptsd_genpop also carries paired CAPS-5 columns (C1..C20); here we
# use only the PCL-5 items, so we select those before standardising.
spec   <- read_combinations(tmp)
genpop <- rename_ptsd_columns(
  simulated_ptsd_genpop[1:120, c("patient_id", "age", "sex", paste0("S", 1:20))],
  id_col = c("patient_id", "age", "sex")
)

applied <- apply_symptom_combinations(genpop, spec$combinations,
                                      n_required = spec$n_required)
summarize_ptsd_changes(applied) %>%
  create_readable_summary()
```

Compare this table with the holdout test performance above, keeping in mind that the clinical test split contains only a handful of non-cases at 94% prevalence, so its specificity and predictive values are coarse. The community sample reveals the shifts external validation exists to expose. NPV rises sharply: non-cases dominate this sample and the rules miss few of them. Sensitivity is lower, because the community sample's symptom profiles are milder than the clinic's, while specificity is high. PPV, in contrast, hardly moves here — the transported rules produce almost no false positives in this sample, so a positive result stays trustworthy even at one-fifth the prevalence; with a less specific rule, the same prevalence drop would pull PPV down instead. None of these shifts is a failure of the rule; all are properties of applying a fixed criterion across settings that differ in prevalence and severity, and revealing them is exactly what external validation is for.

## How special is the winning subset?

Internal and external validation both look at the top of the ranking. A complementary question is how the rest of the candidate set behaves: if thousands of subsets perform nearly as well as the winner, the specific winning items should not be over-interpreted, because many symptom sets are effectively interchangeable. `score_all_combinations()` answers this by scoring every candidate combination — here all $\binom{20}{4} = 4{,}845$ four-symptom subsets — and returning the complete ranked table, the exhaustive companion to `optimize_combinations()`. Plotting the ranking metric against rank typically shows a plateau of near-optimal subsets before performance falls away; the width of that plateau is the interchangeability of the solution.

```{r all-combinations, fig.alt = "Balanced accuracy of every four-symptom combination by rank"}
curve <- score_all_combinations(ptsd, n_symptoms = 4, n_required = 3,
                                show_progress = FALSE)
nrow(curve)
head(curve, 3)

plot(curve$rank, curve$balanced_accuracy, type = "l", log = "x",
     xlab = "Combination rank (log scale)", ylab = "Balanced accuracy")
```

`compute_plateau()` turns that visual impression into a number. It keeps every combination within `delta` percentage points of the best and reports the count both absolutely and as a share of the candidate space that was actually searched. Report the share, not the count: the hierarchical rule searches a much smaller space than the non-hierarchical one, so identical counts describe different degrees of interchangeability. The primary width used in the preprint is one percentage point, with five as a secondary, and both come from one call.

```{r plateau}
plateau <- compute_plateau(curve, delta = c(1, 5))
plateau$summary
```

A plateau measured in one sample is still a single-sample result. `bootstrap_stability()` resamples the respondents, re-runs the entire search in each replicate, and asks how often each combination stays within `delta` of that replicate's winner. The reference point is the winner of the replicate, not of the original data, so every comparison stays on the scale of its own resample. Fifty replicates keep this vignette fast; a reported analysis would use `n_boot = 2000`.

```{r stability}
stability <- bootstrap_stability(curve, ptsd, delta = 1, n_boot = 50,
                                 show_progress = FALSE)
stability$summary

# How often does the sample's winner survive resampling?
head(stability$per_set, 3)
```

`pi` is the proportion of replicates in which a combination remained near-optimal, and `p_argmax` the proportion in which it was the best. A winner with a low `p_argmax` is a winner by accident of sampling. Read these with the selection in mind: for a combination fixed in advance they estimate how often it is near-optimal in samples like this one, but for the combination that won in this particular sample they are optimistic, because bootstrap replicates share most of their observations with the sample that singled it out. The plateau size does not depend on which combination was selected and is unaffected.

Finally, which items are doing the work? Counting how often each item appears across the plateau is not enough on its own, because every combination holds the same number of items and the counts are therefore zero-sum: one item can only gain at another's expense. `symptom_selection()` reports each item's frequency alongside the frequency expected from drawing combinations at random, and their ratio.

```{r selection, fig.alt = "Selection frequency of each PCL-5 item across the plateau, with the chance baseline marked"}
selection <- symptom_selection(curve, delta = 1, stability = stability)
head(selection[order(-selection$enrichment),
               c("symptom", "cluster", "baseline", "plateau_freq", "enrichment")], 5)

plot_symptom_selection(selection)
```

Items with an `enrichment` above 1 are selected more often than chance would produce. The baseline is drawn item by item rather than as a single line because it is not constant under the hierarchical rule: there a combination must take at least one item from every cluster, so each of the two avoidance items appears in 55.8% of the candidate space while each of the seven cognition-and-mood items appears in only 25.2%. Judging hierarchical results against a flat baseline would make the small clusters look preferred when they are simply unavoidable.

To place several rule forms side by side, `compare_rule_forms()` takes a named list of searches and stacks the tables above with a `rule_form` column.

## See also

- [Getting started](getting-started.html) for the single-cohort derivation workflow.
- [Comparing diagnostic criteria](comparing-criteria.html) for the symptom-frequency heatmap and multi-rule comparison.
- [CAPS-5 workflow](caps5-workflow.html) for validation against a clinician-administered reference.