JavaScript での多重比較におけるp値補正の実装

# 多重比較の補正について

## はじめに

統計解析において、3群以上の比較を行う際に各群の組み合わせごとに検定を行うと、検定の回数が増えてしまい、多重比較の問題が生じます。

多重比較とは、複数の仮説を同時に検定する際に生じる問題です。検定を繰り返すことで、第一種の過誤（帰無仮説が正しいにも関わらず棄却してしまう過誤）を犯す確率が高くなってしまいます。

要するに、3群だと3回、4群だと6回、5群だと10回というように、群の数が増えるごとに検定の回数が増加していきます。

一般に、$k$群の比較を行う場合、検定の回数は以下の式で求められます。

検定の回数: $\_{k}C_{2} = \frac{k(k-1)}{2}$

例えば、7群の比較を行う場合、21回の検定が行われるわけですから、そのうちの1回くらいは偶然 p<0.05 になりそうです。

そのため、多重比較の補正を行う必要が生じます。
Bonferroni法やHolm法などの補正方法を用いることで、検定の回数に応じて p値を適切に調整し、第一種の過誤を制御することができます。

ただし、これらの補正方法は保守的になりがちなので、検定のパワーが低下してしまい、有意とすべき違いを検出できない危険性もあります。

そのため、研究の目的や性質に応じて、適切な補正方法を選択する必要があります。
多重比較の問題に対処するため、様々な多重比較の補正方法が提案されていますので紹介します。

## Bonferroni補正

Bonferroni補正は、最も簡単で保守的な補正方法の一つであり、元の p値を比較数（検定の数）で乗算することで補正します。
これで補正した結果が p<0.05 であれば文句の付けようがありません。

- 利点: 簡単で明確。第一種の過誤のリスクが高まることがない
- 欠点: 非常に保守的で、第二種の過誤（帰無仮説が誤っているにも関わらず棄却しない過誤）のリスクが高まる

```html:html
<!DOCTYPE html>
<html>
<body>
<script>
    function bonferroniCorrection(pValues) {
        const n = pValues.length;
        return pValues.map(val => Math.min(1, n * val));
    }
    const pValues = [0.01, 0.03, 0.02, 0.04, 0.05];
    console.log("元のp値:", pValues);
    console.log("Bonferroni補正後のp値:", bonferroniCorrection(pValues));
</script>
</body>
</html>
```

## Holm法

Holm法は、p値を昇順にソートし、最小のp値からBonferroni補正を適用する方法です。
最小のp値が補正後のしきい値を超えた場合、そのp値とその後のすべてのp値を棄却します。

- 利点: Bonferroni補正よりもパワーが高い
- 欠点: まだ保守的な場合がある

```html:html
<!DOCTYPE html>
<html>
<body>
<script>
    function holmCorrection(pValues) {
        const m = pValues.length; // 検定の数

        // p値を昇順にソートし、インデックスも記憶する
        const sortedIndices = Array.from(Array(m).keys()).sort((a, b) => pValues[a] - pValues[b]);
        const sortedPValues = sortedIndices.map(i => pValues[i]);

        const correctedPValues = new Array(m);
        for (let i = 0; i < m; i++) {
            correctedPValues[i] = Math.min(1, sortedPValues[i] * (m - i));
        }

        // 補正後のp値が前のp値より小さくなっていないことを確認
        for (let i = 1; i < m; i++) {
            correctedPValues[i] = Math.max(correctedPValues[i], correctedPValues[i - 1]);
        }

        // 元の順序に戻す
        const finalPValues = new Array(m);
        for (let i = 0; i < m; i++) {
            finalPValues[sortedIndices[i]] = correctedPValues[i];
        }

        return finalPValues;
    }
    const pValues = [0.01, 0.03, 0.02, 0.04, 0.05];
    console.log("元のp値:", pValues);
    console.log("Holm補正後のp値:", holmCorrection(pValues));
</script>
</body>
</html>
```

## その他の補正方法

Reactive statsでは、他にも以下のような補正方法を実装しています。

- Hochbergの方法: p値を降順にソートし、より小さい p値に対してより厳しいBonferroni補正を適用する方法
- Hommelの方法: Bonferroni法よりも柔軟に、複数の検定に対してより効率的に第一種の過誤を制御する補正方法
- Benjamini-Hochberg法: 偽発見率（FDR）をコントロールする方法
- Benjamini-Yekutieli法: Benjamini-Hochberg法よりも保守的なFDRの制御法

### R での実装

R では、以下のコマンドで p.adjust() のソースコードを確認できます。

```R:R
getAnywhere(p.adjust)
```

## 検証用 R コード

[rdrr](https://rdrr.io/snippets/) にコピペして簡単に試せます。

```R:R
# サンプルデータ
p_values <- c(0.01, 0.03, 0.02, 0.04, 0.05)

# Bonferroni補正
bonferroni_correction <- function(p_values) {
  p.adjust(p_values, method = "bonferroni")
}

# Holm法
holm_correction <- function(p_values) {
  p.adjust(p_values, method = "holm")
}

# Hochbergの方法
hochberg_correction <- function(p_values) {
  p.adjust(p_values, method = "hochberg")
}

# Hommelの方法
hommel_correction <- function(p_values) {
  p.adjust(p_values, method = "hommel")
}

# Benjamini-Hochberg法
benjamini_hochberg_correction <- function(p_values) {
  p.adjust(p_values, method = "BH")
}

# Benjamini-Yekutieli法
benjamini_yekutieli_correction <- function(p_values) {
  p.adjust(p_values, method = "BY")
}

# 結果の表示
cat("元のp値:", p_values, "\n")
cat("Bonferroni補正後のp値:", bonferroni_correction(p_values), "\n")
cat("Holm補正後のp値:", holm_correction(p_values), "\n")
cat("Hochberg補正後のp値:", hochberg_correction(p_values), "\n")
cat("Hommel補正後のp値:", hommel_correction(p_values), "\n")
cat("Benjamini-Hochberg補正後のp値:", benjamini_hochberg_correction(p_values), "\n")
cat("Benjamini-Yekutieli補正後のp値:", benjamini_yekutieli_correction(p_values), "\n")
```

## まとめ

多重比較の問題は、統計解析を行う上で重要な問題の一つです。
Bonferroni法やHolm法などの補正方法を適切に用いることで、第一種の過誤を制御しつつ、検定のパワーを保つことができます。

ただし、これらの補正方法は保守的になりがちなので、状況に応じて適切な方法を選択する必要があります。

## 最後に

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

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

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