Rで解析：音韻カテゴリーの分離・重なりを複数指標で評価する「phontrast」パッケージの紹介

音声学や言語学の研究では、母音や子音といった音韻カテゴリーが、フォルマント値やMFCCなどの音響空間上でどれだけ明確に区別されているかを確認したい場面があります。しかし、複数の指標を一括で計算し、話者やグループごとのばらつきを考慮しながら評価するには手間がかかります。

「phontrast」パッケージは、多次元の音響空間における2つのカテゴリー間の分離度や重なり具合を、複数の指標を用いて一括で算出するパッケージです。Pillaiのトレースや Bhattacharyya 距離、Jensen-Shannon 情報量など複数の指標をまとめて計算でき、話者単位の階層ブートストラップによる年齢効果の検定も可能です。また、算出結果にもとづくカテゴリー空間の可視化や、主成分分析による高次元特徴量の投影図のプロットも可能です。本パッケージの利用で、音韻のコントラストを複数の指標から多角的に確認できるのではないかと考えます。

パッケージバージョンは2.4.0。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。

＜おすすめのRに関する書籍です＞

Rによる計量政治学 | 浅野正彦, 矢内勇生

政治現象を統計的に検証するための入門書。各章に練習問題があり、独学でも進めやすい構成。

パッケージのインストール

下記コマンドを実行してください。

# パッケージのインストール
install.packages("phontrast")
# パッケージの読み込み
library("phontrast")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

本パッケージの中心となるのはphontrastコマンドで、2つの音韻カテゴリーについてPillaiのトレースやBhattacharyya距離、Jensen-Shannon情報量など複数の分離・重なり指標を一括で算出します。話者ごとの経年変化を調べたい場合はhier_boot_jsd_modelコマンドで階層ブートストラップモデルを当てはめられます。また、plot_contrastコマンドやplot_category_spaceコマンド、plot_category_pcaコマンドでカテゴリー間の分布を可視化でき、plot_overlap_metricsコマンドで複数の比較対象の指標をまとめて確認できます。配色にはphontrast_paletteコマンドとscale_colour_phontrastコマンドが用意されています。

本記事では、京都市右京区と与謝郡伊根町の高校生が発音した母音「い」「え」のフォルマント測定値を架空のデータとして作成し、記事全体を通して使用します。

# 乱数シードの固定
set.seed(2026)

# 右京区と伊根町の話者情報を作成
chiku <- rep(c("右京区", "伊根町"), each = 4)

wakate <- data.frame(
  speaker = paste0("sp", sprintf("%02d", 1:8)),
  chiku = chiku,
  age = rep(c(16, 17, 18, 17), 2)
)

# 話者ごとに「い」「え」を25回ずつ発音したデータの骨組みを作成
boin <- data.frame(
  speaker = rep(wakate$speaker, each = 50),
  vowel = rep(rep(c("い", "え"), each = 25), nrow(wakate))
)

# 話者情報を結合
boin <- merge(boin, wakate, by = "speaker", sort = FALSE)

# 母音と地区の違いを反映した第1フォルマント（F1）を生成
boin$f1 <- rnorm(nrow(boin),
                 mean = ifelse(boin$vowel == "い", 320, 480) +
                   ifelse(boin$chiku == "伊根町", 25, 0),
                 sd = 45)

# 母音と地区の違いを反映した第2フォルマント（F2）を生成
boin$f2 <- rnorm(nrow(boin),
                 mean = ifelse(boin$vowel == "い", 2400, 2050) +
                   ifelse(boin$chiku == "伊根町", -110, 0),
                 sd = 140)

# 第3フォルマント（F3）を生成
boin$f3 <- rnorm(nrow(boin), mean = ifelse(boin$vowel == "い", 3000, 2900), sd = 150)
# スペクトル傾動を生成

boin$tilt <- rnorm(nrow(boin),
                   mean = ifelse(boin$vowel == "い", -8, -5) +
                     ifelse(boin$chiku == "伊根町", -1, 0),
                   sd = 2)

2つのカテゴリー間の分離・重なりを複数指標で算出：phontrastコマンド

2つの音韻カテゴリーに含まれる音響特徴量から、Pillaiのトレース、Bhattacharyya距離・近似度、Jensen-Shannon情報量とその平方根であるJensen-Shannon距離、Mahalanobis距離、比率重なりの各指標を一括で算出します。metricsで算出する指標を絞り込め、group_colを指定するとグループごとに算出できます。outputを”long”にすると指標ごとに1行の縦持ち形式になり、後述のplotコマンドやplot_overlap_metricsコマンドにそのまま渡せます。

# 右京区・伊根町の母音「い」「え」の分離度をまとめて算出（全列を表示）
print(phontrast(boin, features = c("f1", "f2"), category_col = "vowel"), width = Inf)
# A tibble: 1 × 10
  scope  n_tokens pillai pillai_p_value bhatt_dist bhatt_affinity   jsd
  <chr>     <int>  <dbl>          <dbl>      <dbl>          <dbl> <dbl>
1 global      400  0.814      8.44e-146       2.18          0.113 0.950
  js_distance mahalanobis_dist percent_overlap
        <dbl>            <dbl>           <dbl>
1       0.975             4.18          0.0238

話者ごとの階層ブートストラップでJensen-Shannon情報量の年齢依存性をモデル化：hier_boot_jsd_modelコマンド

話者をグループとしてリサンプリングしながら、話者ごとのJensen-Shannon情報量にformulaで指定したモデルを繰り返し当てはめます。戻り値は再標本化ごとの回帰係数の一覧で、係数のばらつきから年齢の効果の確からしさを確認できます。

# 乱数シードの固定
set.seed(2026)
# 話者ごとのJensen-Shannon情報量に年齢を説明変数とした回帰を階層ブートストラップで当てはめ
hier_boot_jsd_model(
  data = boin,
  group_col = "speaker",
  category_col = "vowel",
  features = "f1",
  formula = jsd_beta ~ age,
  fit_fun = stats::lm,
  n_outer = 5,
  min_tokens = 20,
  progress = FALSE
)
# A tibble: 10 × 3
   boot_id term        estimate
     <int> <chr>          <dbl>
 1       1 (Intercept)   0.584
 2       1 age           0.0226
 3       2 (Intercept)   1.81
 4       2 age          -0.0514
 5       3 (Intercept)   1.76
 6       3 age          -0.0506
 7       4 (Intercept)   1.13
 8       4 age          -0.0137
 9       5 (Intercept)   2.15
10       5 age          -0.0710

2つのカテゴリー間の分布を考慮した対比図のプロット：plot_contrastコマンド

featuresに2つの特徴量を指定すると、密度領域を重ねた特徴量空間の図になります。1つだけ指定した場合は密度曲線になります。overlapで2つのカテゴリー密度が重なる範囲に網掛けし、annotateでJensen-Shannon情報量と比率重なりをパネルに注釈します。

# ggplot2が利用可能であれば対比図を描画
if (requireNamespace("ggplot2", quietly = TRUE)) {
  # 母音空間の慣例に合わせてF2を横軸に反転した対比図
  plot_contrast(boin, c("f2", "f1"), "vowel", reverse_x = TRUE, reverse_y = TRUE)
  # F1のみを使った1次元の対比図
  plot_contrast(boin, "f1", "vowel")
}

＜おすすめのRに関する書籍です＞

東京大学のデータサイエンティスト育成講座 ~Pythonで手を動かして学ぶデ―タ分析~

データ分析の一連の流れを手を動かして学ぶ定番書。ただしコードはRではなくPython中心なので、その点は注意。

www.amazon.co.jp

phontrastの算出結果を直接プロット：plotコマンド

phontrastコマンドの戻り値（output = “long”）にはplotコマンドが用意されています。内部でplot_overlap_metricsコマンドを呼び出し、算出済みの指標をそのままプロットします。

# ggplot2が利用可能であればphontrastの縦持ち結果をそのままプロット
if (requireNamespace("ggplot2", quietly = TRUE)) {
  plot(phontrast(boin, c("f1", "f2"), "vowel", output = "long"))
}

高次元の音響特徴量をPCAで低次元に投影してプロット：plot_category_pcaコマンド

3個以上の音響特徴量を扱う場合に、主成分分析で2次元に投影してカテゴリーの分布を確認できます。componentsで表示する主成分を切り替えられます。

# ggplot2が利用可能であれば4つの音響特徴量をPCAで2次元に投影してプロット
if (requireNamespace("ggplot2", quietly = TRUE)) {
  plot_category_pca(boin, features = c("f1", "f2", "f3", "tilt"), category_col = "vowel")
}

音響空間における音韻カテゴリーのプロット：plot_category_spaceコマンド

featuresに1〜2個の特徴量を指定し、観測トークンと正規楕円で音響空間上のカテゴリー分布をそのまま描画します。plot_contrastコマンドと異なり、密度領域や重なりの網掛けは表示しません。

# ggplot2が利用可能であればカテゴリ空間をプロット
if (requireNamespace("ggplot2", quietly = TRUE)) {
  # F1のみを使った1次元のカテゴリ空間
  plot_category_space(boin, features = "f1", category_col = "vowel")
  # 母音空間の慣例に合わせてF2とF1で2次元のカテゴリ空間
  plot_category_space(boin, features = c("f2", "f1"), category_col = "vowel",
                      reverse_x = TRUE, reverse_y = TRUE)
}

＜おすすめのRに関する書籍です＞

東京大学のデータサイエンティスト育成講座 ~Pythonで手を動かして学ぶデ―タ分析~

データ分析の一連の流れを手を動かして学ぶ定番書。ただしコードはRではなくPython中心なので、その点は注意。

www.amazon.co.jp

地区間で分離・重なり指標を比較する棒グラフの作成：plot_overlap_metricsコマンド

phontrastコマンドの縦持ち出力（複数のグループを含むもの）を渡すと、指標ごとにグループ間の値を比較する棒グラフを描画します。show_ciで信頼区間を、sortで比較対象の並べ替えを制御します。

# 地区ごとの分離・重なり指標を縦持ち形式で算出
chiku_shihyo <- phontrast(boin, features = c("f1", "f2"), category_col = "vowel",
                          group_col = "chiku", output = "long")
# ggplot2が利用可能であれば地区間を比較する棒グラフを描画
if (requireNamespace("ggplot2", quietly = TRUE)) {
  plot_overlap_metrics(chiku_shihyo)
}

phontrastパッケージの標準カラーパレットを取得：phontrast_paletteコマンド

nを指定しない場合はパッケージ標準の8色を返します。色覚多様性に配慮した配色になっています。

# 標準の8色を取得
phontrast_palette()
      blue vermillion      green     orange     purple    skyblue     yellow
 "#0072B2"  "#D55E00"  "#009E73"  "#E69F00"  "#CC79A7"  "#56B4E9"  "#F0E442"
      grey
 "#999999"

# 先頭の2色のみを取得
phontrast_palette(2)
      blue vermillion
 "#0072B2"  "#D55E00"

phontrastの配色をggplot2に適用する色スケール：scale_colour_phontrastコマンド

ggplot2のプロットにphontrastパッケージの配色とテーマをまとめて適用します。scale_colour_phontrastコマンドは離散スケール用で、scale_color_phontrastコマンドは同じ関数の別名です。

# ggplot2が利用可能であればphontrastの配色とテーマを適用
if (requireNamespace("ggplot2", quietly = TRUE)) {
  # 地区ごとに色分け
  ggplot2::ggplot(boin, ggplot2::aes(f2, f1, color = chiku)) +
    ggplot2::geom_point() +
    scale_colour_phontrast() +
    theme_phontrast()
}

＜おすすめのRに関する書籍です＞

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装

時系列データの扱いを、理論とRおよびStanでの実装の両面から解説。状態空間モデルの入口として。

www.amazon.co.jp

この記事が誰かの役に立ちますように。