Package {phontrast}

phontrast: Contrast and Separation Metrics for Phonological Categories

Description

A unified toolkit for quantifying the separation and overlap between phonological categories (e.g., vowels, consonants) in arbitrary n-dimensional acoustic spaces such as formants, MFCCs, spectral features, duration, or learned embeddings. The main entry point, phontrast(), computes and compares multiple contrast metrics in one call: Jensen-Shannon divergence and distance, the Pillai-Bartlett trace, Bhattacharyya distance and affinity, Mahalanobis distance, and proportional overlap, globally or by group, with optional bootstrap intervals. Percent overlap is returned as a 0–1 proportion, not a 0–100 percentage.

Details

phontrast was formerly released as phonJSD (through version 1.2.0); see the package NEWS for migration notes.

Recommended workflow

1. Start with phontrast() to compute one or more contrast metrics for a two-category contrast, globally or by group, in tidy long or wide form.

2. Use estimate_jsd() when Jensen-Shannon divergence or Jensen-Shannon distance is the primary outcome and you need optional bootstrap intervals.

3. Use lower-level helpers such as jsd_kde_nd(), percent_overlap_kde(), pillai_overlap(), and bhattacharyya_mvnorm() when validating methods, debugging one contrast, or reproducing a specific metric.

4. Use plot_contrast() for a distribution-aware, annotated view of one contrast that draws the same density model the metrics use, and plot_overlap_metrics() (also available as plot() / ggplot2::autoplot() on phontrast() results), plot_category_space(), and plot_category_pca() for ggplot2-backed diagnostics and presentation figures, all sharing the colorblind-safe theme_phontrast() visual identity.

5. To rank a set of speakers on one contrast, use rank_contrasts() (see the next section), then plot() the ranking and inspect_contrast() any flagged speaker.

Ranking protocol

rank_contrasts() implements the measurement protocol recommended by the simulation study behind this package (Berry, under review, Sec. VII.A; see citation("phontrast")): compute Jensen-Shannon distance and Pillai on the same tokens and report shared mass beside them; rank speakers by Jensen-Shannon distance on the percentile-rank scale of percentile_rank(); and flag speakers whose Pillai and Jensen-Shannon percentile ranks differ by 0.25 of the ordering or more. Measurements at the Jensen-Shannon ceiling are set apart, sample-size floors (protocol_floors()) decide per speaker whether a rank or a flag may be read, and a bandwidth check marks measurements whose rank depends on the smoothing. recommended_estimator() supplies the kernel settings the study used at each dimensionality; plot_rank_agreement() and inspect_contrast() draw the ranking and a flagged speaker; the bundled vowel_cohort data and the vignette "Ranking speakers by Jensen-Shannon distance and checking Pillai agreement" walk through it.

Choosing metrics

JSD, Jensen-Shannon distance, Pillai trace, Bhattacharyya distance, and Mahalanobis distance increase as categories become more separated. Percent overlap and Bhattacharyya affinity increase as categories overlap more. The long output from phontrast() includes an orientation column and a separation-oriented separation_value column to make these directions explicit. JSD and percent overlap estimate distributional separation/overlap using KDE by default; Pillai and Mahalanobis emphasize mean separation; Bhattacharyya metrics use a multivariate-normal approximation. Opt-in metrics add total variation ("tv"), the matched-kernel Bhattacharyya and Hellinger distances read off the same kernel densities as JSD ("bhattacharyya_kde"), and the Euclidean distance between standardized category means ("euclidean").

Density backends

The distributional metrics (Jensen-Shannon divergence and proportional overlap) are computed from a density estimate for each category. The density argument selects that estimate: "kde" (the default) uses kernel density estimation, and "mvnorm" fits one multivariate normal per category and estimates the metric between the two Gaussians by Monte-Carlo (with mc_n samples and reproducible eval_seed). The "mvnorm" backend matches the estimator behind JSD and overlap to the same multivariate-normal assumptions the Pillai, Bhattacharyya, and Mahalanobis columns already make, and is convenient for higher-dimensional feature spaces where multivariate KDE is impractical. It is available on phontrast(), estimate_jsd(), estimate_overlap(), jsd_summary(), global_boot_jsd(), jsd_kde_nd(), and percent_overlap_kde(); the parametric metrics are unaffected by it.

High-dimensional workflows

Metrics can be estimated in arbitrary n-dimensional numeric feature spaces, including MFCCs and learned embeddings. Use plot_category_pca() for a two-dimensional PCA diagnostic, but report metric estimates from the intended full feature set.

Confidence intervals

Confidence intervals use ci_lower and ci_upper columns. Legacy JSD aliases jsd_low and jsd_high are retained for compatibility.

Author(s)

Maintainer: Grant M. Berry berry.grant@gmail.com

See Also

Useful links:

• https://github.com/berrygrant/phontrast

• Report bugs at https://github.com/berrygrant/phontrast/issues

Bhattacharyya distance and affinity under multivariate normality

Description

Estimates means and covariances per category and computes Bhattacharyya distance and affinity (exp(-distance)) under the assumption of multivariate normality.

Usage

bhattacharyya_mvnorm(data, features, category_col, eps = 1e-06)

Arguments

Value

A list with distance and affinity (exp(-distance)).

Bootstrap JSD for each group

Description

Computes bootstrap mean, SD, and confidence interval for JSD within each group (e.g., speaker), using resampling with replacement.

Usage

boot_jsd(
  data,
  group_col,
  category_col,
  features,
  n_boot = 1000,
  min_tokens = 20,
  est_distance = FALSE,
  conf_level = 0.95,
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  eval_on = c("pooled", "group1", "group2", "pooled_sample"),
  eval_n = NULL,
  eval_seed = NULL,
  engine = c("ks", "fast_diag", "fast_diagonal"),
  chunk_size = 1000L,
  ...
)

Arguments

Value

A tibble with one row per group and columns: group, n_tokens, n_boot, conf_level, jsd_mean, jsd_sd, ci_lower, ci_upper, jsd_low, and jsd_high. jsd_low and jsd_high are retained as legacy aliases for ci_lower and ci_upper.

Compare phonological category overlap metrics (deprecated)

Description

compare_overlap_metrics() was renamed to phontrast() in phontrast 2.0.0 (the package formerly released as 'phonJSD'). It remains as a thin wrapper that calls phontrast() with output = "wide" for backward compatibility and will be removed in a future release. New code should call phontrast().

Usage

compare_overlap_metrics(
  data,
  features,
  category_col,
  group_col = NULL,
  min_tokens = 20,
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  eval_on = c("pooled", "group1", "group2", "pooled_sample"),
  eval_n = NULL,
  eval_seed = NULL,
  engine = c("ks", "fast_diag", "fast_diagonal"),
  chunk_size = 1000L,
  eps = 1e-06,
  output = c("wide", "long"),
  do_boot = FALSE,
  n_boot = 1000,
  conf_level = 0.95,
  progress = TRUE,
  method = c("mc", "legacy"),
  density = c("kde", "mvnorm"),
  mc_n = 10000L,
  bw_scale = 1
)

Arguments

Value

See phontrast(); wide format by default.

See Also

phontrast()

Estimate Bhattacharyya distance, globally or by group

Description

Unified front-end for Bhattacharyya distance (and affinity) under a multivariate normal approximation.

Usage

estimate_bhatt(
  data,
  features,
  category_col,
  group_col = NULL,
  min_tokens = 20,
  eps = 1e-06
)

Arguments

Value

Data frame with either one global row or one row per group.

Estimate Jensen-Shannon divergence or distance between two categories

Description

Use this function when Jensen-Shannon divergence (JSD) or Jensen-Shannon distance is the primary outcome. If you want to compare JSD with Pillai, Bhattacharyya, Mahalanobis, and percent-overlap metrics, start with phontrast() instead.

Usage

estimate_jsd(
  data,
  features,
  category_col,
  group_col = NULL,
  do_boot = FALSE,
  n_boot = 1000,
  min_tokens = 20,
  est_distance = FALSE,
  conf_level = 0.95,
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  eval_on = c("pooled", "group1", "group2", "pooled_sample"),
  eval_n = NULL,
  eval_seed = NULL,
  engine = c("ks", "fast_diag", "fast_diagonal"),
  chunk_size = 1000L,
  method = c("mc", "legacy"),
  density = c("kde", "mvnorm"),
  mc_n = 10000L,
  bw_scale = 1,
  ...
)

Arguments

Details

JSD is bounded from 0 to 1 for two equally weighted distributions. Larger values indicate greater category separation. Jensen-Shannon distance is sqrt(JSD) and has the same direction of interpretation.

Value

A tibble. Global: one row with columns scope, n_tokens, n_boot, conf_level, jsd_point, jsd_mean, jsd_sd, ci_lower, ci_upper, jsd_low, jsd_high. Grouped: one row per group with columns scope, group, n_tokens, n_boot, conf_level, jsd_point, jsd_mean, jsd_sd, ci_lower, ci_upper, jsd_low, jsd_high. ci_lower and ci_upper are the preferred confidence interval columns; jsd_low and jsd_high are retained as legacy aliases.

Examples

set.seed(2026)
vowels <- data.frame(
  speaker = rep(c("s01", "s02"), each = 60),
  vowel = rep(rep(c("ih", "eh"), each = 30), 2),
  f1 = c(
    rnorm(30, 500, 55), rnorm(30, 560, 60),
    rnorm(30, 510, 60), rnorm(30, 575, 65)
  ),
  f2 = c(
    rnorm(30, 1980, 150), rnorm(30, 1880, 155),
    rnorm(30, 1960, 160), rnorm(30, 1840, 165)
  )
)

# Point estimate of JSD (fast), globally and by speaker.
estimate_jsd(
  data = vowels,
  features = c("f1", "f2"),
  category_col = "vowel"
)
estimate_jsd(vowels, c("f1", "f2"), "vowel", group_col = "speaker")

# Bootstrap confidence intervals, shown on a single feature: multivariate
# bootstraps work the same way but repeat multivariate bandwidth selection
# on every resample, so they take correspondingly longer. Increase n_boot
# for real analyses.
estimate_jsd(vowels, "f1", "vowel", do_boot = TRUE, n_boot = 20)

Estimate proportional overlap globally or by group

Description

Unified front-end for KDE-based proportional overlap between two categories. The returned overlap column is a 0–1 proportion, not a 0–100 percentage.

Usage

estimate_overlap(
  data,
  features,
  category_col,
  group_col = NULL,
  min_tokens = 20,
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  eval_on = c("pooled", "group1", "group2", "pooled_sample"),
  eval_n = NULL,
  eval_seed = NULL,
  engine = c("ks", "fast_diag", "fast_diagonal"),
  chunk_size = 1000L,
  method = c("mc", "legacy"),
  density = c("kde", "mvnorm"),
  mc_n = 10000L,
  bw_scale = 1,
  ...
)

Arguments

Value

A tibble (global = one row; grouped = one per group) with overlap as a 0–1 proportion.

Estimate Pillai trace, globally or by group

Description

Unified front-end for Pillai-Bartlett trace. If group_col is NULL, computes a single global Pillai trace for the full dataset. If group_col is provided, computes Pillai per group (e.g., per speaker).

Usage

estimate_pillai(
  data,
  features,
  category_col,
  group_col = NULL,
  min_tokens = 20
)

Arguments

Value

A data frame with either one global row or one row per group.

Extract MFCCs for vowel segments

Description

Adds MFCC feature columns (e.g., mfcc1..mfcc13) to a data frame that contains audio file paths and (optionally) segment boundaries.

Usage

extract_mfcc(
  data,
  file_col,
  start_col = NULL,
  end_col = NULL,
  fs = NULL,
  numcep = 13,
  prefix = "mfcc",
  strict = FALSE,
  warn = TRUE,
  ...
)

Arguments

Details

This function uses tuneR to read WAV files and compute MFCCs. The package is optional; if it is not installed, an informative error is raised.

Value

The input data frame with added MFCC columns.

Examples

if (requireNamespace("tuneR", quietly = TRUE)) {
  wav_path <- tempfile(fileext = ".wav")
  wave <- tuneR::sine(
    freq = 440, duration = 0.25, samp.rate = 16000, xunit = "time"
  )
  tuneR::writeWave(wave, wav_path)

  segments <- data.frame(wav_path = wav_path)
  mfcc <- extract_mfcc(segments, file_col = "wav_path", numcep = 3)
  unlink(wav_path)
  mfcc[, c("mfcc1", "mfcc2", "mfcc3")]
}

Global JSD with bootstrap confidence interval

Description

Computes a single Jensen-Shannon divergence (JSD) value for two categories in an n-dimensional acoustic space, together with bootstrap-based confidence intervals obtained by resampling tokens with replacement.

Usage

global_boot_jsd(
  data,
  features,
  category_col,
  n_boot = 1000,
  min_tokens = 20,
  est_distance = FALSE,
  conf_level = 0.95,
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  eval_on = c("pooled", "group1", "group2", "pooled_sample"),
  eval_n = NULL,
  eval_seed = NULL,
  engine = c("ks", "fast_diag", "fast_diagonal"),
  chunk_size = 1000L,
  method = c("mc", "legacy"),
  density = c("kde", "mvnorm"),
  mc_n = 10000L,
  ...
)

Arguments

Details

This is the "group-wise" version of JSD: it ignores speakers and treats all tokens as coming from a single population for each category.

Value

A one-row data frame with columns:

• n_tokens - total number of tokens used

• n_boot - number of successful bootstrap samples

• conf_level - confidence level used for the interval

• jsd_point - JSD on the full dataset

• jsd_mean - mean JSD across bootstrap samples

• jsd_sd - standard deviation of bootstrap JSD

• ci_lower, ci_upper - bootstrap confidence interval

• jsd_low, jsd_high - legacy aliases for ci_lower and ci_upper

Global Pillai trace (point estimate)

Description

Computes a single Pillai-Bartlett trace for the full dataset in the specified feature space, returning a one-row data frame that includes the total number of tokens used.

Usage

global_pillai(data, features, category_col, min_tokens = 20)

Arguments

Value

A one-row data frame with columns: n_tokens, pillai, and p_value.

Hierarchical bootstrap for JSD-based models

Description

Performs a hierarchical bootstrap: resample groups with replacement, resample tokens within each sampled group, compute JSD per group, fit a model to the bootstrap JSD values, and repeat.

Usage

hier_boot_jsd_model(
  data,
  group_col,
  category_col,
  features,
  formula,
  fit_fun = NULL,
  n_outer = 200,
  min_tokens = 20,
  eps = 1e-06,
  progress = TRUE,
  ...
)

Arguments

Details

This lets you propagate measurement uncertainty in JSD into model parameters (e.g., GAM/LMM coefficients).

Value

A tibble with columns:

• boot_id - bootstrap replicate index

• term - model term

• estimate - estimate for that term in that replicate

Examples

set.seed(2026)
speakers <- paste0("s", 1:4)
dat <- data.frame(
  speaker = rep(speakers, each = 60),
  age = rep(c(22, 35, 48, 61), each = 60),
  vowel = rep(rep(c("ih", "eh"), each = 30), 4)
)
dat$f1 <- rnorm(
  nrow(dat),
  mean = ifelse(dat$vowel == "ih", 500, 560) + dat$age * 0.3,
  sd = 55
)

hier_boot_jsd_model(
  data = dat,
  group_col = "speaker",
  category_col = "vowel",
  features = "f1",
  formula = jsd_beta ~ age,
  fit_fun = stats::lm,
  n_outer = 3,
  min_tokens = 20,
  progress = FALSE
)

Inspect one speaker's contrast across bandwidths

Description

Redraws a single speaker from a rank_contrasts() ranking with plot_contrast()'s distribution-aware layers, one panel per kernel bandwidth, so a flagged or set-aside speaker can be inspected the way the protocol asks: at half, at the selected, and at twice the diagonal Scott bandwidth (the smoothing used by the bandwidth check). Each panel is labelled with the Jensen-Shannon distance and shared mass at that bandwidth and with the speaker's Pillai trace, which takes no bandwidth; the subtitle restates the reported ranks, the rank difference, and the outcome of the bandwidth check. The tokens come from the ranking itself (attr(ranking, "protocol")$data), so no further data is needed.

Usage

inspect_contrast(
  ranking,
  group,
  features = NULL,
  bw_scales = c(0.5, 1, 2),
  levels = c(0.5, 0.8, 0.95),
  points = TRUE,
  overlap = TRUE,
  grid_n = NULL,
  point_alpha = 0.55,
  point_size = 1.6,
  reverse_x = FALSE,
  reverse_y = FALSE
)

Arguments

Value

A ggplot2 plot object. The per-panel metrics are attached as attr(p, "inspect_metrics").

See Also

rank_contrasts(), plot_rank_agreement(), plot_contrast().

Examples

set.seed(2026)
gaps <- c(30, 60, 90, 120, 150, 180, 210, 240)
cohort <- do.call(rbind, lapply(seq_along(gaps), function(i) {
  data.frame(
    speaker = sprintf("s%02d", i),
    vowel = rep(c("ih", "eh"), each = 50),
    f1 = c(rnorm(50, 500, 55), rnorm(50, 500 + gaps[i], 55)),
    f2 = c(rnorm(50, 1900, 120), rnorm(50, 1900 - gaps[i], 120))
  )
}))
ranking <- rank_contrasts(cohort, c("f1", "f2"), "vowel", "speaker",
                          bw_check = FALSE)
if (requireNamespace("ggplot2", quietly = TRUE)) {
  inspect_contrast(ranking, "s04")
}

Jensen-Shannon divergence for discrete distributions

Description

Computes JSD(p, q) in bits for discrete probability vectors. The value is bounded in [0, 1] for equally weighted mixtures.

Usage

jsd(p, q)

Arguments

Value

A single numeric value: the Jensen-Shannon divergence in bits.

n-dimensional JSD via multivariate kernel density estimation

Description

Computes Jensen-Shannon divergence between two categories in an arbitrary n-dimensional acoustic space using multivariate KDE. The default engine uses the ks package; a faster diagonal-Gaussian engine is available for diagonal bandwidths.

Usage

jsd_kde_nd(
  data,
  features,
  group = "category",
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  eval_on = c("pooled", "group1", "group2", "pooled_sample"),
  eval_n = NULL,
  eval_seed = NULL,
  engine = c("ks", "fast_diag", "fast_diagonal"),
  chunk_size = 1000L,
  method = c("mc", "legacy"),
  density = c("kde", "mvnorm"),
  mc_n = 10000L,
  loo = TRUE,
  bw_scale = 1
)

Arguments

Details

By default (method = "mc") JSD is estimated with a Monte-Carlo plug-in: each category's KDE is evaluated at that category's own observations and the log density ratio against the mixture is averaged. This is a consistent estimator of the continuous JSD in any dimension. method = "legacy" reproduces the pre-1.2.0 self-normalized sample-point estimate (a bounded relative separation index rather than the continuous JSD); use it only to reproduce results from phonJSD 1.0.0.

Value

A single numeric JSD value in bits, bounded in [0, 1].

Examples

set.seed(2026)
vowels <- data.frame(
  vowel = rep(c("ih", "eh"), each = 40),
  f1 = c(rnorm(40, 500, 55), rnorm(40, 565, 60)),
  f2 = c(rnorm(40, 1980, 150), rnorm(40, 1870, 155))
)

# One-dimensional JSD, for example a single formant or duration.
jsd_kde_nd(vowels, features = "f1", group = "vowel")

# Two-dimensional JSD in F1/F2 space.
jsd_kde_nd(vowels, features = c("f1", "f2"), group = "vowel")

# Faster high-dimensional path: diagonal Scott bandwidth and sampled
# pooled evaluation points.
jsd_kde_nd(
  vowels,
  features = c("f1", "f2"),
  group = "vowel",
  bw = "scott.diag",
  eval_n = 40,
  eval_seed = 2026,
  engine = "fast_diag"
)

JSD summary: point estimate and optional bootstrap per group

Description

Convenience wrapper that returns both the point-estimate JSD and, optionally, bootstrap-based uncertainty (mean, SD, and CI) for each group.

Usage

jsd_summary(
  data,
  group_col,
  category_col,
  features,
  do_boot = TRUE,
  n_boot = 1000,
  min_tokens = 20,
  conf_level = 0.95,
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  eval_on = c("pooled", "group1", "group2", "pooled_sample"),
  eval_n = NULL,
  eval_seed = NULL,
  engine = c("ks", "fast_diag", "fast_diagonal"),
  chunk_size = 1000L,
  method = c("mc", "legacy"),
  density = c("kde", "mvnorm"),
  mc_n = 10000L,
  ...
)

Arguments

Value

A tibble with one row per group and columns:

• group - group ID (e.g., speaker)

• n_tokens - number of tokens for that group

• jsd_point - single JSD point estimate

• n_boot, conf_level, jsd_mean, jsd_sd, ci_lower, ci_upper, jsd_low, jsd_high - bootstrap summary columns. These are NA (or 0 for n_boot) if do_boot = FALSE.

Kullback-Leibler divergence for discrete distributions

Description

Computes KL(p || q) in bits for discrete probability vectors. Zero-probability events in p contribute zero; positive mass in p where q is zero returns Inf.

Usage

kl_div(p, q)

Arguments

Value

A single numeric value: the KL divergence in bits.

Proportional overlap between two distributions via KDE

Description

Computes the proportional overlap (shared area) between two categories in an n-dimensional acoustic space using multivariate kernel density estimation. Despite the historical function name, the return value is a 0–1 proportion: 0 = no overlap, 1 = identical.

Usage

percent_overlap_kde(
  data,
  features,
  category_col,
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  eval_on = c("pooled", "group1", "group2", "pooled_sample"),
  eval_n = NULL,
  eval_seed = NULL,
  engine = c("ks", "fast_diag", "fast_diagonal"),
  chunk_size = 1000L,
  method = c("mc", "legacy"),
  density = c("kde", "mvnorm"),
  mc_n = 10000L,
  bw_scale = 1,
  ...
)

Arguments

Value

Numeric scalar proportion in [0, 1].

Percentile ranks with averaged ties

Description

Converts a vector of measurements into percentile ranks, (\bar{r} - 1/2) / k, where \bar{r} is the average rank (ties share the mean of the ranks they span) and k is the number of measurements ranked. This is the rank scale on which the agreement check of rank_contrasts() compares Jensen-Shannon distance with Pillai: a difference of 0.25 means the two measures place a speaker a quarter of the ordering apart. Ranks are ascending, so larger values of x receive larger percentile ranks. Measurements that are NA or flagged in exclude are left out of the ranking (they do not count toward k) and come back as NA.

Usage

percentile_rank(x, exclude = NULL)

Arguments

Value

Numeric vector the same length as x with percentile ranks in (0, 1) for the ranked measurements and NA elsewhere.

Examples

percentile_rank(c(0.2, 0.5, 0.5, 0.9))
percentile_rank(c(0.2, 0.5, 0.5, 0.995), exclude = c(0.2, 0.5, 0.5, 0.995) >= 0.99)

Compute and compare phonological contrast metrics

Description

phontrast() is the package's main entry point. It computes one or more category separation and overlap metrics for a two-category phonological contrast in a single call: Jensen-Shannon divergence and distance, the Pillai-Bartlett trace, Bhattacharyya distance and affinity, Mahalanobis distance, and proportional overlap. Choose the metrics you want with metrics; the default computes all of them. Results are returned globally or by group, in a wide format (one column per metric, the default) or a tidy long format (one row per metric per comparison). The percent_overlap values are 0–1 proportions, not 0–100 percentages.

Usage

phontrast(
  data,
  features,
  category_col,
  group_col = NULL,
  metrics = c("jsd", "js_distance", "pillai", "bhattacharyya", "mahalanobis", "overlap"),
  min_tokens = 20,
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  eval_on = c("pooled", "group1", "group2", "pooled_sample"),
  eval_n = NULL,
  eval_seed = NULL,
  engine = c("ks", "fast_diag", "fast_diagonal"),
  chunk_size = 1000L,
  eps = 1e-06,
  output = c("wide", "long"),
  do_boot = FALSE,
  n_boot = 1000,
  conf_level = 0.95,
  progress = TRUE,
  method = c("mc", "legacy"),
  density = c("kde", "mvnorm"),
  mc_n = 10000L,
  bw_scale = 1
)

Arguments

Details

Use estimate_jsd() when Jensen-Shannon divergence is the only outcome of interest, and the lower-level metric helpers when you need direct control over one estimator.

Metric directions differ. JSD, Jensen-Shannon distance, Pillai trace, Bhattacharyya distance, and Mahalanobis distance increase as categories become more separated. Percent overlap and Bhattacharyya affinity increase as categories overlap more. Long output includes orientation, separation_value, and separation_rank columns so all metrics can be read on a separation-oriented scale.

If do_boot = TRUE, each metric is recomputed on n_boot nonparametric bootstrap resamples to estimate uncertainty. This can take substantial time because every resample recomputes KDE, MANOVA, and covariance-based metrics. Progress messages are printed by default while bootstrapping is running; set progress = FALSE to suppress them.

Value

A data frame containing only the requested metrics. Wide output (the default) contains one column per requested metric plus pillai_p_value when Pillai is requested; with do_boot = TRUE it also includes metric-specific *_mean, *_sd, *_ci_lower, *_ci_upper, and *_n_boot columns. Long output contains metric, estimate, orientation, bounded_0_1, separation_value, separation_rank, and p_value (populated for the Pillai row, NA otherwise) columns; with do_boot = TRUE it also includes boot_mean, boot_sd, ci_lower, ci_upper, n_boot, and conf_level. The result carries class "phontrast_contrast", so plot() and ggplot2::autoplot() draw it directly via plot_overlap_metrics().

Examples

set.seed(2026)
vowels <- data.frame(
  speaker = rep(c("s01", "s02"), each = 60),
  vowel = rep(rep(c("ih", "eh"), each = 30), 2),
  f1 = c(
    rnorm(30, 500, 55), rnorm(30, 560, 60),
    rnorm(30, 510, 60), rnorm(30, 575, 65)
  ),
  f2 = c(
    rnorm(30, 1980, 150), rnorm(30, 1880, 155),
    rnorm(30, 1960, 160), rnorm(30, 1840, 165)
  )
)

# All metrics in one wide comparison table (the default), by speaker.
phontrast(
  data = vowels,
  features = c("f1", "f2"),
  category_col = "vowel",
  group_col = "speaker"
)

# A single metric in wide format.
phontrast(
  data = vowels,
  features = c("f1", "f2"),
  category_col = "vowel",
  group_col = "speaker",
  metrics = "pillai",
  output = "wide"
)

# Bootstrapping is useful but slower because every requested metric is
# recomputed on every resample. Use a larger n_boot for real analyses.
phontrast(
  data = vowels,
  features = "f1",
  category_col = "vowel",
  group_col = "speaker",
  metrics = c("jsd", "pillai"),
  do_boot = TRUE,
  n_boot = 5,
  progress = FALSE
)

Okabe-Ito color palette used by phontrast plots

Description

Returns n colors from the Okabe-Ito palette, a qualitative palette designed to be distinguishable under the common forms of color-vision deficiency. The palette is reordered so the first two colors (blue and vermillion) form the highest-contrast pair for two-category contrasts. For more than eight groups the palette is interpolated, with a warning, since interpolated qualitative colors lose their guarantees.

Usage

phontrast_palette(n = NULL)

Arguments

Value

A character vector of n hex colors, named for n <= 8.

Examples

phontrast_palette()
phontrast_palette(2)

Pillai trace for multivariate overlap

Description

Computes the Pillai-Bartlett trace from a MANOVA of features ~ category. Optionally adds separation estimates and a proportion-standardized Pillai score for an exactly two-category, no-covariate design.

Usage

pillai_overlap(data, features, category_col, proportion_standardized = FALSE)

Arguments

Details

With proportion_standardized = FALSE, the function returns the ordinary Pillai trace and its p-value exactly as before. Raw Pillai describes the categories in the realized data; the proportion-standardized score targets the balanced-design score for the same estimated underlying separation.

For two categories with realized counts n_1 and n_2, total N, p features, error degrees of freedom \nu_e=N-2, and H=2n_1n_2/N, the optional estimator chain is

\hat D^2 = 2\nu_e[V/(1-V)]/H,

\tilde D^2 = [(\nu_e-p-1)/\nu_e]\hat D^2 - 2p/H,

followed, when \tilde D^2 >= 0, by

V_{eq}=\tilde D^2/(4+\tilde D^2).

The unbiased estimator is as given by Lachenbruch and Mickey (1968), and V_{eq} implements Becker's (1986) correction and his two-group multivariate generalization.

If \tilde D^2 < 0, pillai_eq is NA, pillai_eq_fallback is TRUE, and d2_fallback contains only the multiplicatively corrected first term. This common near-merger outcome is labelled separately because the fallback retains the split-dependent bias bias_2p_over_H and is not a balanced-design equivalent.

The estimator chain assumes multivariate normality within each category and a common within-category covariance. It supports no covariates. Both categories must contain at least two complete tokens, and the within-class error SSCP must be nonsingular. When \nu_e-p-1 <= 0, all optional fields are returned as typed NA values with a warning. A minority category with fewer than p+1 tokens remains computable but is flagged by fragile_minority.

Because x/(4+x) is strictly concave, pillai_eq is slightly downward biased for the balanced-design target. This bias grows with the sampling variance of \tilde D^2 and is largest for small, imbalanced, near-merged samples.

Value

A list with elements pillai and p_value. When proportion_standardized = TRUE, the list additionally contains n1, n2, H, d2_plugin, d2_unbiased, pillai_eq, pillai_eq_fallback, d2_fallback, bias_2p_over_H, and fragile_minority. Counts follow the category factor-level order used by the model. On the fallback path, pillai_eq is NA and d2_fallback is labelled separately; off that path, d2_fallback is NA.

References

Becker, G. (1986). Correcting the point-biserial correlation for attenuation owing to unequal sample size. Journal of Experimental Education, 55(1), 5-8.

Lachenbruch, P. A., & Mickey, M. R. (1968). Estimation of error rates in discriminant analysis. Technometrics, 10(1), 1-11.

Berry, G. M. (2026). Beyond the null: Calibration, balance, and the interpretation of Pillai scores. Preprint.

Plot a phontrast() result directly

Description

phontrast() results carry the class "phontrast_contrast", so they can be plotted without an explicit call to plot_overlap_metrics(): plot(phontrast(...)) draws the metric comparison, and ggplot2::autoplot(phontrast(...)) returns the same plot unprinted for further composition.

Usage

## S3 method for class 'phontrast_contrast'
plot(x, ...)

## S3 method for class 'phontrast_contrast'
autoplot(object, ...)

Arguments

Value

plot() draws the plot and returns it invisibly; autoplot() returns the ggplot2 object unprinted.

Examples

set.seed(2026)
vowels <- data.frame(
  vowel = rep(c("ih", "eh"), each = 40),
  f1 = c(rnorm(40, 500, 55), rnorm(40, 565, 60)),
  f2 = c(rnorm(40, 1980, 150), rnorm(40, 1870, 155))
)
if (requireNamespace("ggplot2", quietly = TRUE)) {
  plot(phontrast(vowels, c("f1", "f2"), "vowel", output = "long"))
}

Plot a rank_contrasts() result directly

Description

rank_contrasts() results carry the class "phontrast_ranking", so plot(ranking) draws the rank-agreement plot of plot_rank_agreement() and ggplot2::autoplot(ranking) returns it unprinted.

Usage

## S3 method for class 'phontrast_ranking'
plot(x, ...)

## S3 method for class 'phontrast_ranking'
autoplot(object, ...)

Arguments

Value

plot() draws the plot and returns it invisibly; autoplot() returns the ggplot2 object unprinted.

Plot a PCA projection of multidimensional category space

Description

Projects an arbitrary multidimensional acoustic feature set onto two principal components and visualizes the result with ggplot2. This is a diagnostic plot for high-dimensional workflows: metric estimates should still be computed on the full feature set when that is the intended analysis.

Usage

plot_category_pca(
  data,
  features,
  category_col,
  group_col = NULL,
  components = c(1L, 2L),
  center = TRUE,
  scale. = TRUE,
  points = TRUE,
  ellipses = TRUE,
  point_alpha = 0.65,
  point_size = 1.8,
  equal_axes = TRUE,
  facet_scales = c("fixed", "free", "free_x", "free_y")
)

Arguments

Value

A ggplot2 plot object. The fitted prcomp object and variance-explained table are stored as "pca" and "variance_explained" attributes.

Examples

set.seed(2026)
features <- paste0("mfcc", 1:5)
vowels <- data.frame(
  vowel = rep(c("ih", "eh"), each = 50),
  matrix(rnorm(100 * length(features)), ncol = length(features))
)
names(vowels)[-1] <- features
vowels[vowels$vowel == "eh", features[1:2]] <- vowels[vowels$vowel == "eh", features[1:2]] + 0.8

if (requireNamespace("ggplot2", quietly = TRUE)) {
  plot_category_pca(vowels, features = features, category_col = "vowel")
}

Plot phonological categories in acoustic space

Description

Creates a ggplot2 visualization of one or two acoustic dimensions from a token-level table. One-dimensional inputs are plotted as density curves; two-dimensional inputs are plotted as category-colored scatterplots with optional normal ellipses.

Usage

plot_category_space(
  data,
  features,
  category_col,
  group_col = NULL,
  points = TRUE,
  ellipses = TRUE,
  point_alpha = 0.65,
  point_size = 1.8,
  reverse_x = FALSE,
  reverse_y = FALSE,
  equal_axes = FALSE,
  facet_scales = c("fixed", "free", "free_x", "free_y")
)

Arguments

Value

A ggplot2 plot object.

Examples

set.seed(2026)
vowels <- data.frame(
  vowel = rep(c("ih", "eh"), each = 40),
  f1 = c(rnorm(40, 500, 55), rnorm(40, 565, 60)),
  f2 = c(rnorm(40, 1980, 150), rnorm(40, 1870, 155))
)

if (requireNamespace("ggplot2", quietly = TRUE)) {
  plot_category_space(vowels, features = "f1", category_col = "vowel")
  plot_category_space(
    vowels,
    features = c("f2", "f1"),
    category_col = "vowel",
    reverse_x = TRUE,
    reverse_y = TRUE
  )
}

Distribution-aware contrast plot for two categories

Description

The flagship visualization for a two-category contrast. Unlike a decorative scatterplot, plot_contrast() draws the same density model the distributional metrics are computed from: under density = "kde" it shows highest-density regions of each category's kernel density estimate (same bandwidth selection as jsd_kde_nd()); under density = "mvnorm" it shows coverage ellipses of the fitted multivariate normals used by the parametric backend. The pointwise minimum of the two densities – the mass that the proportional-overlap metric integrates – is shaded, so the overlap itself is visible rather than implied.

Usage

plot_contrast(
  data,
  features,
  category_col,
  group_col = NULL,
  density = c("kde", "mvnorm"),
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  bw_scale = 1,
  levels = c(0.5, 0.8, 0.95),
  points = TRUE,
  overlap = TRUE,
  annotate = TRUE,
  n_boot = 0,
  conf_level = 0.95,
  min_tokens = 20,
  mc_n = 10000L,
  eval_seed = NULL,
  grid_n = NULL,
  point_alpha = 0.55,
  point_size = 1.6,
  reverse_x = FALSE,
  reverse_y = FALSE,
  facet_scales = c("fixed", "free", "free_x", "free_y")
)

Arguments

Details

With annotate = TRUE (default) the panel is labelled with the Jensen-Shannon divergence and proportional overlap computed by phontrast() under the same density, bw, mc_n, and eval_seed settings, and the caption records the estimator configuration. The full annotation table is attached to the returned plot as attr(p, "contrast_metrics").

Value

A ggplot2 plot object. When annotate = TRUE, the phontrast() table behind the labels is attached as attr(p, "contrast_metrics").

Examples

set.seed(2026)
vowels <- data.frame(
  vowel = rep(c("ih", "eh"), each = 60),
  f1 = c(rnorm(60, 500, 55), rnorm(60, 565, 60)),
  f2 = c(rnorm(60, 1980, 150), rnorm(60, 1870, 155))
)

if (requireNamespace("ggplot2", quietly = TRUE)) {
  # Two-feature contrast in vowel-space orientation, KDE regions.
  plot_contrast(vowels, c("f2", "f1"), "vowel",
                reverse_x = TRUE, reverse_y = TRUE)

  # One-feature contrast with the overlap ribbon.
  plot_contrast(vowels, "f1", "vowel")

  # The same contrast under the multivariate-normal backend.
  plot_contrast(vowels, c("f2", "f1"), "vowel", density = "mvnorm",
                eval_seed = 2026, reverse_x = TRUE, reverse_y = TRUE)
}

Plot overlap metric comparisons

Description

Visualizes the output of phontrast() with ggplot2. The input may be either wide or long output. By default, values are plotted on a separation-oriented scale so overlap-oriented metrics are transformed as 1 - estimate.

Usage

plot_overlap_metrics(
  metrics,
  value = c("separation", "estimate"),
  metric = NULL,
  group_col = NULL,
  show_ci = TRUE,
  facet = TRUE,
  sort = TRUE
)

Arguments

Value

A ggplot2 plot object.

Examples

set.seed(2026)
vowels <- data.frame(
  speaker = rep(c("s01", "s02"), each = 60),
  vowel = rep(rep(c("ih", "eh"), each = 30), 2),
  f1 = c(rnorm(30, 500, 55), rnorm(30, 560, 60),
         rnorm(30, 510, 60), rnorm(30, 575, 65)),
  f2 = c(rnorm(30, 1980, 150), rnorm(30, 1880, 155),
         rnorm(30, 1960, 160), rnorm(30, 1840, 165))
)

metrics <- phontrast(
  vowels,
  features = c("f1", "f2"),
  category_col = "vowel",
  group_col = "speaker",
  output = "long"
)

if (requireNamespace("ggplot2", quietly = TRUE)) {
  plot_overlap_metrics(metrics)
}

Rank-agreement plot for a speaker ranking

Description

Draws step 3 of the ranking protocol behind rank_contrasts(): each speaker's Pillai percentile rank against its Jensen-Shannon distance percentile rank, the identity line on which the two measures agree, and the inspection band of margin either side of it. Speakers outside the band are flagged; speakers the bandwidth check set aside are crossed; speakers whose flag may not be read at this sample size and dimensionality are hollow. Speakers at the \sqrt{JSD} ceiling have no rank and are listed in the caption instead of drawn.

Usage

plot_rank_agreement(
  ranking,
  label = c("flagged", "all", "none"),
  base_size = 12
)

Arguments

Value

A ggplot2 plot object.

See Also

rank_contrasts(), inspect_contrast().

Examples

set.seed(2026)
gaps <- c(30, 60, 90, 120, 150, 180, 210, 240)
cohort <- do.call(rbind, lapply(seq_along(gaps), function(i) {
  data.frame(
    speaker = sprintf("s%02d", i),
    vowel = rep(c("ih", "eh"), each = 50),
    f1 = c(rnorm(50, 500, 55), rnorm(50, 500 + gaps[i], 55)),
    f2 = c(rnorm(50, 1900, 120), rnorm(50, 1900 - gaps[i], 120))
  )
}))
ranking <- rank_contrasts(cohort, c("f1", "f2"), "vowel", "speaker",
                          bw_check = FALSE)
if (requireNamespace("ggplot2", quietly = TRUE)) {
  plot_rank_agreement(ranking, label = "all")
}

Prepare JSD estimates for beta regression / GAMs

Description

JSD lives in [0,1]. This helper adds a jsd_beta column bounded in (0,1) so it can be used with the Beta family (e.g., in mgcv).

Usage

prepare_jsd_beta(jsd_df, jsd_col = "jsd_mean", eps = 1e-06)

Arguments

Value

A modified data frame with an added jsd_beta column.

Examples

jsd_by_speaker <- data.frame(
  speaker = paste0("s", 1:8),
  age = c(18, 22, 27, 31, 38, 45, 52, 60),
  jsd_mean = c(0.02, 0.05, 0.08, 0.13, 0.18, 0.24, 0.31, 0.39)
)

model_data <- prepare_jsd_beta(jsd_by_speaker)
model_data

if (requireNamespace("mgcv", quietly = TRUE)) {
  fit <- mgcv::gam(
    jsd_beta ~ age,
    data = model_data,
    family = mgcv::betar(),
    method = "REML"
  )
  stats::predict(fit, type = "response")
}

Sample-size floors of the ranking protocol

Description

The token floors per category at which the simulation study behind rank_contrasts() licenses a \sqrt{JSD} ranking and the Pillai–\sqrt{JSD} agreement flag, by number of features d (Berry, under review, Secs. V.D and VII.A). Dimensionalities the study did not calibrate take the floors of the next higher calibrated dimensionality.

Usage

protocol_floors(d)

Arguments

Value

A list with d; rank_floor, the tokens per category from which ranking by \sqrt{JSD} is licensed (50 at d <= 2, 200 at d = 3, 4, 500 at d = 5 to 8, Inf above eight dimensions, where the study gives ordering evidence only); flag_floor, the tokens per category from which the agreement flag may be read (100 at d <= 2, 200 at d = 3, 4, Inf from d = 5, since from eight dimensions the test returns agreement whatever the data contain); calibrated_at, the calibrated dimensionality the floors come from; and ordering_only, TRUE above eight dimensions.

Examples

protocol_floors(2)
protocol_floors(4)$flag_floor

Rank speakers' contrasts by Jensen-Shannon distance and check them against Pillai

Description

One-call implementation of the measurement protocol recommended by the simulation study behind phontrast (Berry, under review, Sec. VII.A). For each speaker (or other unit named by group_col) it

1. computes Jensen-Shannon distance (\sqrt{JSD}) and the Pillai trace on the same tokens, from one shared kernel density estimate, and reports the estimated shared probability mass beside them;

2. ranks the speakers by \sqrt{JSD} on the percentile-rank scale of percentile_rank();

3. flags speakers whose Pillai percentile rank differs from their \sqrt{JSD} percentile rank by margin (0.25 of the ordering) or more, for inspection with inspect_contrast() or plot_contrast().

The conditions the study attaches to those steps are applied and reported rather than left to the user: 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.

Usage

rank_contrasts(
  data,
  features,
  category_col,
  group_col,
  margin = 0.25,
  ceiling = 0.99,
  bw_check = TRUE,
  estimator = recommended_estimator(length(features)),
  min_tokens = 10,
  eval_seed = NULL,
  chunk_size = 1000L
)

Arguments

Value

A tibble of class "phontrast_ranking", one row per measured speaker, sorted by sqrt_jsd (largest first), with columns group, n_tokens, n_min, sqrt_jsd, pillai, shared_mass, at_ceiling, pr_jsd, pr_pillai, rank_diff (pr_jsd - pr_pillai), flag, rank_licensed, flag_licensed, rank_basis, and, with bw_check = TRUE, sqrt_jsd_half, sqrt_jsd_double, bw_shift, sign_change, and set_aside. The protocol settings, the floors, and the cleaned tokens the ranking was computed from are attached as attr(x, "protocol") (so inspect_contrast() can redraw a speaker without the original data); print() summarizes them. plot() and ggplot2::autoplot() draw the rank-agreement plot.

Ceiling

A measurement with \sqrt{JSD} \ge ceiling (default 0.99) is at the ceiling: the kernel estimate can no longer separate degrees of separation. Its shared mass and Pillai value are reported, but it is excluded from both percentile rankings (at_ceiling = TRUE, pr_jsd and pr_pillai are NA), and the ranks of the remaining k speakers are computed over those k alone.

Sample-size floors

Floors are read from the smaller category's token count per speaker (n_min) and from the number of features d (see protocol_floors()). Ranking by \sqrt{JSD} is licensed from 50 tokens per category at d = 2, 200 at d = 3 or 4, and 500 at d = 5 to 8; above eight dimensions the study gives ordering evidence only, so no rank is licensed. Below the floor rank_basis is "pillai": order that speaker by Pillai. The flag has floors of its own: 100 tokens per category at d = 2 and 200 at d = 3 or 4; from eight dimensions up the agreement test returns agreement whatever the data contain, and it is undefined when d >= 2 * n_min. Where the flag may not be read, flag is NA and rank_diff is still reported. The floors are where the simulation recovers the average ordering of its separation levels; they do not guarantee an accurate ranking of an individual speaker. Eight speakers is the smallest useful set (the calibration used sets of 32); fewer draws a warning.

Bandwidth check

With bw_check = TRUE (the default) \sqrt{JSD} is recomputed at half and at twice the diagonal Scott bandwidth (bw = "scott.diag", bw_scale = 0.5 and 2), at two dimensions as well, where the reported estimate uses the plug-in rule. The two re-estimates are ranked over the same speakers and bw_shift is the percentile-rank change between them. A measurement is set aside (set_aside = TRUE) when abs(bw_shift) >= margin, or when a flagged Pillai–\sqrt{JSD} rank difference changes sign between the halved and the doubled bandwidth (sign_change). The check tests sensitivity to smoothing, not the estimate's accuracy.

See Also

percentile_rank(), protocol_floors(), recommended_estimator(), inspect_contrast(), plot_contrast(), phontrast().

Examples

set.seed(2026)
gaps <- c(30, 60, 90, 120, 150, 180, 210, 240)
cohort <- do.call(rbind, lapply(seq_along(gaps), function(i) {
  data.frame(
    speaker = sprintf("s%02d", i),
    vowel = rep(c("ih", "eh"), each = 50),
    f1 = c(rnorm(50, 500, 55), rnorm(50, 500 + gaps[i], 55)),
    f2 = c(rnorm(50, 1900, 120), rnorm(50, 1900 - gaps[i], 120))
  )
}))
ranking <- rank_contrasts(cohort, c("f1", "f2"), "vowel", "speaker",
                          bw_check = FALSE)
ranking

Kernel estimator settings by dimensionality

Description

Returns the kernel-density estimator configuration that the simulation study behind rank_contrasts() used at each dimensionality (Berry, under review, Table III), so that a contrast measured with phontrast is scored with the settings the safe-use envelope was calibrated on. Three settings switch together with the number of features d:

Usage

recommended_estimator(d)

Arguments

Details

d <= 4

Plug-in bandwidth (bw = "Hpi"), the "ks" engine, densities evaluated at every token, and the partial leave-one-out correction on. Calibrated at d = 2 and 4.

5 <= d <= 13

Diagonal Scott bandwidth (bw = "scott.diag"), the "fast_diag" engine, densities evaluated at 200 subsampled tokens per category, leave-one-out on. Calibrated at d = 8 and 13; dimensionalities between the calibrated ones take the settings of the next higher calibrated dimensionality.

d >= 14

As above but with no leave-one-out correction, reproducing the log-space path the study used at d = 32 and 64. phontrast's "fast_diag" engine evaluates kernels in log space, so this tier runs natively. Above eight dimensions the study gives ordering evidence only.

Value

A list with elements d, tier (a label for the row of Table III applied), calibrated_at (the dimensionalities the row was calibrated on), bw, engine, eval_n, loo, and note. The bw, engine, eval_n, and loo elements can be passed straight to jsd_kde_nd(), estimate_jsd(), or phontrast().

Examples

recommended_estimator(2)
recommended_estimator(13)$engine

set.seed(2026)
vowels <- data.frame(
  vowel = rep(c("ih", "eh"), each = 40),
  f1 = c(rnorm(40, 500, 55), rnorm(40, 565, 60)),
  f2 = c(rnorm(40, 1980, 150), rnorm(40, 1870, 155))
)
est <- recommended_estimator(2)
jsd_kde_nd(vowels, c("f1", "f2"), "vowel",
           bw = est$bw, engine = est$engine, loo = est$loo)

Okabe-Ito discrete color and fill scales

Description

Discrete ggplot2 scales built on phontrast_palette(). These are the default scales for all phontrast plotting functions and can be added to any ggplot to match the package's visual identity.

Usage

scale_colour_phontrast(...)

scale_color_phontrast(...)

scale_fill_phontrast(...)

Arguments

Value

A ggplot2 scale object.

Examples

if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::ggplot(iris, ggplot2::aes(Sepal.Length, Sepal.Width,
                                     color = Species)) +
    ggplot2::geom_point() +
    scale_color_phontrast() +
    theme_phontrast()
}

Bhattacharyya distance by group

Description

Computes Bhattacharyya distance and affinity for each group, under a multivariate normal approximation.

Usage

speaker_bhatt(
  data,
  group_col,
  category_col,
  features,
  min_tokens = 20,
  eps = 1e-06
)

Arguments

Value

Data frame with columns: group, n_tokens, bhatt_dist, bhatt_affinity.

Group-level JSD point estimates

Description

Computes JSD for each group (e.g., speaker) comparing two categories (e.g., vowels) in an n-dimensional acoustic space.

Usage

speaker_jsd(
  data,
  group_col,
  category_col,
  features,
  min_tokens = 20,
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  eval_on = c("pooled", "group1", "group2", "pooled_sample"),
  eval_n = NULL,
  eval_seed = NULL,
  engine = c("ks", "fast_diag", "fast_diagonal"),
  chunk_size = 1000L,
  ...
)

Arguments

Value

A tibble with one row per group and columns: group, n_tokens, and jsd.

Group-level Pillai scores

Description

Computes Pillai scores and associated p-values per group (e.g., per speaker).

Usage

speaker_pillai(data, group_col, category_col, features, min_tokens = 20)

Arguments

Value

A tibble with columns: group, n_tokens, pillai, p_value.

Publication theme for phontrast plots

Description

A minimal, publication-oriented ggplot2 theme: quiet major grid, no minor grid, bold plot-aligned title, muted subtitle, left-aligned caption for estimator provenance, and a top-aligned legend. Applied by default in all phontrast plotting functions.

Usage

theme_phontrast(base_size = 12, base_family = "")

Arguments

Value

A ggplot2 theme object.

Examples

if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::ggplot(iris, ggplot2::aes(Sepal.Length, Sepal.Width)) +
    ggplot2::geom_point() +
    theme_phontrast()
}

Simulated twelve-speaker vowel cohort

Description

A simulated cohort for demonstrating the ranking protocol of rank_contrasts(): twelve speakers, two vowel categories ("ih" and "eh"), and two acoustic features (F1 and F2, in Hz). Each speaker was built to exercise one part of the protocol:

spk01–spk08

Gaussian categories whose centroid gap grows from 25 to 200 Hz, 100 tokens per vowel: a graded ordering that both Jensen-Shannon distance and Pillai recover.

spk09

The same centroids for both vowels, but a bimodal "eh" (two variants 190 Hz apart in F1 and 440 Hz apart in F2), 100 tokens per vowel. A mean-based measure sees no contrast while the distributions barely overlap: the planted Pillai / \sqrt{JSD} disagreement the agreement flag should catch.

spk10

Fully separated categories, 100 tokens per vowel: \sqrt{JSD} at the ceiling.

spk11

60 tokens per vowel: ranking by \sqrt{JSD} is licensed at two dimensions (floor 50) but the flag is not readable (floor 100).

spk12

40 tokens per vowel: below the two-dimensional rank floor, so the speaker is ordered by Pillai.

Usage

vowel_cohort

Format

A data frame with 2200 rows and 4 columns:

speaker

Character; speaker identifier "spk01" to "spk12".

vowel

Character; vowel category, "ih" or "eh".

f1

Numeric; first formant frequency in Hz.

f2

Numeric; second formant frequency in Hz.

Source

Simulated with a fixed seed; the generating script is data-raw/vowel_cohort.R in the package repository.

See Also

rank_contrasts(), inspect_contrast().

Examples

head(vowel_cohort)
table(vowel_cohort$speaker, vowel_cohort$vowel)