Reactive stat 内部の反復測定分散分析の実装です

反復測定分散分析（Repeated Measures ANOVA）は、同一の被験者に対して複数の条件での測定を行った場合に、条件間の差を検定するための統計手法です。

この記事では、反復測定分散分析の概要と、JavaScriptでの実装方法を説明します。

はじめに

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

使用したライブラリ

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

反復測定分散分析の概要

反復測定 ANOVA の説明

この手法は、時間経過に伴う変化や異なる条件下での変化を検討する際に特に役立ちます。

同じ対象群の複数の測定: 同じ個体や対象群が複数回の測定を受けることから、「反復測定」と呼ばれます
条件間の変化を分析: 薬の効果、トレーニングの影響、時間経過といった要因の効果を検討するのに適しています
対象群の個体差を考慮: 同じ個体からの複数の測定データを利用するため、個体間の変動を制御することが可能です

注意点:

データの欠損: 全ての時間点や条件でデータが欠損している場合、分析を行うことが難しい場合があります
球面性の仮定: 反復測定ANOVAでは、球面性 (sphericity) の仮定を満たす必要があります。これは、全ての条件間の差の分散が等しいことを意味します。この仮定が満たされない場合、自由度の調整 (例: Greenhouse-Geisser correction) が必要です
正規性: 反復測定ANOVAを適用する前に、データが正規分布に従っているかを確認する必要があります。正規分布に従っていない場合には、ノンパラメトリックな手法の利用を検討してください

計算式:

条件間の平方和 (Between-Conditions Sum of Squares)

この部分は、各条件の平均が全体の平均からどれだけ離れているかを示します。
```
SS_{\text{Between}} = n \sum_{i=1}^{k} (\bar{X}_{i} - \bar{X}_{\text{Grand}})^2
```
ここで、

$n$ は各条件の被験者数、
$k$ は条件の数、
$\bar{X}_{i}$ は第i条件の平均、
$\bar{X}_{\text{Grand}}$ は全体の平均です。
被験者間の平方和 (Between-Subjects Sum of Squares)

この部分は、各被験者の全条件にわたる平均が全体の平均からどれだけ離れているかを示します。
```
SS_{\text{Subjects}} = k \sum_{j=1}^{n} (\bar{X}_{j} - \bar{X}_{\text{Grand}})^2
```
ここで、
$\bar{X}_{j}$ は第j被験者の全条件にわたる平均です。
条件内の平方和 (Within-Conditions Sum of Squares)

この部分は、各被験者のスコアがその条件の平均からどれだけ離れているかを示します。

$$ SS_{\text{Within}} = \sum_{i=1}^{k} \sum_{j=1}^{n} (X_{ij} - \bar{X}_{i})^2 $$

ここで、
$X_{ij}$ は第i条件の第j被験者のスコアです。

最後に、これらの変動を元にF統計量を計算します。

$$ F = \frac{MS_{\text{Between}}}{MS_{\text{Within}}} $$

ここで、

$MS_{\text{Between}}$ は条件間の平均平方 (条件間の平方和を自由度で割ったもの)
$MS_{\text{Within}}$ は条件内の平均平方 (条件内の平方和を自由度で割ったもの)
です。

F統計量とその自由度からp値を計算します。

コード

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

html

<!DOCTYPE html>
<html>
<head>
  <meta charset="utf-8">
  <title>Repeated Measures ANOVA</title>
  <script src="https://cdn.jsdelivr.net/npm/jstat@latest/dist/jstat.min.js"></script>
</head>
<body>
  <h1>Repeated Measures ANOVA Example</h1>
  <p>Check the console for the result.</p>

  <script>
    function repeatedMeasuresANOVA(...groups) {
        // 各グループの値を数値に変換
        groups = groups.map((group) => group.map((x) => Number(x)));
        
        const nGroups = groups.length;
        if (nGroups < 2) {
            throw new Error("At least two groups are required.");
        }
        
        const n = groups[0].length;
        if (!groups.every(group => group.length === n)) {
            throw new Error("All groups must have the same number of subjects.");
        }

        const grandMean = jStat.mean(groups.flat());
        const groupMeans = groups.map(group => jStat.mean(group));

        // 各被験者の平均
        const subjectMeans = Array.from({length: n}, (_, i) => {
            return jStat.mean(groups.map(group => group[i]));
        });

        // 群間の平方和（Between-Groups）
        const ssBetween = groups.reduce((acc, group, index) => {
            return acc + Math.pow(groupMeans[index] - grandMean, 2) * n;
        }, 0);

        // 被験者内の平方和（Within-Subjects）
        const ssWithinSubjects = groups.reduce((acc, group, groupIndex) => {
            return acc + group.reduce((groupAcc, value, subjectIndex) => {
                return groupAcc + Math.pow(value - groupMeans[groupIndex], 2);
            }, 0);
        }, 0);

        // 被験者間の平方和（Between-Subjects）
        const ssSubjects = subjectMeans.reduce((acc, mean) => {
            return acc + Math.pow(mean - grandMean, 2) * nGroups;
        }, 0);

        // 群内の平方和（Within-Groups）
        const ssWithin = ssWithinSubjects - ssSubjects;

        // 群間の自由度
        const dfBetween = nGroups - 1;
        // 群内の自由度
        const dfWithin = (n - 1) * (nGroups - 1);
        // 被験者間の自由度
        const dfSubjects = n - 1;

        // 平均平方
        const msBetween = ssBetween / dfBetween;
        const msWithin = ssWithin / dfWithin;
        const msSubjects = ssSubjects / dfSubjects;

        // F統計量の計算
        const fStatistic = msBetween / msWithin;

        // p値の計算
        const pValue = 1 - jStat.centralF.cdf(fStatistic, dfBetween, dfWithin);

        return pValue;
    }

    // サンプルデータ
    const group1 = [254, 269, 292, 248, 258, 232, 332, 215, 285, 202, 168, 341, 329, 269, 232, 284, 213, 237, 231, 221, 286, 233, 177, 270, 286, 183, 341, 308, 270, 270, 176, 283, 276, 176, 231, 301, 173, 169, 344, 192, 188, 223, 287, 218, 280, 152, 334, 179, 197, 179, 200, 219, 249, 246, 256, 214, 295, 253, 317, 364, 258, 326, 201, 297, 227, 353, 302, 290, 309, 263, 261, 168, 328, 250, 254, 375, 258, 249, 204, 218, 92, 267, 277, 245, 174, 246, 279, 247, 147, 320, 255, 315, 171, 177, 150, 292, 187, 243, 217, 175];
    const group2 = [252, 255, 263, 272, 262, 237, 353, 223, 263, 200, 170, 331, 348, 275, 233, 276, 256, 248, 232, 230, 285, 236, 164, 313, 314, 174, 319, 294, 281, 266, 158, 278, 268, 180, 229, 270, 182, 153, 359, 192, 170, 224, 274, 248, 260, 168, 318, 195, 263, 213, 196, 224, 290, 287, 281, 227, 291, 230, 343, 359, 246, 288, 209, 288, 205, 345, 305, 291, 315, 247, 262, 179, 338, 223, 242, 409, 251, 230, 233, 210, 78, 312, 263, 259, 180, 231, 294, 260, 157, 328, 275, 325, 209, 184, 165, 309, 181, 203, 208, 154];
    const group3 = [266, 297, 314, 225, 282, 223, 337, 220, 304, 224, 157, 339, 349, 269, 255, 316, 198, 294, 224, 206, 266, 215, 233, 285, 310, 163, 334, 282, 279, 277, 174, 289, 269, 212, 225, 324, 159, 160, 393, 195, 187, 185, 295, 213, 266, 169, 351, 205, 209, 219, 246, 225, 258, 263, 271, 203, 289, 278, 339, 344, 236, 346, 212, 293, 203, 363, 275, 293, 311, 272, 272, 200, 340, 212, 241, 391, 244, 239, 197, 237, 64, 285, 321, 257, 190, 261, 264, 256, 180, 320, 241, 299, 142, 206, 185, 325, 207, 219, 222, 167];
    const group4 = [259, 297, 296, 263, 318, 245, 326, 224, 279, 205, 171, 365, 334, 304, 259, 276, 207, 228, 201, 206, 255, 233, 151, 285, 326, 177, 384, 316, 303, 261, 155, 298, 295, 199, 267, 303, 157, 167, 339, 172, 198, 226, 301, 213, 301, 182, 345, 166, 244, 203, 188, 226, 258, 303, 267, 264, 309, 263, 323, 359, 293, 346, 232, 313, 274, 365, 351, 269, 322, 257, 286, 165, 350, 243, 263, 380, 271, 288, 209, 184, 111, 305, 277, 229, 170, 253, 291, 296, 141, 328, 291, 326, 194, 188, 132, 306, 195, 262, 212, 176];

    const pValue = repeatedMeasuresANOVA(group1, group2, group3, group4);
    console.log("Repeated Measures ANOVA p-value: ", pValue);
  </script>
</body>
</html>

検証用 R コード

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

R

# サンプルデータ
data <- data.frame(
  Subject = factor(rep(1:100, times=4)),
  Time = factor(rep(1:4, each=100), , labels=c("A","B","C","D")),
  Value = c(254,269,292,248,258,232,332,215,285,202,168,341,329,269,232,284,213,237,231,221,286,233,177,270,286,183,341,308,270,270,176,283,276,176,231,301,173,169,344,192,188,223,287,218,280,152,334,179,197,179,200,219,249,246,256,214,295,253,317,364,258,326,201,297,227,353,302,290,309,263,261,168,328,250,254,375,258,249,204,218,92,267,277,245,174,246,279,247,147,320,255,315,171,177,150,292,187,243,217,175,252,255,263,272,262,237,353,223,263,200,170,331,348,275,233,276,256,248,232,230,285,236,164,313,314,174,319,294,281,266,158,278,268,180,229,270,182,153,359,192,170,224,274,248,260,168,318,195,263,213,196,224,290,287,281,227,291,230,343,359,246,288,209,288,205,345,305,291,315,247,262,179,338,223,242,409,251,230,233,210,78,312,263,259,180,231,294,260,157,328,275,325,209,184,165,309,181,203,208,154,266,297,314,225,282,223,337,220,304,224,157,339,349,269,255,316,198,294,224,206,266,215,233,285,310,163,334,282,279,277,174,289,269,212,225,324,159,160,393,195,187,185,295,213,266,169,351,205,209,219,246,225,258,263,271,203,289,278,339,344,236,346,212,293,203,363,275,293,311,272,272,200,340,212,241,391,244,239,197,237,64,285,321,257,190,261,264,256,180,320,241,299,142,206,185,325,207,219,222,167,259,297,296,263,318,245,326,224,279,205,171,365,334,304,259,276,207,228,201,206,255,233,151,285,326,177,384,316,303,261,155,298,295,199,267,303,157,167,339,172,198,226,301,213,301,182,345,166,244,203,188,226,258,303,267,264,309,263,323,359,293,346,232,313,274,365,351,269,322,257,286,165,350,243,263,380,271,288,209,184,111,305,277,229,170,253,291,296,141,328,291,326,194,188,132,306,195,262,212,176)
)

# 反復測定分散分析の実行
result <- aov(Value ~ Time + Error(Subject/Time), data=data)

# 結果の表示
summary(result)

注意点

このコードは、jStatライブラリに依存しています。jStatを読み込まないと動作しません
浮動小数点数の計算では、わずかな誤差が生じる可能性があります
反復測定分散分析を行う前に、以下の仮定を確認する必要があります
- 正規性：各条件の測定値が正規分布に従うこと
- 球面性：条件間の分散が等しいこと。球面性の仮定が満たされない場合、自由度を調整するなどの補正が必要になります

最後に

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

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

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