Ranking speakers by Jensen-Shannon distance and checking Pillai agreement

This vignette walks through the measurement protocol that phontrast implements in rank_contrasts() (since 2.5.0; with the study’s calibrated estimator since 2.5.1). The protocol comes from a simulation study that scored seven vowel-overlap measures against a known ground truth (Berry, under review, Sec. VII.A). Its recommendations, in three steps:

Compute Jensen-Shannon distance ($\sqrt{JSD}$) and Pillai on the same tokens, and report the estimated shared probability mass beside them.
Rank speakers by $\sqrt{JSD}$.
Flag speakers whose Pillai percentile rank differs from their $\sqrt{JSD}$ percentile rank by 0.25 of the ordering or more, inspect them, and plot them.

Each step carries conditions: measurements at the $\sqrt{JSD}$ ceiling are set apart, sample-size floors decide whether a rank or a flag may be read, and a bandwidth check marks measurements whose rank depends on the smoothing. rank_contrasts() applies and reports all of them.

The example cohort

vowel_cohort is a simulated twelve-speaker cohort built so that each part of the protocol has something to do: eight speakers with a graded centroid gap, one speaker whose two vowels share a centroid but differ in shape (a bimodal “eh”), one speaker with fully separated vowels, and two speakers with fewer tokens than the others.

library(phontrast)

head(vowel_cohort)
#>   speaker vowel    f1     f2
#> 1   spk01    ih 549.2 1864.0
#> 2   spk01    ih 507.4 1966.1
#> 3   spk01    ih 515.8 1598.3
#> 4   spk01    ih 500.9 2007.8
#> 5   spk01    ih 553.0 1721.5
#> 6   spk01    ih 606.3 1666.8
table(vowel_cohort$speaker, vowel_cohort$vowel)
#>        
#>          eh  ih
#>   spk01 100 100
#>   spk02 100 100
#>   spk03 100 100
#>   spk04 100 100
#>   spk05 100 100
#>   spk06 100 100
#>   spk07 100 100
#>   spk08 100 100
#>   spk09 100 100
#>   spk10 100 100
#>   spk11  60  60
#>   spk12  40  40

Steps 1 to 3 in one call

ranking <- rank_contrasts(
  data = vowel_cohort,
  features = c("f1", "f2"),
  category_col = "vowel",
  group_col = "speaker"
)
ranking
#> <phontrast ranking: 12 speakers, 2 features (f1, f2)>
#> sqrt(JSD), Pillai, and shared mass from the same tokens; kernel: legacy estimator, Hpi / ks, all pooled tokens, no leave-one-out.
#> Ranked: 11 (1 at ceiling, sqrt(JSD) >= 0.99). Floors at d = 2: rank from 50, flag from 100 tokens per category.
#> Below rank floor (rank with Pillai): 1. Flag not readable: 2. Flagged (|rank_diff| >= 0.25): 1 -> spk09.
#> Bandwidth check (pooled scott.diag x0.5 / x2): 2 set aside -> spk09, spk06.
#> 
#> # A tibble: 12 × 19
#>    group n_tokens n_min sqrt_jsd   pillai shared_mass at_ceiling  pr_jsd
#>    <chr>    <int> <int>    <dbl>    <dbl>       <dbl> <lgl>        <dbl>
#>  1 spk10      200   100    1     0.993      2.23e-307 TRUE       NA     
#>  2 spk08      200   100    0.968 0.811      3.98e-  2 FALSE       0.955 
#>  3 spk07      200   100    0.967 0.794      2.80e-  2 FALSE       0.864 
#>  4 spk09      200   100    0.905 0.000248   7.73e-  2 FALSE       0.773 
#>  5 spk12       80    40    0.829 0.689      1.79e-  1 FALSE       0.682 
#>  6 spk05      200   100    0.788 0.637      2.03e-  1 FALSE       0.591 
#>  7 spk06      200   100    0.788 0.665      1.98e-  1 FALSE       0.5   
#>  8 spk11      120    60    0.625 0.545      3.87e-  1 FALSE       0.409 
#>  9 spk04      200   100    0.582 0.494      4.42e-  1 FALSE       0.318 
#> 10 spk03      200   100    0.560 0.490      4.71e-  1 FALSE       0.227 
#> 11 spk02      200   100    0.326 0.206      6.67e-  1 FALSE       0.136 
#> 12 spk01      200   100    0.192 0.0513     8.12e-  1 FALSE       0.0455
#> # ℹ 11 more variables: pr_pillai <dbl>, rank_diff <dbl>, flag <lgl>,
#> #   rank_licensed <lgl>, flag_licensed <lgl>, rank_basis <chr>,
#> #   sqrt_jsd_half <dbl>, sqrt_jsd_double <dbl>, bw_shift <dbl>,
#> #   sign_change <lgl>, set_aside <lgl>

The header states what was computed and under which conditions; the table has one row per speaker, sorted by $\sqrt{JSD}$. Reading it column by column:

sqrt_jsd, pillai, shared_mass: step 1. $\sqrt{JSD}$ and shared mass are read off one kernel density estimate per speaker, Pillai off the same tokens. The kernel estimator is the one the protocol was calibrated on (see “The estimator behind the numbers” below).
at_ceiling: TRUE where $\sqrt{JSD} \ge 0.99$. spk10’s vowels are fully separated, so its kernel estimate can no longer grade separation; its shared mass and Pillai are still reported, but it leaves both rankings.
pr_jsd, pr_pillai, rank_diff: step 2, on the percentile-rank scale $(\bar r - \tfrac12)/k$ over the $k$ speakers below the ceiling (ties averaged).
flag: step 3, TRUE where abs(rank_diff) >= 0.25. spk09 is the planted disagreement: Pillai, a mean-based measure, sees almost no contrast between two vowels with the same centroid, while $\sqrt{JSD}$ sees two distributions that barely overlap.

ranking[ranking$flag %in% TRUE, c("group", "sqrt_jsd", "pillai", "pr_jsd", "pr_pillai", "rank_diff")]
#> # A tibble: 1 × 6
#>   group sqrt_jsd   pillai pr_jsd pr_pillai rank_diff
#>   <chr>    <dbl>    <dbl>  <dbl>     <dbl>     <dbl>
#> 1 spk09    0.905 0.000248  0.773    0.0455     0.727

Licensing: when a rank or a flag may be read

The simulation licenses these readings only from certain token counts per vowel per speaker, and only up to certain dimensionalities. The floors are applied per speaker from the smaller vowel’s count (n_min):

protocol_floors(2)
#> $d
#> [1] 2
#> 
#> $rank_floor
#> [1] 50
#> 
#> $flag_floor
#> [1] 100
#> 
#> $calibrated_at
#> [1] 2
#> 
#> $ordering_only
#> [1] FALSE
ranking[, c("group", "n_min", "rank_licensed", "flag_licensed", "rank_basis", "flag")]
#> # A tibble: 12 × 6
#>    group n_min rank_licensed flag_licensed rank_basis flag 
#>    <chr> <int> <lgl>         <lgl>         <chr>      <lgl>
#>  1 spk10   100 TRUE          TRUE          sqrt_jsd   NA   
#>  2 spk08   100 TRUE          TRUE          sqrt_jsd   FALSE
#>  3 spk07   100 TRUE          TRUE          sqrt_jsd   FALSE
#>  4 spk09   100 TRUE          TRUE          sqrt_jsd   TRUE 
#>  5 spk12    40 FALSE         FALSE         pillai     NA   
#>  6 spk05   100 TRUE          TRUE          sqrt_jsd   FALSE
#>  7 spk06   100 TRUE          TRUE          sqrt_jsd   FALSE
#>  8 spk11    60 TRUE          FALSE         sqrt_jsd   NA   
#>  9 spk04   100 TRUE          TRUE          sqrt_jsd   FALSE
#> 10 spk03   100 TRUE          TRUE          sqrt_jsd   FALSE
#> 11 spk02   100 TRUE          TRUE          sqrt_jsd   FALSE
#> 12 spk01   100 TRUE          TRUE          sqrt_jsd   FALSE

At two dimensions, ranking by $\sqrt{JSD}$ is licensed from 50 tokens per vowel and the flag is readable from 100. spk11 (60 tokens) can be ranked by $\sqrt{JSD}$ but its flag is NA: a single speaker’s rank difference spans about 0.28 of the ordering across redraws at 50 tokens, so a 0.25 margin cannot be read there. spk12 (40 tokens) is below the rank floor, so rank_basis says to order that speaker by Pillai. rank_diff is still reported for both; only the inference is withheld. From eight dimensions the agreement test returns agreement whatever the data contain, and above eight the study offers ordering evidence only, so rank_contrasts() withholds the flag and the licence there too.

The agreement plot

plot() on the ranking draws step 3: Pillai percentile rank against $\sqrt{JSD}$ percentile rank, the identity line, and the inspection band of $\pm 0.25$ around it. Flagged speakers are coloured and labelled, set-aside speakers crossed, and speakers whose flag is not readable are hollow.

plot(ranking)

Inspecting a flagged speaker across bandwidths

Before interpreting a flag, the protocol asks for one more check: recompute $\sqrt{JSD}$ at half and at twice the diagonal Scott bandwidth (selected once on the pooled tokens of the pair, as in the study), re-rank, and set the measurement aside if its rank moves by 0.25 of the ordering or more, or if a flagged rank difference changes sign. rank_contrasts() runs this check by default (bw_check = TRUE):

ranking[, c("group", "sqrt_jsd_half", "sqrt_jsd_double", "bw_shift", "sign_change", "set_aside")]
#> # A tibble: 12 × 6
#>    group sqrt_jsd_half sqrt_jsd_double bw_shift sign_change set_aside
#>    <chr>         <dbl>           <dbl>    <dbl> <lgl>       <lgl>    
#>  1 spk10         1               0.979  NA      FALSE       NA       
#>  2 spk08         0.974           0.714  -0.0909 FALSE       FALSE    
#>  3 spk07         0.972           0.673  -0.0909 FALSE       FALSE    
#>  4 spk09         0.981           0.547   0.273  FALSE       TRUE     
#>  5 spk12         0.898           0.496   0.182  FALSE       FALSE    
#>  6 spk05         0.846           0.515   0      FALSE       FALSE    
#>  7 spk06         0.824           0.559  -0.273  FALSE       TRUE     
#>  8 spk11         0.715           0.398   0      FALSE       FALSE    
#>  9 spk04         0.663           0.388   0.0909 FALSE       FALSE    
#> 10 spk03         0.637           0.391  -0.0909 FALSE       FALSE    
#> 11 spk02         0.441           0.220   0      FALSE       FALSE    
#> 12 spk01         0.311           0.120   0      FALSE       FALSE

Here the check sets aside spk09 and spk06. spk09’s separation lives in the shape of its bimodal “eh”, which twice the Scott bandwidth smooths away, so its $\sqrt{JSD}$ rank falls when the bandwidth doubles. The flag stands, but the protocol asks for the distributions to be inspected before it is interpreted.

inspect_contrast() then redraws one speaker with plot_contrast()’s layers, one panel per bandwidth, so the flag can be judged against the distributions that produced it. The panels are labelled with $\sqrt{JSD}$ and shared mass at that bandwidth and with the speaker’s Pillai, and the subtitle restates the reported ranks and the outcome of the bandwidth check.

inspect_contrast(ranking, "spk09", reverse_x = TRUE, reverse_y = TRUE)

The same bandwidth multiplier is available throughout the kernel path as bw_scale (on jsd_kde_nd(), estimate_jsd(), phontrast(), plot_contrast(), and the overlap functions), so any kernel estimate can be bracketed by hand.

The estimator behind the numbers

The floors, the margin, the ceiling, and the bandwidth check above were calibrated on one kernel estimator: phontrast’s method = "legacy", which evaluates both vowels’ kernel densities at the pooled tokens and reads $\sqrt{JSD}$ off the two self-normalized density vectors, with no leave-one-out correction. Its bandwidth rule, engine, and evaluation subsample change with dimensionality. recommended_estimator() returns the settings, and rank_contrasts() uses them by default:

recommended_estimator(2)[c("tier", "method", "bw", "engine", "eval_on", "eval_n", "bracket_bw")]
#> $tier
#> [1] "d <= 4"
#> 
#> $method
#> [1] "legacy"
#> 
#> $bw
#> [1] "Hpi"
#> 
#> $engine
#> [1] "ks"
#> 
#> $eval_on
#> [1] "pooled"
#> 
#> $eval_n
#> NULL
#> 
#> $bracket_bw
#> [1] "scott.pooled"
recommended_estimator(8)[c("tier", "method", "bw", "engine", "eval_on", "eval_n", "bracket_bw")]
#> $tier
#> [1] "5 <= d <= 13"
#> 
#> $method
#> [1] "legacy"
#> 
#> $bw
#> [1] "scott.diag"
#> 
#> $engine
#> [1] "fast_diag"
#> 
#> $eval_on
#> [1] "pooled"
#> 
#> $eval_n
#> [1] 200
#> 
#> $bracket_bw
#> [1] "scott.pooled"

phontrast’s default estimator elsewhere, method = "mc", targets the continuous divergence and reads higher on overlapping vowels, so it moves contrasts across the ceiling and through the bandwidth check. Rankings computed with it are not covered by the study’s calibration.

To run the protocol with other settings, pass a modified list:

rank_contrasts(
  vowel_cohort, c("f1", "f2"), "vowel", "speaker",
  estimator = modifyList(recommended_estimator(2), list(bw = "scott.diag", engine = "fast_diag"))
)

Measures and estimators side by side

The study’s central distinction is between a measure (the quantity) and its estimator. phontrast() makes that distinction visible: the closed-form Gaussian Bhattacharyya columns and the matched-kernel ones estimate the same quantity under different models, and the whole kernel family (Jensen-Shannon, overlap, total variation, kernel Bhattacharyya, Hellinger) is scored on one shared density estimate per comparison.

one_speaker <- vowel_cohort[vowel_cohort$speaker == "spk05", ]
phontrast(
  one_speaker, c("f1", "f2"), "vowel",
  metrics = c("js_distance", "pillai", "overlap", "tv",
              "bhattacharyya", "bhattacharyya_kde", "euclidean")
)
#> # A tibble: 1 × 13
#>   scope  n_tokens pillai pillai_p_value bhatt_dist bhatt_affinity js_distance
#>   <chr>     <int>  <dbl>          <dbl>      <dbl>          <dbl>       <dbl>
#> 1 global      200  0.637       4.52e-44      0.881          0.415       0.826
#> # ℹ 6 more variables: percent_overlap <dbl>, total_variation <dbl>,
#> #   bhatt_kde_dist <dbl>, bhatt_kde_affinity <dbl>, hellinger <dbl>,
#> #   euclidean_dist <dbl>

What to report

For a set of speakers measured on one contrast, report per speaker the three step-1 quantities (sqrt_jsd, pillai, shared_mass), the token count per vowel (n_min) and the number of features, and then the $\sqrt{JSD}$ ordering with its licensing. Name the flagged speakers, say whether the flag was readable at that sample size, and show the inspected distributions for the ones you discuss. State the estimator settings (attr(ranking, "protocol")$estimator) and the bandwidth-check outcome. print(ranking) gives all of this in one block; the plots above give the picture.