Reactive stat 内部の Kruskal-Wallis検定の実装です

Kruskal-Wallis検定は、3つ以上の独立なグループ間の差異を検定するノンパラメトリック手法です。この記事では、Kruskal-Wallis検定の概要と、JavaScriptでの実装方法を説明します。

はじめに

弊社では ブラウザだけで使える無料統計ソフト Reactive stat を提供しています。
信頼性の高い R で統計解析し、その結果を AI が解説します！
その背景には、統計に苦労している医療者の助けになりたい、という気持ちがあります。
最終的な統計解析は R で行うのですが、レスポンスとサーバー負荷の改善のため、一部の統計計算はブラウザ内で行っています。
そのうち、汎用性の高い部分を共有させていただきたいと思います。
できるだけ R の出力と整合性を持つように javascript をインプリメントしてあるので、参考になれば幸いです。

使用したライブラリ

このコードでは、jStatライブラリを使用しています。jStatは、JavaScriptで統計計算を行うための便利なライブラリです。
jStatの詳細については、公式ドキュメントを参照してください。
全関数の一覧は jStat v1.9.3 Documentation を参照してください。

Kruskal-Wallis検定の概要

Kruskal-Wallis検定は以下の手順で行います：

ランクの計算

$\text{allValues} = \bigcup_{i=1}^k \text{groups}_i$
ここで、$k$ はグループの数です。

ランク合計

$R_i = \sum_{j=1}^{n_i} \text{rank}_{i,j}$
ここで、$R_i$ は第$i$グループのランク合計で、$n_i$ は第$i$グループのデータ点の数です。

Kruskal-Wallis 統計量 $H$ の計算

$H = \frac{12}{N(N+1)} \sum_{i=1}^k \frac{R_i^2}{n_i} - 3(N+1)$
ここで、$N$ は全データ点の数です。

同順位の補正

$T = \sum_{i=1}^{n_{\text{ties}}} (t_i^3 - t_i)$
ここで、$t_i$ は同順位のデータ点の数です。

補正値： $\text{correctionT} = 1 - \frac{T}{N^3 - N}$

補正済み Kruskal-Wallis 統計量

$\text{kruskalWallisH} = \frac{H}{\text{correctionT}}$

p値の計算

最後に、$\chi^2$ 分布を使用してp値を計算します。

コード

ライブラリ込みの html にしてありますので、
jsfiddle などにコピペして簡単に試せます。

html

<!DOCTYPE html>
<html>
<head>
  <meta charset="utf-8">
  <title>Kruskal-Wallis Test</title>
  <script src="https://cdn.jsdelivr.net/npm/jstat@latest/dist/jstat.min.js"></script>
</head>
<body>
  <h1>Kruskal-Wallis Test Example</h1>
  <p>Check the console for the result.</p>
  
  <script>
    //Kruskal-Wallis 検定
    function kruskalWallis(...groups) {
        // すべてのデータを連結し、ランクを計算
        const allValues = groups.flat();
        const ranks = jStat.rank(allValues);
        const N = allValues.length;

        // 各グループのランクの合計を計算
        let rankIndex = 0;
        const rankSums = groups.map(group => {
            return group.reduce((sum, _) => sum + ranks[rankIndex++], 0);
        });

        // Kruskal-Wallis 統計量 H を計算
        const sumOfSquares = rankSums.reduce((sum, rankSum, idx) => {
            return sum + (rankSum ** 2) / groups[idx].length;
        }, 0);
        const H = (12 / (N * (N + 1))) * sumOfSquares - 3 * (N + 1); // H value

        // 同順位が多い場合の補正
        // https://real-statistics.com/one-way-analysis-of-variance-anova/kruskal-wallis-test/
        // 全体の tied ranks の数を取得
        const T = jStat.unique(ranks).map(rank =>
            ranks.filter(r => r === rank).length
        ).filter(count => count > 1).map(t => Math.pow(t, 3) - t).reduce((acc, t) => acc + t, 0) / 12;
        //補正値を計算
        const correctionT = 1 / (1 - (T / (N ** 3 - N)));

        const kruskalWallisH = H / correctionT;
        return 1 - jStat.chisquare.cdf(kruskalWallisH, groups.length - 1); //  p value
    }

    // サンプルデータ
    const group1 = [1, 2, 3, 4, 5];
    const group2 = [2, 4, 6, 8, 10];
    const group3 = [3, 6, 9, 12, 15];

    const pValue = kruskalWallis(group1, group2, group3);
    console.log("Kruskal-Wallis Test p-value: ", pValue);
  </script>
</body>
</html>

検証用 R コード

rdrr にコピペして簡単に試せます。

R

# Kruskal-Wallis Test

# サンプルデータ
group1 <- c(1, 2, 3, 4, 5)
group2 <- c(2, 4, 6, 8, 10)
group3 <- c(3, 6, 9, 12, 15)

# Kruskal-Wallis 検定を実行
result <- kruskal.test(list(group1, group2, group3))

# 結果の表示
print(result)

注意点

このコードは、jStatライブラリに依存しています。jStatを読み込まないと動作しません
浮動小数点数の計算では、わずかな誤差が生じる可能性があります

最後に

コードの内容はよく吟味し、R との食い違いが極力ないように気を付けていますが、もし間違いに気づかれた場合には指摘していただけますとありがたいです。

ただ、Chat GPT に読ませてその出力を検証もせず「問題があるので修正します」として貼り付けられると、他の方を混乱させてしまうのと、弊社の製品に対する信頼を損ねることにもつながりかねませんので、どうかそのようなことのないようにお願いいたします。

検証のための R のコードもつけてありますので、ご活用いただけますと幸いです。