# 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 is
- `perch()` — quantifies match count, link rate, and clerical burden at each candidate cutoff, annotated with AIHW, WA Data Linkage Unit, and PHRN reference benchmarks

## Australian linkage authority benchmarks

The three key Australian and international reference points:

Threshold range | Authority | Basis
--- | --- | ---
10–20 (clerical zone) | AIHW / WA Data Linkage Unit (WADLU) | Standard two-threshold practice: confirmed matches above ~20, confirmed non-matches below ~10, marginal pairs sent for human adjudication
~15–20 | Population Health Research Network (PHRN) | Operational target: false-match rate < 0.5% with a full variable set (Medicare + 2 names + DOB)
17 | starling default | Balanced starting point for SCPHU routine surveillance linkage

These are reference points, not rules. The correct threshold for your dataset is the one that sits in the valley of your score distribution.

## 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