Rで試す「新しい」相関尺度

作者

所属

鈴木 了太

株式会社 ef-prime

はじめに

本稿はR言語 Advent Calendar 2025の12日目の記事です。現在進行中の研究に関する内容を含み、カレンダー登録日（2025/12/12）以降も更新される可能性があります。

先日、「新しい相関係数」についてという記事が話題になりました。いわゆる「Chatterjee相関」を紹介する記事で、手法の概要がわかりやすくまとまっています。

偶然にもその翌日、Japan.R 2025にて「パッケージ環境の再現ツールと新しい相関指標」というタイトルで、筆者が開発している別の「新しい」相関尺度について発表させていただきました。

よい機会ですので、これらの手法をRで利用する方法をご紹介したいと思います。また具体例を用いて値を算出し、MICや距離相関、HSICといった関連指標との比較についても試みます。

⚠️本稿は進行中の研究に関する暫定的な内容を含みます。記載内容は現時点のもので、将来的に変更の可能性があることにご注意ください（例えば記載のコード例が、そのままでは動かなくなる可能性があります）。

本稿の構成

まず具体例として用いるサンプルデータを概観します。続いてChatterjee相関などを計算するためのRパッケージを紹介し、具体例を用いて算出方法と結果を確認します。最後にさまざまな関連指標を比較してまとめ、参考文献や公式サイトへのリンクを整理します。

サンプルデータ

サンプルデータとして、Wikipediaの「ピアソンの積率相関係数」に掲載されている例からいくつか借用したいと思います。

Wikimedia Commonsにおいてパブリックドメインとして公開されているコードをもとに、サンプルデータを生成する関数gen_example()および散布図を描画する関数plot_examples()を作成しました。コードはやや長くなるため、本稿の補遺セクションに掲載しておきます。

再現性を担保するため、乱数のシードをSEEDオブジェクトで指定します。plot_examples()関数で散布図を描画してみましょう。サンプルサイズはn = 1000に設定し、2次元正規分布を用いるケースでは相関係数をrho = 0.8に設定しました：

# 乱数シードを保持するオブジェクト
SEED <- 12345

# 散布図の描画
plot_examples(n = 1000, rho = 0.8, seed = SEED)

ここでtype: /は2変量正規分布に従うパターンで、rhoで指定した相関係数をもつ母集団からデータがサンプリングされます。それ以外のtype: U, W, X, O, 88は、元ネタとなったコードに含まれる特徴的なパターンについて、見た目に応じたラベルを付けたものです。

各タイプに対応するサンプルデータはgen_example()関数で得ることができます。ひとまずtype: Wのデータを作成しておきます：

# データの生成
d <- gen_example(type = "W", n = 1000, seed = SEED)

# 散布図の描画
plot(y ~ x, data = d, col = "darkblue", pch = 16, cex = 0.4)

古典的な積率相関係数について、標本相関係数を算出してみましょう：

# 標本相関係数（相関行列）
cor(d)

            x           y
x  1.00000000 -0.03175413
y -0.03175413  1.00000000

散布図からは強い非線形関係が確認できますが、xとyの標本相関係数はほぼ0に近い値となります。

XICORパッケージ

RパッケージXICORはChatterjee相関の実装を提供します。Chatterjee先生の論文「A new coefficient of correlation」やスライド資料においてギリシア文字\(\xi\)が用いられていることから、このような名称になったものと想像されます。

XICORパッケージはCRANに登録されており、install.packages()関数でインストールすることができます。値の算出にはxicor()関数を用います。さっそく試してみましょう：

# install.packages("XICOR")
library(XICOR)

# xとyの相関を評価
xicor(d$x, d$y)

[1] 0.5213735

Chatterjee相関の理論値は0から1の値をとり、このデータにおいては非線形関係を評価できているようです（ただし推定値は負の値をとる場合もあります）。pvalue = TRUEを指定すれば、独立性の仮説検定を行うことも可能です：

# p値の算出
xicor(d$x, d$y, pvalue = TRUE)

$xi
[1] 0.5213735

$sd
[1] 0.02000004

$pval
[1] 0

推定値（xi）に加えて標準誤差（sd）および\(p\)値（pval）が出力されました。この例では独立性の帰無仮説は棄却されるようです。

引数としてデータフレームを与えると、すべての変数の組み合わせについて値を算出し、行列形式で返すことができます：

# すべての組み合わせを算出し、行列形式で返す
xicor(d)

           x         y
x 0.99700300 0.5213735
y 0.09165609 0.9970030

この結果から、ふたつの興味深いことがわかります：

• 自分自身との「相関」がぴったり1になるとは限らない
• 対称性をもたず、xicor(x,y)とxicor(y,x)は異なる値を取り得る

Chatterjee相関は、母集団における理論値について「\(Y = f(X)\)という関数関係が成り立つとき（またそのときに限り）1」という性質を持つことが示されています。したがって前者については、「理論値は1であることがいえるものの、有限サンプルによる推定値が1となる保証はない」ということのようです。

Chatterjee相関の推定バイアスについては、こちらの記事「新しい相関係数の理論的性質」でも触れられています。

後者については、\(Y = f(X)\)という関数関係が成り立つ場合でも、\(X = g(Y)\)という関係が成り立つとは限らない、と考えれば納得できると思います。実際このtype: Wのデータでは、散布図をみてもxからyを予測することはうまくできそうですが、逆にyからxを予測することは難しいことが示唆されます。

以上のことからChatterjee相関は非線形な関係を捉えることができ、また「予測の方向」を加味した評価を与える指標と考えることができます。そのうえ独立性の仮説検定を行うことができ、計算も速い（計算量のオーダーは\(O(n \log n)\)）ということで、非常に強力かつ実用的な手法といえそうです。

gcorパッケージ

Rパッケージgcorは、筆者が提案する分割相関（contingency correlation）の実装を提供します。

ご利用にあたっては、本パッケージは開発版であることにご注意ください。本稿の執筆時点（2025年12月）では、手法自体の開発が進行中の段階にあります。ただし計算実験において一定の性能を確認できているため、早期に開発版を公開しています。

分割相関について

分割相関は数値だけでなく、カテゴリ値にも適用できる統一的な相関尺度（の推定値）として設計されています。古典的な相関係数を情報理論的アプローチで一般化したもので、非線形な関係についても評価することができます。ここまでの議論を踏まえて整理すると、次のような特徴をもちます：

1. 独立のとき0に近く、関数関係\(Y = f(X)\)が成り立つとき1に近い
2. \((X,Y)\)が2変量正規分布に従うとき、相関係数の絶対値\(|\rho|\)の近似を与える
3. 変数の単調変換に対して不変

1についてはChatterjee相関と類似の性質ですが、理論上「1のとき関数関係\(Y = f(X)\)が成り立つ」という逆の関係は成り立ちません。分割相関は「独立性からの乖離」を測る手法であり、決定論的でなくても、強い関係があれば1に近い値をとります。本稿では理論面には深入りしませんが、のちほど数値例で確認したいと思います。

2についてはあくまで近似的ですが、慣れ親しんだ相関係数と同じ感覚で値を評価することができます。また3の特徴と組み合わせると、「もし適当な単調変換\(u, v\)によって\(\big(u(X), v(Y)\big)\)が2変量正規分布に従うならば、その相関係数の絶対値を近似する」といえます。

利用例

では実際のデータに適用してみましょう。gcorパッケージの開発版はGitHubで公開されており、pakパッケージを用いて

# install.packages("pak")
pak::pak("r-suzuki/gcor")

とするか、またはdevtoolsパッケージを用いてインストールすることができます：

# install.packages("devtools")
devtools::install_github("r-suzuki/gcor")

値の算出にはgcor()関数を用います。さっそく試してみましょう：

library(gcor)

# xとyの相関を評価
gcor(d$x, d$y)

[1] 0.613481

xicor()と同様、こちらも非線形関係を評価できているようです。現状では仮説検定には対応しておらず、今後の課題として検討を進めています。

引数としてデータフレームを与えると、すべての変数の組み合わせについて値を算出し、行列形式で返すことができます：

# すべての組み合わせを算出し、行列形式で返す
gcor(d)

         x        y
x 1.000000 0.613481
y 0.613481 1.000000

こちらはxicor()と異なり、結果はcor()と同様に対称的となります。実は分割相関にも「予測の方向」を加味したバージョンがあるのですが、現状この例では同じ結果を返すこともあり、本稿では省略します。

計算時間

計算量のオーダーについては評価を進めているところですが、理論上はChatterjee相関と同等の\(O(n \log n)\)という評価になるものと考えています。実際の計算時間については、簡易的には次のように比較することができます：

system.time(for(i in 1:1000) xicor(d$x, d$y))

   user  system elapsed 
  0.246   0.024   0.273 

system.time(for(i in 1:1000) gcor(d$x, d$y))

   user  system elapsed 
  0.982   0.013   1.016 

この例ではxicor()のほうが高速ですが、行列形式の出力を得る場合は同等または逆転します：

system.time(for(i in 1:1000) xicor(d))

   user  system elapsed 
  1.478   0.070   1.657 

system.time(for(i in 1:1000) gcor(d))

   user  system elapsed 
  1.131   0.021   1.153 

gcor()では変数の組み合わせごとに1回の計算で済むところ、xicor()では変数の並び（xicor(x,y)とxicor(y,x)で異なる）ごとに計算が必要となります。またgcor()では「自分自身との相関は常に1」といった性質があり、計算対象となる組み合わせそのものが少なくなることが理由として考えられます。

離散データへの適用

前述のとおり、分割相関はカテゴリ変数にも適用することができます。例としてirisデータに対して、「Sepal.LengthとSpeciesの相関」などを評価することができます：

gcor(iris)

             Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
Sepal.Length    1.0000000   0.2349075    0.8846517   0.8741873 0.7623968
Sepal.Width     0.2349075   1.0000000    0.3143301   0.2669031 0.6510740
Petal.Length    0.8846517   0.3143301    1.0000000   0.9503289 0.8221674
Petal.Width     0.8741873   0.2669031    0.9503289   1.0000000 0.8237429
Species         0.7623968   0.6510740    0.8221674   0.8237429 1.0000000

この結果をみると、Speciesとの相関がもっとも強いのはPetal.Width、もっとも弱いのはSepal.Widthということになります。ボックスプロットを描画して確認してみましょう：

boxplot(Petal.Width ~ Species, data = iris)

boxplot(Sepal.Width ~ Species, data = iris)

確かにそのような関係が窺われます。このように、分割相関は数値とカテゴリが混在したデータを自然に扱えるように設計されています。

xicor()の場合も、カテゴリ変数を整数に変換することで（少なくとも形式的には）相関を評価することができます：

library(dplyr)
iris %>%
    # Speciesを「形式的に」整数に変換
    mutate(Species = as.integer(Species)) %>%
    xicor()

             Sepal.Length Sepal.Width Petal.Length Petal.Width Species
Sepal.Length    0.9799826   0.1651616    0.5515098   0.4910766  0.5650
Sepal.Width     0.1580596   0.9798165    0.2218503   0.1994221  0.1975
Petal.Length    0.5856799   0.3119999    0.9796016   0.7021707  0.8950
Petal.Width     0.5064163   0.3436975    0.6606205   0.9791778  0.8950
Species         0.4259436   0.2662145    0.6492576   0.6797689  0.9850

factor型における内部表現としての整数が、カテゴリの自然な順序を表す場合（すなわち順序尺度として解釈できる場合）、この方法は有効に働くものと期待されます。ただしタイが多数発生することから乱数を用いた処理が適用されるとみられ、（乱数シードを設定しない限り）計算のたびに結果が異なる様子もみられました。xicor()をカテゴリ値に対して適用する場合、これらの点に留意することが望ましいと思われます。

算出結果の比較

それでは他のサンプルデータについても値を算出し、比較してみましょう。いずれも乱数シードを設定して同じデータを生成しています。算出結果は小数点以下2桁で丸め、値が\(2/3\)以上のとき赤、\(1/3\)未満のとき緑、それ以外は黒で表示しています。

標本相関係数

コード例

cor(x, y)

相関係数\(\rho = 0.8\)の2次元正規分布であるtype: /については高い値をとり、真の\(\rho\)を精度よく推定できています。それ以外のパターンは0に近い値が算出されました。

plot_examples(title = "cor", n_digit = 2, seed = SEED)

Chatterjee相関

コード例

XICOR::xicor(x, y)

type: /, U, Wについては比較的高い値を示しました。関数関係とはいえないX, Oは想定通り低く、88についてはほぼ0に近い値が得られています。実際88は\(X\)と\(Y\)が独立なケースであるため、これについては望ましい結果といえます。

plot_examples(title = "xicor", n_digit = 2, seed = SEED)

分割相関

コード例

gcor::gcor(x, y)

type: 88を除くすべてのケースで比較的高い値が得られました。独立性からの乖離という観点からは、X, Oといったパターンも高く評価されることがわかります。

/については理論上、相関係数の絶対値である\(|\rho| = 0.8\)の近似となっており、近似誤差はあるものの比較的近い値が得られています。88については独立であるため理論値は0であり、こちらも近似誤差はありますが、望ましい傾向といえそうです。

plot_examples(title = "gcor", n_digit = 2, seed = SEED)

MIC

コード例

minerva::mine(x, y)$MIC

ここまで本稿では触れてきませんでしたが、同様に非線形相関を検出可能な手法であるMIC (maximal information coefficient)についても比較してみましょう。type: Wが特に高く評価され、/は原著論文にもある通り\(\rho^2 = 0.64\)に比較的近い値をとります。U, X, Oはそこそこの評価で、88はやや上振れしていますが小さめの値が得られています。

plot_examples(title = "mic", n_digit = 2, seed = SEED)

距離相関

コード例

energy::dcor(x, y)

同様に距離相関 (distance correlation)についても比較してみます。「\(X\)軸において近いとき、\(Y\)軸でも近い」といった観点から評価するため、Chatterjee相関と同様に関数的な関係を高く評価することが期待されます。

ここではtype: /がもっとも高い値をとり、\(\rho = 0.8\)に近い値をとります。U, Wのような「関数的」なものが比較的上位にあがるものの/よりは低く、X, Oがそれに続きます。独立である88についてはほぼ0となっており、望ましい結果といえます。

plot_examples(title = "dcor", n_digit = 2, seed = SEED)

HSICと独立性検定

HSIC

コード例

dHSIC::dhsic(x, y)$dHSIC

HSIC (Hilbert-Schmidt independence criterion)についても試してみましょう。ほかの指標と異なり「0から1の範囲にスケーリングされた指標」ではないものの、カーネル法を用いた独立性検定に関係する指標として利用することができます。

まずはHSICを算出してみます。既定では正規カーネルが設定され、バンド幅は自動選択されます。ほかの指標より小さな値をとるので、文字の色分けは行わず、小数点以下3桁で丸めて表示します：

plot_examples(title = "hsic", n_digit = 3, seed = SEED)

HSICによる独立性検定

コード例

dHSIC::dhsic.test(x, y)$p.value

HSICに基づく解釈しやすい値として、独立性検定の\(p\)値を算出してみましょう。値が0.01未満のとき赤、0.05未満のとき黒、0.05以上のとき緑とし、小数点以下2桁で丸めて表示します：

plot_examples(title = "hsic_test", n_digit = 2, seed = SEED)

独立であるtype: 88のパターンを除き、見事に独立性の帰無仮説が棄却されました。

距離相関による独立性検定

コード例

energy::dcor.test(x, y, R = 1000)$p.value

話題が前後しますが、距離相関についてもdcor.test()関数を用いて独立性検定の\(p\)値を算出することができます：

plot_examples(title = "dcor_test", n_digit = 2, seed = SEED)

独立である88パターンの値はHSICのほうが高く出る傾向が見られます。ただし、すでに見たように距離相関は「0から1の値をとる、相関係数のような指標」として解釈しやすいという利点があり、目的に応じた使い分けが有効と考えられます。

Chatterjee相関による独立性検定

コード例

XICOR::xicor(x, y, pvalue = TRUE)$pval

Chatterjee相関についても、前述のとおりxicor()関数にpvalue = TRUEを指定することで独立性検定の\(p\)値が得られます：

plot_examples(title = "xicor_test", n_digit = 2, seed = SEED)

既定では漸近理論に基づいた\(p\)値が計算されます。xicor()関数にmethod = "permutation"を指定することで、並べ替え検定による\(p\)値が計算されます：

コード例

XICOR::xicor(x, y, pvalue = TRUE, method = "permutation")$pval

plot_examples(title = "xicor_test_perm", n_digit = 2, seed = SEED)

なおヘルプ（?XICOR::xicor）には「漸近理論が十分優れているので、通常は並べ替え検定は必要ない」旨が記載されています。この例では、確かに両者の違いはあまりないようです。ヘルプから以下に引用します：

Usually, there is no need for the permutation test. The asymptotic theory is good enough.

まとめ

本稿では「Chatterjee相関」および「分割相関」を中心に、Rにおいてさまざまな相関尺度を算出し、結果を比較しました。

Chatterjee相関は\(Y = f(X)\)のような関数関係の検出に強く、予測の方向を加味した評価を与えます。独立性の仮説検定を行うことができ、計算も高速であることから、特に予測の観点に興味がある場面において強力かつ実用的な手法といえそうです。

分割相関は独立性からの乖離を評価し、数値だけでなくカテゴリ値にも適用できる統一的な尺度であることが特徴です。計算は比較的高速で、関数関係に限らず幅広い相関パターンを検出することができます。現行の実装に仮説検定は含まれないものの、大規模データにおける探索的解析において役立つことが期待されます。

比較対象として標本相関係数、MIC、距離相関およびHSICについても算出しました。本稿では詳しく触れられませんでしたが、手法によって高く評価されるパターンが異なることが確認できました。あわせて各種の独立性検定についても結果を比較しました。

非線形相関の世界は奥が深く、目的やデータ型、サンプルサイズなどを考慮して手法を選択することが重要と考えられます。本稿が皆様のデータ解析においてお役に立てれば幸いです。

参考リンク

本稿で参照した指標について、論文やスライドなど主要な一次情報、およびパッケージの公式サイトへのリンクをまとめます。

Chatterjee相関

• Chatterjee (2021) A new coefficient of correlation
• スライド資料
• XICORパッケージ

分割相関

• R研究集会2025 スライド資料
• gcorパッケージ

MIC

• Reshef et al. (2011) Detecting Novel Associations in Large Data Sets
• minervaパッケージ

距離相関

• Székely et al. (2007) Measuring and testing dependence by correlation of distances
• energyパッケージ

HSIC

• Gretton et al. (2005) Measuring Statistical Dependence with Hilbert-Schmidt Norms
• dHSICパッケージ

補遺

本稿の執筆にあたって作成したソースコードを以下に掲載します。無保証ですが、0BSDライセンス（日本語訳）のもとで自由にご利用いただけます。

サンプルデータ生成関数

# Original: Wikimedia Commons (CC0 public domain)
# https://commons.wikimedia.org/wiki/File:Correlation_examples2.svg
gen_example <- function(type = c("/", "U", "W", "X", "O", "88"),
                        n = 1000, rho = 0.8, seed = NULL) {
  requireNamespace("mvtnorm")

  type <- match.arg(type)

  if(!is.null(seed)) {
    set.seed(seed)
  }

  if (type == "/") {
    mat <- mvtnorm::rmvnorm(n, c(0, 0), matrix(c(1, rho, rho, 1), nrow = 2))
  } else if(type == "88") {
    nn <- table(rep(1:4, length.out = n))

    xy1 <- mvtnorm::rmvnorm(nn[1], c(3, 3))
    xy2 <- mvtnorm::rmvnorm(nn[2], c(-3, 3))
    xy3 <- mvtnorm::rmvnorm(nn[3], c(-3, -3))
    xy4 <- mvtnorm::rmvnorm(nn[4], c(3, -3))

    mat <- rbind(xy1, xy2, xy3, xy4)

  } else {
    x <- runif(n, -1, 1)

    if (type == "O") {
      xx <- sin(x * pi) + rnorm(n, 0, 1 / 8)
      yy <- cos(x * pi) + rnorm(n, 0, 1 / 8)
      mat <- cbind(xx, yy)

    } else {
      y <- switch(
        type,
        "U" = 2 * x^2 + runif(n, -1, 1),
        "W" = 4 * (x^2 - 1 / 2)^2 + runif(n, -1, 1) / 3,
        "X" = (x^2 + runif(n, 0, 1 / 2)) *
          sample(seq(-1, 1, 2), n, replace = TRUE),
        stop("unimplemented type")
      )
      mat <- cbind(x, y)
    }
  }

  colnames(mat) <- c("x", "y")
  data.frame(mat)
}

サンプルデータの散布図を描画する関数

plot_examples <- function(n = 1000, rho = 0.8, seed = NULL,
                          title = c("type", "cor", "gcor",
                                    "xicor", "xicor_test", "xicor_test_perm",
                                    "mic",
                                    "dcor", "dcor_test",
                                    "hsic", "hsic_test"),
                          n_digit = 2) {
  COL_STRONG <- "brown3"
  COL_WEAK   <- "green4"
  
  title <- match.arg(title)

  mfrow_orig <- par()$mfrow
  par(mfrow = c(2, 3))
  on.exit(par(mfrow = mfrow_orig))

  for(t in c("/", "U", "W", "X", "O", "88")) {
    d <- gen_example(type = t, n = n, rho = rho, seed = seed)

    col_main <- "black"

    if(title == "type") {
      main <- paste0("type: ", t)
    } else {
      val <- switch(title, 
        "cor"   = cor(d$x, d$y),
        "gcor"  = gcor::gcor(d$x, d$y),
        "xicor" = XICOR::xicor(d$x, d$y),
        "xicor_test" = XICOR::xicor(d$x, d$y, pvalue = TRUE)$pval,
        "xicor_test_perm" = XICOR::xicor(d$x, d$y, pvalue = TRUE, method = "permutation")$pval,
        "mic"   = minerva::mine(d$x, d$y)$MIC,
        "dcor"  = energy::dcor(d$x, d$y),
        "dcor_test" = energy::dcor.test(d$x, d$y, R = 1000)$p.value,
        "hsic"  = dHSIC::dhsic(d$x, d$y)$dHSIC,
        "hsic_test" = dHSIC::dhsic.test(d$x, d$y)$p.value
      )

      val_round <- round(val, digits = n_digit)

      if(title %in% c("xicor_test", "hsic_test", "dcor_test")) {
        if(val <  0.01) col_main <- COL_STRONG
        if(val >= 0.05) col_main <- COL_WEAK
      } else if(title != "hsic") {
        if(val_round <  1/3) col_main <- COL_WEAK
        if(val_round >= 2/3) col_main <- COL_STRONG
      }

      main <- val_round
    }

    plot(y ~ x, data = d, col = "darkblue", pch = 16, cex = 0.3,
      cex.main = 1.6, xaxt = "n", yaxt = "n", bty = "n",
      main = main, xlab = "", ylab = "", col.main = col_main)
  }
}