Reactive stat 内部の t検定と Welch の t検定の実装 です
t検定は、2つのグループの平均値の違いを検証する際の基本的な手法です。
この記事では、t検定と Welch の t検定の概要と、JavaScriptでの実装方法を説明します。
はじめに
弊社では ブラウザだけで使える無料統計ソフト Reactive stat を提供しています。
信頼性の高い R で統計解析し、その結果を AI が解説します！
その背景には、統計に苦労している医療者の助けになりたい、という気持ちがあります。
最終的な統計解析は R で行うのですが、レスポンスとサーバー負荷の改善のため、一部の統計計算はブラウザ内で行っています。
そのうち、汎用性の高い部分を共有させていただきたいと思います。
できるだけ R の出力と整合性を持つように javascript をインプリメントしてあるので、参考になれば幸いです。
使用したライブラリ
このコードでは、jStatライブラリを使用しています。jStatは、JavaScriptで統計計算を行うための便利なライブラリです。
jStatの詳細については、公式ドキュメント を参照してください。
全関数の一覧は jStat v1.9.3 Documentation を参照してください。
t検定と Welch の t検定の概要
t検定
t検定の前提条件:
- 正規性：データが正規分布に従っていること。
- 分散の等質性：2つの群の分散が等しいこと。
t値の計算:
$t = \frac{\text{標本平均の差}}{\text{標準誤差}}$
自由度:
- 自由度は、分析で利用できるデータの「自由な」数を示します。
- 2群のt検定の場合、自由度は $df = n_1 + n_2 - 2$ です。ここで、$n_1$ と $n_2$ はそれぞれのグループのサンプルサイズです。
効果量:
- t検定は有意性を検定しますが、実際の効果の大きさ (効果量) も重要です
- 一般的な効果量の尺度には Cohenのd があります
Welch の t検定
Welch の t検定の前提条件:
- 正規性：データが正規分布に従っていること。
- 分散の不等質性：2つの群の分散が等しくない場合に適用されます。
Welch の t値の計算:
$t = \frac{\text{標本平均の差}}{\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}}$
ここで、$s_1^2$ と $s_2^2$ はそれぞれのグループの不偏分散、$n_1$ と $n_2$ はそれぞれのグループのサンプルサイズです。
自由度 (Welch-Satterthwaite の自由度):
$df = \frac{(\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2})^2}{\frac{(s_1^2/n_1)^2}{n_1-1} + \frac{(s_2^2/n_2)^2}{n_2-1}}$
注意点:
- t検定とWelchのt検定は、どちらも2つのグループ間の平均値の差を検定します。
- 分散の等質性が満たされる場合はt検定を、満たされない場合はWelchのt検定を使用します。
- 3つ以上のグループを比較する場合は、ANOVAや多重比較の手法を検討してください。
- データが正規分布に従わない場合は、Mann-Whitney U検定などの非パラメトリック検定を検討してください。
コード
ライブラリ込みの html にしてありますので、
jsfiddle などにコピペして簡単に試せます。
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8">
<title>t-test and Welch's t-test</title>
<script src="https://cdn.jsdelivr.net/npm/jstat@latest/dist/jstat.min.js"></script>
</head>
<body>
<h1>t-test and Welch's t-test Example</h1>
<p>Check the console for the result.</p>
<script>
// t検定
function ttest(groupA, groupB) {
const mean1 = jStat.mean(groupA);
const mean2 = jStat.mean(groupB);
const sd1 = jStat.stdev(groupA, true); // true は不偏標準偏差を取得
const sd2 = jStat.stdev(groupB, true);
// サンプルサイズ
const n1 = groupA.length;
const n2 = groupB.length;
// t統計量の計算
const tscore = (mean1 - mean2) / Math.sqrt((sd1 ** 2 / n1) + (sd2 ** 2 / n2));
// p値の計算
const df = n1 + n2 - 2; // 自由度
const p = 2 * (1 - jStat.studentt.cdf(Math.abs(tscore), df));
return p;
}
// Welchのt検定
function welchTTest(groupA, groupB) {
const meanA = jStat.mean(groupA);
const meanB = jStat.mean(groupB);
const varA = jStat.variance(groupA, true); // true は不偏分散を取得
const varB = jStat.variance(groupB, true);
// サンプルサイズ
const nA = groupA.length;
const nB = groupB.length;
// t統計量の計算
const tscore = (meanA - meanB) / Math.sqrt((varA / nA) + (varB / nB));
// 自由度の計算（Welch-Satterthwaite equation）
const df = ((varA / nA) + (varB / nB)) ** 2 /
((varA / nA) ** 2 / (nA - 1) + (varB / nB) ** 2 / (nB - 1));
// p値の計算
const p = 2 * (1 - jStat.studentt.cdf(Math.abs(tscore), df));
return p;
}
// サンプルデータ
const group1 = [85, 90, 92, 87, 99];
const group2 = [80, 82, 78, 83, 79];
const pValueTtest = ttest(group1, group2);
console.log("t-test p-value: ", pValueTtest);
const pValueWelch = welchTTest(group1, group2);
console.log("Welch's t-test p-value: ", pValueWelch);
</script>
</body>
</html>
検証用 R コード
rdrr にコピペして簡単に試せます。
# t-test and Welch's t-test
# サンプルデータ
group1 <- c(85, 90, 92, 87, 99)
group2 <- c(80, 82, 78, 83, 79)
# t検定を実行 (var.equal = TRUE で等分散を仮定)
result_ttest <- t.test(group1, group2, var.equal = TRUE)
# Welchのt検定を実行 (var.equal = FALSE で等分散を仮定しない)
result_welch <- t.test(group1, group2, var.equal = FALSE)
# 結果の表示
print("t-test result:")
print(result_ttest)
print("Welch's t-test result:")
print(result_welch)
注意点
- このコードは、jStatライブラリに依存しています。jStatを読み込まないと動作しません
- 浮動小数点数の計算では、わずかな誤差が生じる可能性があります
- R の t.test 関数は、デフォルトで Welch の t検定を行います。等分散を仮定する場合は、var.equal = TRUE を指定する必要があります
- サンプルサイズが十分大きく、等分散性の検定の結果が有意でない場合は、通常のt検定を使用できます。
- サンプルサイズが小さい場合や等分散性の検定の結果が有意である場合は、Welchのt検定を使用するのが安全です。
- ただし、サンプルサイズが非常に大きい場合、わずかな分散の違いでも等分散性の検定で有意になることがあります。そのような場合、実質的な分散の差が小さければ、通常のt検定を使用しても問題ありません。
最後に
コードの内容はよく吟味し、R との食い違いが極力ないように気を付けていますが、もし間違いに気づかれた場合には指摘していただけますとありがたいです。
ただ、Chat GPT に読ませてその出力を検証もせず「問題があるので修正します」として貼り付けられると、他の方を混乱させてしまうのと、弊社の製品に対する信頼を損ねることにもつながりかねませんので、どうかそのようなことのないようにお願いいたします。
検証のための R のコードもつけてありますので、ご活用いただけますと幸いです。