Reactive stat 内部の Wilcoxon の符号付順位検定 の実装 です
Wilcoxon の符号付順位検定（Wilcoxon Signed-Rank Test）は、対応のある2群のデータに対して、中央値に差があるかどうかを検定する非パラメトリック検定です。データが正規分布に従わない場合や、サンプルサイズが小さい場合に使用されます。
この記事では Wilcoxon の符号付順位検定 の概要と、JavaScriptでの実装方法を説明します。
はじめに
最近、弊社より ブラウザだけで使える無料統計ソフト Reactive stat をリリースしました。
信頼性の高い R で統計解析し、その結果を AI が解説します！
その背景には、統計に苦労している医療者の助けになりたい、という気持ちがあります。
最終的な統計解析は R で行うのですが、レスポンスとサーバー負荷の改善のため、一部の統計計算はブラウザ内で行っています。
そのうち、汎用性の高い部分を共有させていただきたいと思います。
できるだけ R の出力と整合性を持つように javascript をインプリメントしてあるので、参考になれば幸いです。
使用したライブラリ
このコードでは、jStatライブラリを使用しています。jStatは、JavaScriptで統計計算を行うための便利なライブラリです。
jStatの詳細については、公式ドキュメント を参照してください。
全関数の一覧は jStat v1.9.3 Documentation を参照してください。
Wilcoxon の符号付順位検定 の概要
Wilcoxon の符号順位検定は、対応のある2群間の中央値の差を検定する非パラメトリック手法です。
2つの群の間で繰り返し測定が行われたデータやペアのデータに対して検定を行います。
例えば、治療前後の変化を調査したい場合などに使用されます。
データが正規分布をしていない場合や、サンプルサイズが小さい場合でも使用可能です。また、外れ値の影響を受けにくいという特徴も持っています。
注意事項
この検定は中央値の差を検定するものであるため、平均値の差に関する情報は提供しません。
対応のないデータには適用できません。対応のある群間の検定手法をご利用ください。
計算式
- ２つのデータセット $A$ と $B$ から、差の値を計算します
$$d_i = A_i - B_i$$
- これらの差 $d_i$ の絶対値に順位をつけます
$$\text{ranks} = \text{rank}(|d_i|)$$
- 差の符号に応じて、順位に正または負の符号をつけます
$$
\text{signedRanks}_i =
\begin{cases}
-\text{ranks}_i & \text{if } d_i < 0 \
\text{ranks}_i & \text{otherwise}
\end{cases}
$$
- 正と負の符号順位の和を計算します
$$\text{posSum} = \sum_{i=1}^{n} \text{signedRanks}_i^{+}$$
$$\text{negSum} = \sum_{i=1}^{n} \text{signedRanks}_i^{-}$$
- $W$ 統計量は、正と負の符号順位の和の絶対値の小さい方として計算します
$$W = \min(|\text{posSum}|, |\text{negSum}|)$$
- $W$ の期待値と分散を計算します
$$E(W) = \frac{n(n + 1)}{4}$$
$$\text{Var}(W) = \frac{n(n + 1)(2n + 1)}{24}$$
- サンプルサイズが大きい場合や同順位が存在する場合には、正規近似を使用してp値を計算します。この際、**連続性補正（continuity correction）**を適用します。
連続性補正値は以下のように計算されます。
$$
\text{連続性補正値} = \begin{cases}
-0.5, & \text{if } W < E(W) \
0.5, & \text{if } W \geq E(W)
\end{cases}
$$
この連続性補正値を用いて、標準化されたW値（Z）を次のように計算します。
$$
Z = \frac{W - E(W) - \text{連続性補正値}}{\sqrt{\text{Var}(W)}}
$$
直観的には、$W$が期待値より小さい場合は補正値として$-0.5$を減算することで$Z$値を小さくし、$W$が期待値以上の場合は補正値として$0.5$を減算することで$Z$値を大きくします。これにより、$W$が極端な値をとる確率を適切に調整し、正規近似による計算結果を正確な分布に近づけることができます。
- 最後に、標準正規分布を用いてZ値からp値を求めます
コード
ライブラリ込みの html にしてありますので、
jsfiddle などにコピペして簡単に試せます。
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8">
<title>Wilcoxon Signed-Rank Test</title>
<script src="https://cdn.jsdelivr.net/npm/jstat@latest/dist/jstat.min.js"></script>
</head>
<body>
<h1>Wilcoxon Signed-Rank Test Example</h1>
<p>Check the console for the result.</p>
<script>
function wilcoxonSignedRank(dataA, dataB) {
if (dataA.length !== dataB.length) {
throw new Error("Both arrays should have the same length");
}
const differences = dataA.map((value, index) => value - dataB[index]);
const ranks = jStat.rank(differences.map(Math.abs));
let signedRanks = ranks.map((rank, index) => differences[index] < 0 ? -rank : rank);
let posSum = signedRanks.reduce((sum, rank) => (rank > 0 ? sum + rank : sum), 0);
let negSum = signedRanks.reduce((sum, rank) => (rank < 0 ? sum + rank : sum), 0);
let W = Math.min(Math.abs(posSum), Math.abs(negSum));
const n = dataA.length;
const expectedW = n * (n + 1) / 4;
const rankCounts = Array.from(new Set(ranks)).map(rank => ranks.filter(r => r === rank).length);
const sumOfSquaredRanks = rankCounts.map(count => Math.pow(count, 3) - count).reduce((acc, val) => acc + val, 0);
const varianceW = (n * (n + 1) * (2 * n + 1)) / 24 - sumOfSquaredRanks / (2 * n); // Correct for tied ranks
// 連続性補正の適用
const correction = W < expectedW ? -0.5 : 0.5;
const Z = (W - expectedW - correction) / Math.sqrt(varianceW);
return jStat.ztest(Z, 2); // p value
}
// サンプルデータ
const before = [31, 32, 43, 54, 65, 72, 80, 90, 92];
const after = [22, 26, 27, 38, 49, 60, 65, 78, 84];
const pValue = wilcoxonSignedRank(before, after);
console.log("Wilcoxon Signed-Rank Test p-value: ", pValue);
</script>
</body>
</html>
検証用 R コード
rdrr にコピペして簡単に試せます。
# Wilcoxon Signed-Rank Test
# サンプルデータ
before <- c(31, 32, 43, 54, 65, 72, 80, 90, 92)
after <- c(22, 26, 27, 38, 49, 60, 65, 78, 84)
# Wilcoxon の符号付順位検定を実行
result <- wilcox.test(before, after, paired = TRUE)
# 結果を表示
print(result)
注意点
- このコードは、jStatライブラリに依存しています。jStatを読み込まないと動作しません
- 浮動小数点数の計算では、わずかな誤差が生じる可能性があります
- Wilcoxon の符号付順位検定は、中央値の差を検定するため、平均値の差を検出することはできません
- データが対応していることを確認してください。対応のないデータに対しては、Wilcoxon の順位和検定（Mann-Whitney U検定）を使用してください
- このコードでは、常に正規近似を使用し、連続補正を適用しています
- 正確な分布を使用するコードは含まれていません
- R でのデフォルトの動作は、サンプルサイズが大きい（両方のグループの合計が50以上）か同順位（ties）が存在する場合、正規近似が使用され、連続性補正が適用されます。そうでない場合、正確な方法が使用され、連続性補正は適用されません
- そのため、得られる結果は、R の wilcox.test() とわずかに異なっています
最後に
コードの内容はよく吟味し、R との食い違いが極力ないように気を付けていますが、もし間違いに気づかれた場合には指摘していただけますとありがたいです。
ただ、Chat GPT に読ませてその出力を検証もせず「問題があるので修正します」として貼り付けられると、他の方を混乱させてしまうのと、弊社の製品に対する信頼を損ねることにもつながりかねませんので、どうかそのようなことのないようにお願いいたします。
検証のための R のコードもつけてありますので、ご活用いただけますと幸いです。