Advanced Capstone: Compare Families, Not Just Scores

You already completed the main binary-model capstone in Chapter 18. This extended capstone is the last chapter of Part VI and of the course. It draws on the method chapters of Parts III to V, and on the selection rules of chapters 29 and 30, to ask you to demonstrate methods literacy across supervised and unsupervised tasks. All datasets are synthetic or built into R; never import real patient information.

1. Run the method labs

Start from the project folder in RStudio:

source("scripts/01_install_packages.R")
source("scripts/03_smoke_test.R")
source("scripts/05_supervised_comparison.R")
source("scripts/06_unsupervised_comparison.R")
source("scripts/07_multiclass_demo.R")

The scripts create output files in an outputs/ folder inside your project. They should be reproducible with the same R/package versions and seed.

Example: Before sourcing anything, check that the lab scripts are where the project root expects them. file.exists() gives TRUE or FALSE for each name.

labs <- c("scripts/01_install_packages.R", "scripts/03_smoke_test.R",
          "scripts/05_supervised_comparison.R", "scripts/06_unsupervised_comparison.R",
          "scripts/07_multiclass_demo.R")
file.exists(labs)

Output

[1] TRUE TRUE TRUE TRUE TRUE

All five scripts exist (TRUE). A FALSE means R is not running from the project root, so fix the working folder (RStudio: open the .Rproj file) before sourcing anything.

2. Project A: supervised numeric prediction

Question: Use synthetic EHR-like measurements to predict known fictional systolic BP (a continuous outcome).

Describe what each row means and which rows have an observed SBP target.
Make a train/test split before learning imputation statistics; keep test outcomes locked.
Compare the training mean, linear regression, a regression tree, and random forest using the same folds.
Explain RMSE and MAE units, whether scaling is necessary, and whether errors differ across sites.
Select a model on training CV only and record what independent validation is still missing.

Use scripts/05_supervised_comparison.R as the worked training-CV starting point.

Example: Steps 1 and 2 in code. Each row is one fictional patient; keep rows with an observed SBP and BMI (so every model is scored on the same rows), split first, and count how many rows land in each part.

library(tidymodels)
ehr <- readr::read_csv("data/synthetic_ehr.csv", show_col_types = FALSE)
a_data <- ehr |> filter(!is.na(sbp), !is.na(bmi)) |> select(sbp, age, bmi)
set.seed(31)
a_split <- initial_split(a_data, prop = 0.75)
a_train <- training(a_split)
a_test <- testing(a_split)
c(complete_rows = nrow(a_data), dropped_rows = nrow(ehr) - nrow(a_data),
  train = nrow(a_train), test = nrow(a_test))

Output

complete_rows  dropped_rows         train          test 
          546            54           409           137

546 of the 600 rows have both SBP and BMI; the other 54 are left out so that every model is scored on the same rows. The split gives 409 training rows and 137 test rows. The test rows stay locked: every comparison below uses a_train only.

3. Project B: supervised binary classification

Question: Predict the fully known synthetic diabetes outcome (not the real-world proxy “diagnosis recorded”).

Report positive support and the majority-class baseline.
Compare logistic regression, tree, random forest, and kNN with identical training folds.
Explain which preprocessor parameters are learned inside each fold.
Report ROC AUC / PR AUC, not just accuracy.
Explain how you would select a threshold and assess calibration before final testing.

Do not pick whichever algorithm has the highest test score from repeated tries; this demo intentionally focuses on development-stage comparison.

Example: Step 1: positive support and the accuracy of always answering “No”.

table(ehr$diabetes)

Output

 No Yes 
546  54

round(mean(ehr$diabetes == "No"), 3)

Output

[1] 0.91

The file has 54 diabetes cases (9%) and 546 non-cases, so always answering “No” is right 91% of the time. Any model has to be judged against that 0.91, which is why accuracy alone is not enough and ROC AUC and PR AUC are asked for.

4. Project C: multiclass classification

Use the built-in iris flower dataset to avoid falsely suggesting that the synthetic EHR has validated three-way labels.

Explain why this is multiclass, not multilabel.
Compare fold-wise LDA, QDA, kNN, decision tree, and multinomial logistic regression.
Interpret which categories are confused, not only overall accuracy.
Explain how you would modify evaluation for an ordinal or multilabel outcome.

Example: Each iris flower has exactly one of three species, which makes it multiclass. Here is leave-one-out LDA, so every flower is predicted by a model that never saw it, shown as a confusion table.

table(iris$Species)

Output

    setosa versicolor  virginica 
        50         50         50

lda_pred <- MASS::lda(Species ~ ., data = iris, CV = TRUE)$class
table(predicted = lda_pred, actual = iris$Species)

Output

            actual
predicted    setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         48         1
  virginica       0          2        49

Setosa is predicted perfectly. All three mistakes are between versicolor and virginica: two versicolor flowers were called virginica and one virginica was called versicolor. That is the confusion to describe, not just an overall accuracy.

5. Project D: unsupervised methods

Use only the fake predictor columns age, BMI, SBP, pulse.

Impute and scale the columns (descriptive exercise).
Compare K-means, PAM, and hierarchical clustering. Describe the shapes/distances each method assumes.
Choose k for K-means with both the elbow plot (tot.withinss) and the average silhouette width for k = 2 to 6, and report both. If they disagree or are weak, say so and explain which k you use and why. Set a seed and use nstart of 10 or more.
Run PCA. Report the share of variance for the first two components, read the loadings in words (which variables drive PC1 and PC2). Optionally, draw a scree plot and the clusters on the first two components.
Run arules::apriori() on binned columns (use cut()) with stated supp, conf and maxlen. Report the number of rules and the top rules by lift with their count, and add a caution: lift near 1 means nothing beyond the base rate, and a high lift on a small count may be chance. Also compute one rule’s support/confidence/lift by hand to check the inspect() output.
Identify potential extreme values with a robust score. Do not call clusters “high diabetes risk” or anomalies “disease.”

Example: Step 1 in code. Impute missing values with the median, standardize, and check the result has mean 0 and standard deviation 1 in every column. This matrix x is reused in the Try it below.

x_raw <- ehr |> select(age, bmi, sbp, pulse)
prep_plan <- recipe(~ age + bmi + sbp + pulse, data = x_raw) |>
  step_impute_median(all_numeric_predictors()) |>
  step_normalize(all_numeric_predictors())
x <- as.matrix(bake(prep(prep_plan, training = x_raw), new_data = x_raw))
dim(x)

Output

[1] 600   4

round(colMeans(x), 3)

Output

  age   bmi   sbp pulse 
    0     0     0     0

round(apply(x, 2, sd), 3)

Output

  age   bmi   sbp pulse 
    1     1     1     1

All 600 rows are kept (missing values were filled with the median) and each column now has mean 0 and standard deviation 1, so no column dominates a distance just because of its units.

6. Choose one specialized extension

Write a one-page plan (no large software build required) for one of:

Right-censored return-to-clinic survival analysis.
Temporal forecasting of fictional aggregate visits.
Semi-supervised / positive–unlabeled learning when diagnoses are incomplete.
Anomaly detection for fictional data-quality audits.
Multilabel classification where several fake codes can be present.
Self-supervised representation learning on approved synthetic text.

Describe the observation unit, target (if any), required labels or assumptions, split, metrics, baseline, limitations, and who must review the work before clinical use.

Example (scenario): A learner picks the first option. Unit: one fictional patient after a hospital stay. Target: days until the first clinic visit, with patients who never returned by day 90 marked censored (event = 0, time = 90). Split: by patient, with the last three months of discharges as the test period. Baseline: a Kaplan–Meier curve with no predictors. Metric: concordance and a calibration check at 30 days. Limitation: patients who moved away look censored. Reviewer: a clinician and a data steward, before anyone uses it. Each part of the plan is a concrete sentence a reviewer can accept or challenge.

7. Submission format and rubric

A complete submission includes:

A 1–2 page plan with one main question, data dictionary, allowable predictors, and validation strategy.
Reproducible R code for Projects A–D (can reuse and annotate provided scripts).
Two plots, a numeric evaluation table with denominators, and one paragraph interpreting mistakes.
A risk memo addressing missing labels, patient/time leakage, uncertainty, and ethics.

Criterion	Points
Correct choice of task, paradigm, and label definition	20
Reproducible imports, audit, and leak-safe preprocessing	20
Fair supervised comparison across shared training folds	20
Accurate unsupervised interpretation and limitations, including k chosen with elbow and silhouette, PCA scree and loadings read in words, and apriori thresholds and lift reported with a caution	15
Proper metrics, positive/class support, and uncertainty	15
Clear writing, safe synthetic-data handling, and peer review	10

Example: Check your own tally before submitting. A submission that earns 18, 20, 15, 12, 15 and 8 on the six rows:

earned <- c(18, 20, 15, 12, 15, 8)
possible <- c(20, 20, 20, 15, 15, 10)
c(total = sum(earned), out_of = sum(possible))

Output

 total out_of 
    88    100

The maximum is 100. This submission has 88, and the largest gap is on the shared-folds row (5 points lost), so that is where a rewrite should start.

Self-check before you submit

Before you submit, check: can you say yes to each question?

Do you know that classification is a type of supervised learning?
Does every modeling task have a correct reference label (or do you say there is none)?
Do you describe unsupervised clusters as similarity groups rather than diagnoses?
Does a new run from the project root reproduce your outputs without real patient data?
Are preprocessing and model selection confined to training folds?
Do you keep final test, calibration, threshold and uncertainty separate from training metrics?
Can you describe what would change for a new site, later year, or repeated patient visits?

8. What a successful completion means

You should be able to explain when to use regression, binary/multiclass classification, clustering, PCA, anomaly detection, association rules, forecasting, or survival and how each would be evaluated. You are not expected to train every existing algorithm or to deploy anything clinically.

Example (scenario): Three questions, three methods. “Will this patient’s SBP be above 140 next week?” is binary classification, judged by ROC AUC on held-out patients. “Which clinics see the same mix of patients?” is clustering, judged by silhouette width and by whether the groups are stable, never by a diagnosis. “How many visits will we have in week 30?” is forecasting, judged by MAE on later weeks. Matching the question to the method comes first; the algorithm choice comes second.

Try it

These five questions practise the capstone steps on the synthetic data. The solutions reuse objects built in the examples of sections 2 to 5.

Project A: compare the training mean, linear regression, a regression tree and a random forest on the same training folds. Which has the lowest cross-validated RMSE, and what does RMSE measure in this problem?

Show solution

Answer: Linear regression has the lowest cross-validated RMSE, 7.12 mmHg, against 8.28 for the tree, 8.58 for the forest and 13.0 for the mean baseline. RMSE is the typical size of a prediction error in the same unit as the outcome (mmHg), with large misses counted more heavily.

Step by step

Build five folds on the training rows only, once, and use them for every model.
Score each model with RMSE and MAE (both in mmHg). A small helper scores one model spec.

set.seed(31)
a_folds <- vfold_cv(a_train, v = 5)
cv_a <- function(spec) {
  workflow() |> add_formula(sbp ~ age + bmi) |> add_model(spec) |>
    fit_resamples(resamples = a_folds, metrics = metric_set(rmse, mae)) |>
    collect_metrics() |> select(.metric, mean)
}
bind_rows(
  mean_baseline = cv_a(null_model() |> set_mode("regression")),
  linear = cv_a(linear_reg()),
  tree = cv_a(decision_tree() |> set_mode("regression")),
  forest = cv_a(rand_forest(trees = 200) |> set_engine("ranger") |> set_mode("regression")),
  .id = "model"
) |> pivot_wider(names_from = .metric, values_from = mean)

Output

# A tibble: 4 x 3
  model           mae  rmse
  <chr>         <dbl> <dbl>
1 mean_baseline 10.6  13.0 
2 linear         6.41  7.12
3 tree           7.05  8.28
4 forest         7.54  8.58

3. Read the table: every model beats the mean baseline, and the simple linear model is best on both RMSE and MAE (RMSE 7.12, MAE 6.41). The tree and forest do not gain from the extra flexibility with only age and BMI. Choose on these training-CV numbers; the 137 test rows stay locked, and a new site or later year would still need its own check. (Rows with a missing bmi were removed before the split, so every model is scored on the same rows.)

Project B: report the positive support and majority baseline for diabetes, then compare logistic regression, a tree, a random forest and kNN on the same folds using ROC AUC and PR AUC.

Show solution

Answer: After dropping rows with a missing predictor there are 517 rows, 44 of them positive (8.5%), so always answering “No” is right about 0.915 of the time and a model with no skill has a PR AUC of about 0.085. Ranked by ROC AUC on identical folds, kNN is highest (0.795), then logistic regression (0.770), the forest (0.748) and the tree (0.674). PR AUC gives a slightly different order: kNN 0.386, logistic 0.377, tree 0.317, forest 0.307.

Step by step

Clean the data (drop rows with a missing predictor), count the positives and the majority baseline on those rows, then build one set of folds with the predictors standardized inside each fold by the recipe.

b_data <- ehr |> select(diabetes, age, bmi, sbp, pulse) |>
  mutate(diabetes = factor(diabetes, levels = c("Yes", "No"))) |> drop_na()
table(b_data$diabetes)

Output

Yes  No 
 44 473

round(c(majority_accuracy = mean(b_data$diabetes == "No"),
        no_skill_pr_auc = mean(b_data$diabetes == "Yes")), 3)

Output

majority_accuracy   no_skill_pr_auc 
            0.915             0.085

b_rec <- recipe(diabetes ~ ., data = b_data) |> step_normalize(all_numeric_predictors())
set.seed(31)
b_folds <- vfold_cv(b_data, v = 5, strata = diabetes)
cv_b <- function(spec) {
  workflow() |> add_recipe(b_rec) |> add_model(spec) |>
    fit_resamples(resamples = b_folds, metrics = metric_set(roc_auc, pr_auc)) |>
    collect_metrics() |> select(.metric, mean)
}
bind_rows(
  logistic = cv_b(logistic_reg()),
  tree = cv_b(decision_tree() |> set_mode("classification")),
  forest = cv_b(rand_forest(trees = 200) |> set_engine("ranger") |> set_mode("classification")),
  knn = cv_b(nearest_neighbor(neighbors = 15) |> set_engine("kknn") |> set_mode("classification")),
  .id = "model"
) |> pivot_wider(names_from = .metric, values_from = mean)

Output

# A tibble: 4 x 3
  model    pr_auc roc_auc
  <chr>     <dbl>   <dbl>
1 logistic  0.377   0.770
2 tree      0.317   0.674
3 forest    0.307   0.748
4 knn       0.386   0.795

2. Read the table: all four models beat the no-skill PR AUC of 0.085 and the 0.915 accuracy is the number any model must be judged against. The answer’s order is by ROC AUC (the printed table is in model-list order); by PR AUC the tree (0.317) is slightly ahead of the forest (0.307). The gap between kNN and logistic regression is small next to the noise of 5 folds on 44 positives, so it does not settle which is better; the logistic model is simpler and nearly as good. No test-set score was used to choose.

Project C: why are the iris species a multiclass problem and not multilabel, and which species does LDA confuse?

Show solution

Answer: Each flower has exactly one species, so the classes are mutually exclusive: multiclass, not multilabel. LDA (leave-one-out) is 98% accurate and its only confusions are between versicolor and virginica.

Step by step

A flower has exactly one species, so the three classes are exclusive (multiclass). In a multilabel problem a record could have several labels at once.
Read the confusion table from section 4 for the off-diagonal cells, and compute the overall accuracy for reference.

cm <- table(predicted = lda_pred, actual = iris$Species)
cm

Output

            actual
predicted    setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         48         1
  virginica       0          2        49

round(sum(diag(cm)) / sum(cm), 3)

Output

[1] 0.98

3. Read it: 147 of 150 flowers are right (0.98). The three errors are all in the versicolor/virginica block, and setosa is never confused. Those two species overlap in their measurements, which is the part to explain. 4. For an ordinal outcome (for example mild, moderate, severe), use a metric that counts a two-step miss worse than a one-step miss, such as weighted kappa or mean absolute class distance. For a multilabel outcome, score each label separately (one ROC AUC per label) and report the average.

Project D: choose k for K-means with the elbow (tot.withinss) and the average silhouette width for k = 2 to 6. Report both and say which k you would use.

Show solution

Answer: I would use k = 2, the only k with an average silhouette above 0.25 (0.276), and flag it as weak structure. Every value is low (0.23 to 0.28), and the within-group spread falls smoothly from 1637 to 860 with no sharp elbow, so the data show no clearly separated groups. Report both rows and justify the k you keep by its usefulness, not by a score.

Step by step

Run K-means for each k with a seed and nstart = 10, and record the total within-cluster sum of squares and the average silhouette width.

set.seed(31)
ks <- 2:6
d_x <- dist(x)
fits <- lapply(ks, function(k) kmeans(x, centers = k, nstart = 10))
tw <- sapply(fits, function(f) f$tot.withinss)
sw <- sapply(fits, function(f) mean(cluster::silhouette(f$cluster, d_x)[, "sil_width"]))
round(rbind(k = ks, tot_withinss = tw, avg_silhouette = sw), 3)

Output

                   [,1]     [,2]     [,3]    [,4]    [,5]
k                 2.000    3.000    4.000   5.000   6.000
tot_withinss   1637.170 1362.894 1140.348 969.630 859.774
avg_silhouette    0.276    0.232    0.235   0.246   0.247

2. Read the rows: tot_withinss always falls as k grows, so the elbow is a judgment call and there is no obvious bend here. avg_silhouette peaks at k = 2 and is only 0.232 to 0.247 for k = 3 to 6. On the Kaufman and Rousseeuw scale, an average of 0.25 or below means no substantial structure, and 0.26 to 0.50 means weak structure that may be artificial. So k = 2 (0.276) is weak at best and k = 3 to 6 (0.23 to 0.25) show none. That is why k = 2 is the pick, with the caveat written next to it. If you prefer another k for a description, say so and say why.

Project D: run PCA and report the share of variance for the first two components and the loadings in words. Then run apriori() with stated thresholds, report the top rules by lift with their count, check one rule by hand, and flag extreme SBP values with a robust score.

Show solution

Answer: PC1 explains 45.5% and PC2 25.4% of the variance (70.9% together). PC1 is driven by SBP and age together, PC2 mostly by BMI. With supp = 0.05, conf = 0.5, minlen = 2 and maxlen = 3, apriori finds 142 rules; the top by lift is {sbp=140+, family_history=No} => {age=60+} (lift 2.63, count 42). The hand check of {sbp=140+} => {diabetes=Yes} gives support 0.038, confidence 0.397 and lift 4.41, the same as inspect() (count 23). No SBP value is flagged by the robust score (largest 2.62).

Step by step

Run PCA on the standardized matrix and print the share of variance of every component and the loadings of the first two.

pc <- prcomp(x)
share <- pc$sdev^2 / sum(pc$sdev^2)
round(share, 3)

Output

[1] 0.455 0.254 0.246 0.045

round(pc$rotation[, 1:2], 2)

Output

        PC1   PC2
age   -0.64 -0.37
bmi   -0.28  0.88
sbp   -0.70  0.04
pulse -0.12 -0.31

2. Read the printed shares: the first two are 0.455 and 0.254, then 0.246 and 0.045, so the variance drops sharply only after the third component. PC1 has large loadings of the same sign for sbp (-0.70) and age (-0.64): it describes older people with higher SBP (the overall sign of a component is arbitrary). PC2 is mostly bmi (0.88). These describe how the columns vary together, not diseases. 3. Bin the columns with cut() and run apriori() with the thresholds written next to the table (supp = 0.05, conf = 0.5, minlen = 2, maxlen = 3).

library(arules)
bins <- data.frame(
  age = cut(ehr$age, c(0, 40, 60, Inf), labels = c("<40", "40-59", "60+"), right = FALSE),
  bmi = cut(ehr$bmi, c(0, 25, 30, Inf), labels = c("<25", "25-29", "30+"), right = FALSE),
  sbp = cut(ehr$sbp, c(0, 120, 140, Inf), labels = c("<120", "120-139", "140+"), right = FALSE),
  family_history = factor(ehr$family_history, levels = 0:1, labels = c("No", "Yes")),
  diabetes = factor(ehr$diabetes, levels = c("No", "Yes"))
)
tx <- as(bins, "transactions")
quiet <- list(verbose = FALSE)
rules <- apriori(tx, parameter = list(supp = 0.05, conf = 0.5, minlen = 2, maxlen = 3),
                 control = quiet)
length(rules)

Output

[1] 142

inspect(head(sort(rules, by = "lift"), 3))

Output

    lhs                              rhs       support    confidence coverage  
[1] {sbp=140+, family_history=No} => {age=60+} 0.07000000 0.9545455  0.07333333
[2] {sbp=140+}                    => {age=60+} 0.09166667 0.9482759  0.09666667
[3] {bmi=30+, sbp=140+}           => {age=60+} 0.06500000 0.9285714  0.07000000
    lift     count
[1] 2.627189 42   
[2] 2.609934 55   
[3] 2.555701 39

4. Read it: 142 rules pass the thresholds. The top three by lift all end in age=60+ with lift 2.56 to 2.63 and confidence about 0.93 to 0.95: older age is far more common among rows with SBP of 140 or more than in the file overall. Report the counts (42, 55 and 39 rows) beside the lift. Caution: lift near 1 means nothing beyond the base rate, and a high lift on a small count may be chance. This is an association, not a cause. 5. Check one rule by hand: {sbp=140+} => {diabetes=Yes}. Support is the share of rows with both items, confidence is that share divided by the share with sbp=140+, and lift is confidence divided by the share with diabetes=Yes.

n_rows <- nrow(bins)
both <- sum(bins$sbp == "140+" & bins$diabetes == "Yes", na.rm = TRUE) / n_rows
lhs <- sum(bins$sbp == "140+", na.rm = TRUE) / n_rows
rhs <- sum(bins$diabetes == "Yes", na.rm = TRUE) / n_rows
round(c(support = both, confidence = both / lhs, lift = both / lhs / rhs), 3)

Output

   support confidence       lift 
     0.038      0.397      4.406

inspect(apriori(tx, parameter = list(supp = 0.01, conf = 0.01, minlen = 2, maxlen = 2),
                appearance = list(lhs = "sbp=140+", rhs = "diabetes=Yes"), control = quiet))

Output

    lhs           rhs            support    confidence coverage   lift    count
[1] {sbp=140+} => {diabetes=Yes} 0.03833333 0.3965517  0.09666667 4.40613 23

6. Read it: support 0.038 is 23 rows out of 600, confidence 0.397 is 23 of the 58 rows with sbp=140+, and lift 4.41 is that confidence divided by the 9% of rows with diabetes. The inspect() row shows the same numbers (count 23). The hand calculation divides by all 600 rows, as apriori() does; dividing by only the 580 rows with an SBP value would give slightly different numbers. 7. A robust score for extreme SBP uses the median and the median absolute deviation (MAD): |sbp - median| / mad. Count values with a score above 3.5 (a common rule of thumb), and treat them as rows to check, not as disease.

robust_z <- abs(ehr$sbp - median(ehr$sbp, na.rm = TRUE)) / mad(ehr$sbp, na.rm = TRUE)
c(flagged = sum(robust_z > 3.5, na.rm = TRUE),
  largest_score = round(max(robust_z, na.rm = TRUE), 2))

Output

      flagged largest_score 
         0.00          2.62

8. Read it: flagged is 0 and the largest robust score is 2.62, below 3.5. The simulated SBP values contain no extreme points by this rule. If any had been flagged, they would be rows to check for entry errors, not patients with a disease.

Final checkpoint

Teach another volunteer this sentence, with examples: “We choose the learning task from the question, select the algorithm from the data and constraints, and judge the result using independent evidence rather than an attractive score.”

This is the end of the course. For reference material, see the glossary, the governance appendix, and the exercise solutions index pages in the appendices.