--- title: "Linked Cohort Analysis: Linking Case Records to Vaccination History" subtitle: "A complete starling workflow from preflight to perch" author: "Dr Nicolas Smoll, SCPHU, Sunshine Coast Hospital and Health Service" date: "`r Sys.Date()`" output: html_document: toc: true toc_depth: 3 toc_float: true theme: flatly pdf_document: toc: true toc_depth: 3 vignette: > %\VignetteIndexEntry{Linked Cohort Analysis} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE, eval = FALSE, # chunks are illustrative: reclin2 not available at build time fig.width = 8, fig.height = 5 ) ``` ## Overview This vignette demonstrates the complete `starling` probabilistic record linkage workflow: from pre-linkage data quality assessment through blocking variable construction, threshold sensitivity analysis, linkage, and post-linkage validation. The scenario mirrors a routine SCPHU task: linking a notifiable disease linelist (EDIS extracts) to the Australian Immunisation Register (AIR) to determine vaccination status at the time of disease onset. The datasets used (`cases_notifiable` and `vax_air`) are synthetic — no real person data. They include deliberate data quality issues (name typos, corrupted Medicare numbers) to demonstrate how the starling toolkit handles real-world messiness. ```{r load} library(starling) data(cases_notifiable) data(vax_air) cat("Cases linelist: ", nrow(cases_notifiable), "records\n") cat("Vaccination register:", nrow(vax_air), "records\n") cat("True matches (known):", sum(!is.na(cases_notifiable$true_link_id)), "\n") ``` --- ## Step 1: Pre-linkage audit with `preflight()` Before generating a single candidate pair, `preflight()` runs a structured battery of checks across both datasets: completeness of linkage variables, duplicate identifiers, date plausibility, Medicare validity, name field quality, and factor-level consistency. ```{r preflight} audit <- preflight( data1 = cases_notifiable, data2 = vax_air, linkage_vars = c("lettername1", "lettername2", "dob", "medicare10"), id_col1 = "id_var", id_col2 = "id_var", date_cols = c("dob", "onset_date"), medicare_col = "medicare10" ) ``` The audit flags include: - Any linkage variables with missingness above 10% - Duplicate ID values in either dataset - Medicare numbers that fail the Modulus 10 checksum - Date values before 1900 or after today --- ## Step 2: Medicare checksum validation with `check_medicare()` The `preflight()` report includes Medicare validity, but `check_medicare()` can also be called standalone for a more detailed summary and to add the validation flag column for downstream use. ```{r medicare} # Validate cases cases_checked <- check_medicare(cases_notifiable, medicare_col = "medicare10", output_col = "medicare_valid", verbose = TRUE) # Confirm AIR Medicare numbers are all valid vax_checked <- check_medicare(vax_air, medicare_col = "medicare10", output_col = "medicare_valid", verbose = TRUE) # Replace corrupted Medicare numbers with NA before linkage # so they don't negatively score a true match cases_checked$medicare10 <- ifelse( cases_checked$medicare_valid == 1L, cases_checked$medicare10, NA_character_ ) ``` The cases dataset has ~10% corrupted Medicare numbers by design. Setting those to `NA` before linkage is better than passing an invalid number, because the EM algorithm treats `NA` as "not observed" (no contribution to the score, positive or negative), whereas an invalid number that happens to match a wrong AIR record would add spurious positive weight. --- ## Step 3: Blocking variable construction with `flock()` `flock()` creates blocking keys that partition both datasets into candidate comparison groups. `murmuration()` only compares pairs within the same block, making the search tractable for large datasets. ```{r flock} # Extract birth year for composite blocking cases_blocked <- flock(cases_checked, block1_vars = "gender", block2_vars = "gender", block3_vars = "postcode", birth_year_col = "dob") vax_blocked <- flock(vax_checked, block1_vars = "gender", block2_vars = "gender", block3_vars = "postcode", birth_year_col = "dob") # Summary of blocking key distributions cat("block1 (gender) — unique values in cases:", dplyr::n_distinct(cases_blocked$block1), "\n") cat("block3 (postcode) — unique values in cases:", dplyr::n_distinct(cases_blocked$block3), "\n") ``` For this small synthetic dataset we use `block1` (gender only). For large production datasets (> 100 000 records), use multi-pass blocking: run `murmuration()` separately with `block1` and `block3`, then union the results. --- ## Step 4: Threshold sensitivity analysis with `perch()` Before committing to a threshold, we can use `perch()` standalone to understand the score landscape. Alternatively, `murmuration(perch_before_linking = TRUE)` calls `perch()` automatically mid-linkage after the EM model fits. ```{r perch-standalone, eval = FALSE} # This would require running the EM model first — # see the murmuration() call below which does this in one step. # For standalone use on a pre-scored pairs object: # pairs <- reclin2::pair_blocking(cases_blocked, vax_blocked, "block1") # reclin2::compare_pairs(pairs, # on = c("lettername1", "lettername2", "dob", "medicare10"), # default_comparator = reclin2::jaro_winkler(0.9), inplace = TRUE) # m <- reclin2::problink_em( # ~ lettername1 + lettername2 + dob + medicare10, data = pairs) # pairs_pred <- predict(m, pairs = pairs, add = TRUE) # # perch(pairs_pred, n_records_df1 = nrow(cases_blocked), # thresholds = seq(8, 25, by = 1)) ``` The threshold guidance from Australian linkage authorities: | Range | Source | Meaning | |---|---|---| | 10–20 | AIHW / WA Data Linkage Unit | Clerical review zone | | 15–20 | PHRN | Operational target for <0.5% false-match rate | | 17 | starling default | Balanced for routine surveillance | --- ## Step 5: Probabilistic linkage with `murmuration()` `murmuration()` runs the complete Fellegi-Sunter EM linkage pipeline in one call. We use `perch_before_linking = TRUE` to inspect the score distribution before the threshold is applied. ```{r murmuration} linked <- murmuration( df1 = cases_blocked, df2 = vax_blocked, linkage_type = "v2c", event_date = "onset_date", id_var = "id_var", blocking_var = "block1", compare_vars = c("lettername1", "lettername2", "dob", "medicare10"), threshold_value = 17, perch_before_linking = FALSE, # set TRUE in interactive sessions to inspect days_allowed_before_event = 14, clean_eggs = TRUE ) cat("Linked records: ", nrow(linked), "\n") cat("With vaccination: ", sum(!is.na(linked$vax_date_1)), "\n") cat("Without vaccination: ", sum( is.na(linked$vax_date_1)), "\n") ``` --- ## Step 6: Visualise the score distribution with `murmuration_plot()` Even if `perch_before_linking = FALSE` during the linkage call, we can still inspect the weight distribution afterwards by accessing the `weights` column on the linked output. ```{r plot, fig.cap = "Linkage weight distribution. The threshold (dashed line) should sit in the valley between the two score clusters."} # The linked output retains the weights column when clean_eggs = TRUE # For the visualisation, we use the weights from the linked data frame if ("weights" %in% names(linked)) { murmuration_plot(linked, threshold = 17, show_density = FALSE, palette = "sch") } ``` --- ## Step 7: Post-linkage validation Because `cases_notifiable` contains `true_link_id` (the ground-truth match identifier), we can compute recall and precision on the synthetic data. This step is only possible with synthetic data — in production, post-linkage validation requires a clerical review sample. ```{r validate} # Recall: proportion of true matches recovered true_positives <- sum( !is.na(linked$true_link_id) & !is.na(linked$id_var_df2) & linked$true_link_id == linked$id_var_df2, na.rm = TRUE ) total_true_matches <- sum(!is.na(cases_notifiable$true_link_id)) recall <- true_positives / total_true_matches # Precision: proportion of accepted links that are true matches total_links <- sum(!is.na(linked$id_var_df2)) precision <- true_positives / total_links cat(sprintf("Recall: %.1f%% (%d / %d true matches recovered)\n", recall * 100, true_positives, total_true_matches)) cat(sprintf("Precision: %.1f%% (%d / %d links are true matches)\n", precision * 100, true_positives, total_links)) cat(sprintf("F1 score: %.3f\n", 2 * precision * recall / (precision + recall))) ``` --- ## Summary: the complete starling workflow ```{r workflow-summary, eval = FALSE} library(starling) data(cases_notifiable); data(vax_air) # 1. Pre-linkage audit preflight(cases_notifiable, vax_air, linkage_vars = c("lettername1", "lettername2", "dob", "medicare10"), medicare_col = "medicare10") # 2. Medicare validation — replace invalid numbers with NA cases <- check_medicare(cases_notifiable) cases$medicare10 <- ifelse(cases$medicare_valid == 1L, cases$medicare10, NA_character_) # 3. Blocking variables cases <- flock(cases, block1_vars = "gender", birth_year_col = "dob") vax <- flock(vax_air, block1_vars = "gender", birth_year_col = "dob") # 4. Link (perch_before_linking = TRUE in interactive sessions) linked <- murmuration(cases, vax, linkage_type = "v2c", event_date = "onset_date", id_var = "id_var", blocking_var = "block1", compare_vars = c("lettername1", "lettername2", "dob", "medicare10"), threshold_value = 17) # 5. Pass to mudnester or bowerbird for downstream analysis ``` --- ## Session information ```{r session, eval = TRUE} sessionInfo() ```