Threshold Selection: Using perch() and murmuration_plot()
Understanding and calibrating the Fellegi-Sunter linkage cutoff
Dr Nicolas Smoll, SCPHU, Sunshine Coast Hospital and Health Service
2026-07-10
Why the threshold matters
In the Fellegi-Sunter probabilistic linkage framework, every candidate pair of records receives a composite log-likelihood ratio score. Pairs scoring above a chosen threshold are classified as matches; those below are classified as non-matches. This threshold is dataset-specific — there is no universal correct value, because the score scale depends entirely on which variables are compared, their quality, and the composition of the two datasets being linked.
Choosing the threshold too low increases recall (more true matches recovered) at the cost of precision (more false matches accepted). Choosing it too high does the reverse. The optimal value sits in the valley between the bimodal score distribution — the gap between the non-match cluster at low scores and the match cluster at high scores.
starling provides two complementary tools:
murmuration_plot()— shows the full weight distribution visually so you can see where the valley isperch()— quantifies match count, link rate, and clerical burden at each candidate cutoff, annotated with AIHW, WA Data Linkage Unit, and PHRN reference benchmarks
Example: running perch() on scored pairs
library(starling)
# Simulate a scored pairs object (bimodal: non-matches low, matches high)
set.seed(20260624L)
n_nonmatch <- 800
n_match <- 200
pairs_pred <- data.frame(
weights = c(
rnorm(n_nonmatch, mean = 5, sd = 3), # non-match cluster
rnorm(n_match, mean = 20, sd = 3) # match cluster
)
)
# Run the sensitivity sweep
results <- perch(
pairs_pred = pairs_pred,
n_records_df1 = 250L, # size of the primary dataset
thresholds = seq(5, 28, by = 1),
report = TRUE,
plot = FALSE # set TRUE in interactive session
)
#>
#> == starling::perch() ==========================================
#>
#> Total pairs : 1,000
#> Threshold range : 5.0 – 28.0 (step 1.00)
#> Clerical window : ±3.0 units (6-unit band)
#>
#> Threshold n_above pct (%) n_clerical clerical% link_rt% Reference
#> ------------------------------------------------------------------------------------
#> 5.0 597 59.7% 558 55.8% 238.8%
#> 6.0 486 48.6% 522 52.2% 194.4%
#> 7.0 398 39.8% 455 45.5% 159.2%
#> 8.0 313 31.3% 375 37.5% 125.2%
#> 9.0 268 26.8% 280 28.0% 107.2%
#> [*] 10.0 237 23.7% 199 19.9% 94.8% <- AIHW/WADLU: lower bound of clerical review zone
#> 11.0 222 22.2% 118 11.8% 88.8%
#> 12.0 206 20.6% 78 7.8% 82.4%
#> 13.0 199 19.9% 56 5.6% 79.6%
#> 14.0 195 19.5% 53 5.3% 78.0%
#> [*] 15.0 190 19.0% 54 5.4% 76.0% <- PHRN: lower bound for <0.5% false-match rate (full variable set)
#> 16.0 181 18.1% 72 7.2% 72.4%
#> [*] 17.0 169 16.9% 88 8.8% 67.6% <- starling default (balanced: Medicare + 2 names + DOB)
#> 18.0 152 15.2% 113 11.3% 60.8%
#> 19.0 127 12.7% 125 12.5% 50.8%
#> [*] 20.0 107 10.7% 129 12.9% 42.8% <- AIHW/WADLU: upper bound of clerical zone / confirmed matches above
#> 21.0 77 7.7% 129 12.9% 30.8%
#> 22.0 56 5.6% 111 11.1% 22.4%
#> 23.0 40 4.0% 100 10.0% 16.0%
#> 24.0 23 2.3% 77 7.7% 9.2%
#> 25.0 16 1.6% 56 5.6% 6.4%
#> 26.0 7 0.7% 40 4.0% 2.8%
#> 27.0 0 0.0% 23 2.3% 0.0%
#> 28.0 0 0.0% 16 1.6% 0.0%
#>
#> [*] Highlighted threshold (Australian/international reference)
#>
#> Key benchmarks:
#> 10–20 AIHW / WA Data Linkage Unit clerical review zone
#> Confirmed matches above ~20, non-matches below ~10,
#> marginal pairs sent for human adjudication.
#> 15–20 PHRN operational target for <0.5% false-match rate
#> (full variable set: Medicare + 2 names + DOB).
#> 17 starling default — balanced for SCPHU routine surveillance.The [*] markers in the table show the four reference
threshold values from AIHW, PHRN, and the starling default.
Example: visualising the score distribution
murmuration_plot(pairs_pred, threshold = 17, show_density = FALSE,
palette = "sch")
Linkage weight distribution. The threshold line should sit in the valley between the two clusters.
A well-behaved distribution is clearly bimodal. If yours is not:
- Unimodal at low scores: almost all pairs are non-matches. Check that your blocking variable is not too restrictive, excluding true match candidates from ever being compared.
- Unimodal at high scores: unusual. Suggests a highly specific blocking variable that only passes near-certain matches.
- Flat / no clear valley: the comparison variables
may lack discriminating power. Review
preflight()output for high missingness or poor name quality.
Using perch_before_linking = TRUE inside
murmuration()
For the common one-step workflow, set
perch_before_linking = TRUE in the
murmuration() call. The sensitivity table is printed
immediately after the EM model fits — before the threshold is applied —
giving you a chance to review before the linked dataset is produced.
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 = TRUE # <-- prints the perch() table mid-linkage
)In a Quarto document or batch job
(interactive() == FALSE), the table is printed to the
console but the plot is suppressed and execution continues
automatically. In an interactive session, the plot is displayed.
If the table suggests a better threshold (say 19), re-run
murmuration() with the new value. The EM model must be
re-fitted — the scored pairs are not stored between calls. This is by
design: murmuration() is a single-step function.
Clerical review zone
In a formal two-threshold design (AIHW/WADLU practice), pairs in the
10–20 zone are neither automatically accepted nor rejected — they are
sent for human review. The n_clerical column in
perch() output shows how many pairs fall in the review zone
around each candidate threshold (default window: ±3 units).
A large n_clerical relative to n_above
means the threshold sits in a high-density region — small threshold
changes would reclassify many pairs. This is a signal that you are
sitting on a slope, not in a valley.
Session information
sessionInfo()
#> R version 4.5.2 (2025-10-31 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26100)
#>
#> Matrix products: default
#> LAPACK version 3.12.1
#>
#> locale:
#> [1] LC_COLLATE=C LC_CTYPE=English_Australia.utf8
#> [3] LC_MONETARY=English_Australia.utf8 LC_NUMERIC=C
#> [5] LC_TIME=English_Australia.utf8
#>
#> time zone: Australia/Brisbane
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] starling_1.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 jsonlite_2.0.0 dplyr_1.2.0
#> [4] compiler_4.5.2 Rcpp_1.1.1 tidyselect_1.2.1
#> [7] stringr_1.6.0 parallel_4.5.2 jquerylib_0.1.4
#> [10] scales_1.4.0 yaml_2.3.12 fastmap_1.2.0
#> [13] ggplot2_4.0.3 R6_2.6.1 labeling_0.4.3
#> [16] generics_0.1.4 knitr_1.51 tibble_3.3.1
#> [19] lubridate_1.9.5 bslib_0.10.0 pillar_1.11.1
#> [22] RColorBrewer_1.1-3 rlang_1.2.0 stringi_1.8.7
#> [25] cachem_1.1.0 reclin2_0.6.0 xfun_0.56
#> [28] sass_0.4.10 S7_0.2.1 otel_0.2.0
#> [31] timechange_0.4.0 cli_3.6.5 withr_3.0.2
#> [34] magrittr_2.0.4 stringdist_0.9.17 digest_0.6.39
#> [37] grid_4.5.2 rstudioapi_0.18.0 lifecycle_1.0.5
#> [40] vctrs_0.7.1 lpSolve_5.6.23 data.table_1.18.2.1
#> [43] evaluate_1.0.5 glue_1.8.0 farver_2.1.2
#> [46] rmarkdown_2.31 tools_4.5.2 pkgconfig_2.0.3
#> [49] htmltools_0.5.9