この回のコードは,すべて tidyverse
を読み込んだ状態で動きます。
RStudioを起動したら,まずこれを実行してください。
毎回必要です。読み込まないまま
filter()やggplot()を使うと, 「関数が見つかりません」というエラーになります。
3群(A, B, C)を比較するとき,t検定を繰り返すと:
第一種の過誤(偽陽性)の確率が膨らむ!
\[1 - (1 - 0.05)^3 = 1 - 0.857 = 0.143\]
→ 3回の検定では,帰無仮説が正しいのに少なくとも1回有意になる確率が14.3%!
1回の検定で偽陽性を出す確率は \(\alpha = 0.05\) です。裏を返せば,偽陽性を出さない確率は \(1 - \alpha = 0.95\) です。
\(k\) 回の検定すべてで偽陽性を出さない確率を考えます。ここで「\(k\) 個の確率を掛け算してよい」と言えるのは,\(k\) 回の検定が互いに独立なときだけです。
\[P(k\text{回とも偽陽性なし}) = (1-\alpha)^k\]
求めたいのは「少なくとも1回は偽陽性」の確率なので,余事象を取ります。
\[P(\text{少なくとも1回偽陽性}) = 1 - (1-\alpha)^k\]
\(\alpha = 0.05\),\(k = 3\) を入れると \((0.95)^3 = 0.857375\) なので \(1 - 0.857375 = 0.142625\),つまり約 \(14.3\%\) です。
ただし,この式は目安です。 3群の総当たり比較 A対B,A対C,B対C は,同じ群のデータを使い回しています。たとえば A対B と A対C は,どちらも A のデータを共有しています。だから3つの検定は独立ではなく,正に相関します。
実際に確かめてみます。3群の母平均がまったく同じ(差がゼロ)という状況を2万回作り,そのたびに総当たりのt検定3回と,分散分析1回を実行します。
set.seed(2026)
R <- 20000; n <- 50
g <- factor(rep(1:3, each = n))
p_min <- numeric(R) # t検定3回のうち最小のp値
p_aov <- numeric(R) # 分散分析のp値
for (r in 1:R) {
y <- rnorm(3 * n) # 3群すべて同じ分布(差はゼロ)
a <- y[1:n]; b <- y[(n+1):(2*n)]; c <- y[(2*n+1):(3*n)]
p_min[r] <- min(t.test(a, b, var.equal = TRUE)$p.value,
t.test(a, c, var.equal = TRUE)$p.value,
t.test(b, c, var.equal = TRUE)$p.value)
p_aov[r] <- anova(aov(y ~ g))[["Pr(>F)"]][1]
}
c(t検定3回 = mean(p_min < .05),
分散分析1回 = mean(p_aov < .05),
独立を仮定した式 = 1 - 0.95^3)## t検定3回 分散分析1回 独立を仮定した式
## 0.122950 0.048400 0.142625
総当たりのt検定3回では約 \(12\%\),独立を仮定した式の \(14.3\%\) より少し小さくなりました。検定どうしが連動しているぶん,「どれかが当たる」機会が減るからです。
この式は正確な偽陽性率ではなく,比較の数を増やすとどれくらいの勢いで誤りが膨らむかを示す目安として読んでください。
一方,分散分析を1回だけ使えば約 \(5\%\) で,狙いどおりに保たれています。これが「1回の検定で済ませる」ことの意味です。
このコードは2万回ループするので,実行に20秒ほどかかります。
set.seed(2026)を付けてあるので,手元で実行しても上とまったく同じ数値になります(set.seed()については第11回)。
使った数学:余事象の確率と,独立な事象の積。掛け算をした瞬間に独立を仮定しています。
Analysis of Variance(ANOVA): - 3群以上を1回の検定で同時に比較 - 全体の第一種の過誤率を α = .05 に保つ
\[SS_{全体} = SS_{群間} + SS_{群内}\]
| 変動の種類 | 意味 |
|---|---|
| 群間変動 (SS_between) | 群の平均値のばらつき(群の効果) |
| 群内変動 (SS_within) | 各群内のばらつき(誤差) |
第 \(j\) 群の第 \(i\) 番目の観測値を \(x_{ij}\),第 \(j\) 群の平均を \(\bar{x}_j\),全体平均を \(\bar{x}\) と書きます。
各観測値の全体平均からのズレを,2つに分けます。あいだに \(\bar{x}_j\) を挟んで,足して引くだけです。
\[x_{ij} - \bar{x} = \underbrace{(x_{ij} - \bar{x}_j)}_{\text{群内のズレ}} + \underbrace{(\bar{x}_j - \bar{x})}_{\text{群の位置のズレ}}\]
両辺を二乗して全部足します。\((a+b)^2 = a^2 + 2ab + b^2\) なので,交差項が出てきます。
\[\sum_j \sum_i (x_{ij}-\bar{x})^2 = \sum_j \sum_i (x_{ij}-\bar{x}_j)^2 + 2\sum_j \sum_i (x_{ij}-\bar{x}_j)(\bar{x}_j-\bar{x}) + \sum_j \sum_i (\bar{x}_j-\bar{x})^2\]
この交差項が0になるのが要点です。\(j\) を固定すると \((\bar{x}_j - \bar{x})\) は定数なので,内側の和の外に出せます。
\[2\sum_j (\bar{x}_j-\bar{x}) \sum_i (x_{ij}-\bar{x}_j)\]
ところが \(\sum_i (x_{ij} - \bar{x}_j)\) は,その群の中での偏差の合計です。平均からのズレを足し合わせれば,プラスとマイナスが打ち消し合って必ず0になります。だから交差項は丸ごと消えます。
残るのが
\[SS_{全体} = SS_{群内} + SS_{群間}\]
です。
数値で確かめる:iris のがく片の長さでは
\(SS_{全体} = 102.1683\),\(SS_{群間} = 63.2121\),\(SS_{群内} = 38.9562\) です。足すと \(63.2121 + 38.9562 = 102.1683\)
となり,ぴったり一致します。Rで差を取ると \(10^{-13}\)
程度の値が残りますが,これは小数を2進法で扱うときの計算機の丸め誤差で,数学的には厳密に0です。
使った数学:二乗の展開と,偏差の和が0になる性質。分散分析の「分散を分解する」という言い方は,この恒等式のことを指しています。
\[F = \frac{MS_{群間}}{MS_{群内}} = \frac{群間の平均二乗}{群内の平均二乗}\]
平方和をそのまま比べることはできません。群間は \(k-1\) 個,群内は \(N-k\) 個の自由な情報から作られていて,足し合わせた個数が違うからです。そこでそれぞれを自由度で割り,1個あたりのばらつきに直します。
\[MS_{群間} = \frac{SS_{群間}}{k-1}, \qquad MS_{群内} = \frac{SS_{群内}}{N-k}\]
ここで帰無仮説「すべての群の母平均が等しい」が正しいとしましょう。このとき群の違いは存在しないので,群平均のばらつきも,群内のばらつきも,どちらも同じ母分散 \(\sigma^2\) を映しているだけになります。つまり \(MS_{群間}\) と \(MS_{群内}\) は,同じ \(\sigma^2\) を2通りのやり方で推定したものです。
同じものを2通りに推定して比を取れば,値は1のあたりに来ます。
\[F = \frac{MS_{群間}}{MS_{群内}} \approx 1 \quad (\text{帰無仮説が正しいとき})\]
逆に群の平均が本当に違えば,\(SS_{群間}\) だけが大きくなります。群内のばらつきは群の違いとは無関係なので変わりません。だから \(F\) は1より大きい側へ動きます。\(F\) が1から離れて大きいほど,群の違いは誤差だけでは説明しにくいということです。
数値で確かめる:iris では \(MS_{群間} = 63.2121/2 = 31.6061\),\(MS_{群内} = 38.9562/147 = 0.2650\) なので
\(F = 31.6061/0.2650 =
119.26\)。1からは大きく離れています。
同じ分解から \(\eta^2\) も読めます。
\[\eta^2 = \frac{SS_{群間}}{SS_{全体}} = \frac{63.2121}{102.1683} = 0.619\]
全体のばらつきのうち \(61.9\%\) が「種の違い」で説明でき,残りが群内の誤差だという意味です。\(F\) が「誤差と比べて大きいか」を見るのに対し,\(\eta^2\) は「全体の何割を説明するか」を見ています。同じ2つの平方和から,別の問いに答えているわけです。
使った数学:自由度で割って平均に直すことと,比を取ること。なお \(F\) 分布の密度関数そのものを導くのは,この科目の範囲外です。
# データの確認
iris |>
group_by(Species) |>
summarise(
n = n(),
mean = mean(Sepal.Length),
sd = sd(Sepal.Length)
)# 箱ひげ図で確認
ggplot(iris, aes(x = Species, y = Sepal.Length, fill = Species)) +
geom_boxplot() +
geom_jitter(width = 0.1, alpha = 0.3) +
scale_fill_brewer(palette = "Set2") +
labs(title = "種ごとのがく片の長さ",
x = "種", y = "がく片の長さ (cm)") +
theme_bw() +
theme(legend.position = "none")## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 63.21 31.606 119.3 <2e-16 ***
## Residuals 147 38.96 0.265
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 63.21 31.606 119.3 <2e-16 ***
Residuals 147 38.96 0.265
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
| 項目 | 意味 |
|---|---|
Df |
自由度(群間 = k-1,群内 = N-k) |
Sum Sq |
平方和 |
Mean Sq |
平均二乗(平方和÷自由度) |
F value |
F統計量 |
Pr(>F) |
p値 |
p < .05 → 帰無仮説を棄却(少なくとも1組の群に有意差がある)
ただし,どの群とどの群が異なるかはこれだけではわかりません。 → 次回の多重比較へ
| 前提条件 | 確認方法 |
|---|---|
| 各群が正規分布に従う | Shapiro-Wilk 検定 |
| 各群の分散が等しい(等分散性) | Bartlett 検定 / Levene 検定 |
| 観測値の独立性 | データの収集方法を確認(検定では確かめられない) |
# 各群ごとに正規性を確認する
iris |>
group_by(Species) |>
summarise(
W = shapiro.test(Sepal.Length)$statistic,
p = shapiro.test(Sepal.Length)$p.value
)3群ともp値が大きく,正規分布から外れているとは言えません。
##
## Bartlett test of homogeneity of variances
##
## data: Sepal.Length by Species
## Bartlett's K-squared = 16.006, df = 2, p-value = 0.0003345
Levene 検定はcarパッケージのleveneTest()で実行できます。 標準のRには入っていないので,この授業では標準で使える Bartlett 検定を使います。 (Levene検定のほうが正規性からのズレに強い,という違いがあります。)
結果は p = .0003 でした。各群の標準偏差を見ると, setosa 0.35/versicolor 0.52/virginica 0.64 と,分散で3倍以上の違いがあります。
等分散を仮定しない Welch の ANOVA も実行できます。
##
## One-way analysis of means (not assuming equal variances)
##
## data: Sepal.Length and Species
## F = 138.91, num df = 2.000, denom df = 92.211, p-value < 2.2e-16
\(F = 138.91\)(分母の自由度 \(92.21\)),\(p
< .001\) です。 aov() の \(F = 119.26\)
とは値が違いますが,結論は同じです。
第10回と同じ方針で考えます。 第10回では,F検定の結果を見てからt検定の方法を選ぶ(2段階の手続き)のは 推奨されないと述べました。分散分析でも同じです。 「Bartlett検定が有意だったからWelch」と機械的に切り替えるのではなく, 群ごとのばらつきの違いを自分で見て判断します。
この授業では,考え方を学ぶために
aov()を主に使います(次回のTukey HSDもaov()の結果を使うため)。ただし実際の研究では,群ごとのSDがこのくらい違うならoneway.test()の結果を併せて報告するのが安全です。 このデータではどちらでも \(p < .001\) で結論は変わりません。
\[\eta^2 = \frac{SS_{群間}}{SS_{全体}}\]
# summary() の表から読み取って計算する
ss_between <- 63.21
ss_total <- 63.21 + 38.96
eta_sq <- ss_between / ss_total
round(eta_sq, 3)## [1] 0.619
出力に表示された丸めた値(63.21・38.96)を打ち込んでいます。 手で確かめるにはこれでよいのですが,桁を落としているので誤差が入ります。
aov()の結果から直接取り出せば,その心配がありません。
## [1] 63.21213 38.95620
## [1] 0.619
同じ \(0.619\) になりました。レポートではこちらの書き方を使ってください。
| η² | 効果の大きさ |
|---|---|
| 0.01 | 小 |
| 0.06 | 中 |
| 0.14以上 | 大 |
iris_summary <- iris |>
group_by(Species) |>
summarise(
mean = mean(Sepal.Length),
sd = sd(Sepal.Length),
n = n()
)
ggplot(iris_summary, aes(x = Species, y = mean, fill = Species)) +
geom_col(width = 0.6) +
geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.15) +
scale_fill_brewer(palette = "Set2") +
labs(title = "種ごとのがく片の長さ(平均±SD)",
x = "種", y = "がく片の長さ (cm)") +
theme_bw() +
theme(legend.position = "none")この形のグラフには注意点が2つあります。
① エラーバーが何を表すかを必ず書く。 上の図は平均 \(\pm\) SD(標準偏差)ですが, 平均 \(\pm\) SE(標準誤差)や95%信頼区間を描くこともあり,見た目では区別できません。 SEはSDを \(\sqrt{n}\) で割った値なので,同じデータでも \(n=50\) なら7分の1ほどの短さになります。 図のタイトルか凡例に「平均±SD」と明記してください。
② 分布の形が見えない。 棒グラフは平均という1つの値に要約してしまうので, 外れ値や分布の偏り,群ごとのばらつきの違いが隠れます。 第3回で学んだ箱ひげ図や,
geom_jitter()で生データを重ねたほうが情報量は多いです。 実際,上の §3 の箱ひげ図と見比べると,virginica のばらつきが大きいことは 箱ひげ図のほうがはっきりわかります。棒グラフは分野の慣習で求められることも多いので「使ってはいけない」わけではありません。 何が見えなくなるかを知って選んでください。
iris データで,3種の
Petal.Length(花びらの長さ)
を一元配置分散分析で比較してください。
group_by() + summarise() で各群の M と SD
を計算するaov() で分散分析を実行し,summary()
で結果を確認する演習1で使ったデータについて,Bartlett 検定で等分散性を確認してください。 あわせて,群ごとのSDも見てください。
bartlett.test(Petal.Length ~ Species, data = iris)
iris |>
group_by(Species) |>
summarise(sd = sd(Petal.Length))このデータでも等分散は成り立ちません
Bartlett 検定は \(p\) < .001 になり,SDも 0.17 / 0.47 / 0.55 と3倍以上の開きがあります。 4節で決めたとおり,検定の結果で手法を機械的に切り替えるのではなく, 結論が変わらないかを自分で確かめます。
Welch法では \(F\) = 1828.09,通常の分散分析では \(F\) = 1180.16 で, どちらも \(p\) < .001。「3種の花びらの長さには差がある」という結論は変わりません。
演習1の summary() の出力から SS_between と SS_total
を読み取り, η² を計算してください。
演習1の結果を,以下のフォーマットで記述してください。
3種のアヤメの花びらの長さを一元配置分散分析によって比較した結果, 種の主効果が有意であった,F(, ) = ____.__, p .001, η² = .。
F値の桁数に注意してください。 本文で使った
Sepal.Lengthは119.26でしたが, 演習のPetal.Lengthは4桁になります。下線の数ではなく,実際の値に合わせて書いてください。
aov() + summary()
で一元配置分散分析を実行