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.

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.

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
)

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,
  ...
)

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).
estimate_jsd(
  data = vowels,
  features = c("f1", "f2"),
  category_col = "vowel"
)

# Bootstrap confidence intervals, global and by speaker.
# Increase n_boot for real analyses.
estimate_jsd(vowels, c("f1", "f2"), "vowel",
             do_boot = TRUE, n_boot = 10)
estimate_jsd(vowels, c("f1", "f2"), "vowel", group_col = "speaker",
             do_boot = TRUE, n_boot = 10)

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,
  ...
)

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
)

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
)

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,
  ...
)

Arguments

Value

Numeric scalar proportion in [0, 1].

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
)

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. This is a convenience wrapper for comparison with JSD.

Usage

pillai_overlap(data, features, category_col)

Arguments

Value

A list with elements pillai and p_value.

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 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"),
  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)
}

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")
}

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()
}