title: JavaScript での Brunner-Munzel検定 の実装
tags: JavaScript R 統計 検定 Brunner-Munzel
Reactive stat 内部の Brunner-Munzel検定 の実装 です
この記事では、Brunner-Munzel検定 の概要と、JavaScriptでの実装方法を説明します。
はじめに
弊社では ブラウザだけで使える無料統計ソフト Reactive stat を提供しています。
信頼性の高い R で統計解析し、その結果を AI が解説します！
その背景には、統計に苦労している医療者の助けになりたい、という気持ちがあります。
最終的な統計解析は R で行うのですが、レスポンスとサーバー負荷の改善のため、一部の統計計算はブラウザ内で行っています。
そのうち、汎用性の高い部分を共有させていただきたいと思います。
できるだけ R の出力と整合性を持つように javascript をインプリメントしてあるので、参考になれば幸いです。
使用したライブラリ
このコードでは、jStatライブラリを使用しています。jStatは、JavaScriptで統計計算を行うための便利なライブラリです。
jStatの詳細については、公式ドキュメント を参照してください。
全関数の一覧は jStat v1.9.3 Documentation を参照してください。
Brunner-Munzel検定 の概要
Brunner-Munzel検定は、2つの独立した群を比較するためのノンパラメトリック検定です。この検定は、従来のMann-Whitney U検定よりも柔軟で、分布の形状や分散が異なる場合にも適用できる特徴があります。
主な特徴
-
分布の仮定が緩やか
Brunner-Munzel検定は、2つの群の分布の形状や分散が異なる場合でも適用できます。これは、Mann-Whitney U検定が2群の分布形状と分散が同じであることを前提としているのとは対照的です。 -
順位に基づく検定
この検定は、データの実際の値ではなく、順位情報を使用します。具体的には、全てのデータを合わせて順位付けし、各群のデータがどの順位に位置するかを分析します。 -
高い検出力
データがMann-Whitney U検定の前提条件を満たさない場合、Brunner-Munzel検定は誤った結論を導き出すリスクが低く、真の差を検出する能力（検出力）が高くなります。 -
外れ値への耐性
順位を用いるため、極端な値（外れ値）の影響を受けにくい特徴があります。これは実際のデータ分析で重要な利点となります。
Mann-Whitney U検定との比較
-
Mann-Whitney U検定：2群の分布形状と分散が同じであることを前提とします。この条件が満たされる場合、計算が簡単で効率的な検定方法です。
-
Brunner-Munzel検定：分布の形状や分散に関する前提条件がより緩やかです。そのため、より幅広い状況で適用可能ですが、計算はやや複雑になります。
適用場面
Brunner-Munzel検定は以下のような状況で特に有効です：
- データの分布形状が2群で明らかに異なる場合
- 2群の分散（ばらつき）が大きく異なる場合
- データに外れ値が含まれている可能性がある場合
- サンプルサイズが小さい場合
小括
Brunner-Munzel検定は、Mann-Whitney U検定よりも柔軟なノンパラメトリック検定です。分布の形状や分散の違い、外れ値の存在などに対して頑健性があり、より広範なデータに適用できます。
従来のMann-Whitney U検定の前提条件が満たされない場合に、Brunner-Munzel検定の使用を検討すべきです。
コード
ライブラリ込みの html にしてありますので、
jsfiddle などにコピペして簡単に試せます。
<!DOCTYPE html>
<html lang="ja">
<head>
<meta charset="UTF-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Brunner-Munzel Test</title>
<script src="https://cdnjs.cloudflare.com/ajax/libs/jstat/1.9.5/jstat.min.js"></script>
</head>
<body>
<h1>Brunner-Munzel Test</h1>
<button onclick="runTest()">テストを実行</button>
<pre id="result"></pre>
<script>
// Brunner-Munzel 検定
// このテストは、2つの独立したサンプル間の確率的平等性を検定するために使用されます。
const brunnerMunzel = (x, y, alternative = "twoSided", alpha = 0.05) => {
// Remove NA values
x = x.filter(val => val !== null && !isNaN(val));
y = y.filter(val => val !== null && !isNaN(val));
const n1 = x.length;
const n2 = y.length;
// Calculate ranks
const allValues = [...x, ...y];
const ranks = jStat.rank(allValues);
const r1 = ranks.slice(0, n1);
const r2 = ranks.slice(n1);
// Calculate means
const m1 = jStat.mean(r1);
const m2 = jStat.mean(r2);
// Calculate pst (corrected as per R implementation)
const pst = (m2 - (n2 + 1) / 2) / n1;
// Calculate variances
const v1 = jStat.sum(r1.map((r, i) => Math.pow(r - jStat.rank(x)[i] - m1 + (n1 + 1) / 2, 2))) / (n1 - 1);
const v2 = jStat.sum(r2.map((r, i) => Math.pow(r - jStat.rank(y)[i] - m2 + (n2 + 1) / 2, 2))) / (n2 - 1);
// Calculate test statistic
const statistic = n1 * n2 * (m2 - m1) / (n1 + n2) / Math.sqrt(n1 * v1 + n2 * v2);
// Calculate degrees of freedom
const dfbm = Math.pow(n1 * v1 + n2 * v2, 2) / (Math.pow(n1 * v1, 2) / (n1 - 1) + Math.pow(n2 * v2, 2) / (n2 - 1));
// Calculate p-value
let pValue;
if (alternative === "greater") {
pValue = 1 - jStat.studentt.cdf(statistic, dfbm);
} else if (alternative === "less") {
pValue = jStat.studentt.cdf(statistic, dfbm);
} else {
pValue = 2 * Math.min(jStat.studentt.cdf(Math.abs(statistic), dfbm), 1 - jStat.studentt.cdf(Math.abs(statistic), dfbm));
}
// Calculate confidence interval
const tValue = jStat.studentt.inv(1 - alpha / 2, dfbm);
const ciLower = pst - tValue * Math.sqrt(v1 / (n1 * n2 * n2) + v2 / (n2 * n1 * n1));
const ciUpper = pst + tValue * Math.sqrt(v1 / (n1 * n2 * n2) + v2 / (n2 * n1 * n1));
const result = {
statistic,
df: dfbm,
pValue,
estimate: pst,
ci: [ciLower, ciUpper],
confLevel: 1 - alpha
};
return result;
};
function runTest() {
// サンプルデータ
const groupA = [1, 2, 3, 4, 5];
const groupB = [2, 3, 4, 5, 6];
const result = brunnerMunzel(groupA, groupB);
document.getElementById('result').textContent = JSON.stringify(result, null, 2);
}
</script>
</body>
</html>
検証用 R コード
rdrr にコピペして簡単に試せます。
library(lawstat)
# サンプルデータ
groupA <- c(1, 2, 3, 4, 5)
groupB <- c(2, 3, 4, 5, 6)
# Brunner-Munzel検定の実行
result <- brunner.munzel.test(groupA, groupB)
# 結果の表示
print(result)
注意点
- このコードは、jStatライブラリに依存しています。jStatを読み込まないと動作しません
- 浮動小数点数の計算では、わずかな誤差が生じる可能性があります
最後に
コードの内容はよく吟味し、R との食い違いが極力ないように気を付けていますが、もし間違いに気づかれた場合には指摘していただけますとありがたいです。
ただ、Chat GPT に読ませてその出力を検証もせず「問題があるので修正します」として貼り付けられると、他の方を混乱させてしまうのと、弊社の製品に対する信頼を損ねることにもつながりかねませんので、どうかそのようなことのないようにお願いいたします。
検証のための R のコードもつけてありますので、ご活用いただけますと幸いです。