t検定・分散分析・線形回帰は
同じ結果を返す
t検定・分散分析・線形回帰の関係?
いままで、t検定、分散分析、そして線形回帰を独立に取り扱ってきた。説明の端々に関連性があることをにおわせてきたが、ここで明示的に、Studentのt検定、2標本の一元配置分散分析、および回帰分析の結果が完全に一致することを紹介しよう。
具体例
以下のようなデータを解析し、結果が一致することを確認する。
------------------------------------------------------
library(plotn)
yA <- c(12, 19, 30, 19, 24, 21, 15, 17, 16, 24,
20, 25, 26, 18, 13, 23, 12, 21, 21, 24,
21, 23, 19, 24, 16, 27, 15, 21, 13, 19)
yB <- c(25, 18, 27, 22, 18, 21, 26, 22, 22, 21,
10, 33, 27, 18, 26, 25, 15, 21, 17, 27,
13, 20, 27, 23, 23, 18, 22, 23, 24, 28,
21, 33, 21, 18, 31, 23, 25, 29, 27, 24)
d <- data.frame(g = c(rep("A",30), rep("B",40)), y = c(yA, yB))
head(d)
## g y
## 1 A 12
## 2 A 19
## 3 A 30
## 4 A 19
## 5 A 24
## 6 A 21
boxplotn(y ~ g, data = d, jitter.method = "swarm")#図1の描画
------------------------------------------------------
図1 データの様子
まず、Studentのt検定を行う。
------------------------------------------------------
fitT <- t.test(d[d$g == "A", "y"], d[d$g == "B", "y"], var.equal = T)
fitT
##
## Two Sample t-test
##
## data: d[d$g == "A", "y"] and d[d$g == "B", "y"]
## t = -2.5013, df = 68, p-value = 0.01479
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -5.2434663 -0.5898671
## sample estimates:
## mean of x mean of y
## 19.93333 22.85000
------------------------------------------------------
P = 0.01479となった。続いて、一元配置分散分析をしてみよう。
------------------------------------------------------
fitA <- aov(y ~ g, data = d)
summary(fitA)
## Df Sum Sq Mean Sq F value Pr(>F)
## g 1 145.8 145.83 6.257 0.0148 *
## Residuals 68 1585.0 23.31
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
------------------------------------------------------
P = 0.0148となった。これは表示の問題で、取り出して確認してみるとぴったり一致していることが確認できる。最後に線形回帰で確認してみる。
------------------------------------------------------
fitL <- lm(y ~ g, data = d)
summary(fitL)
##
## Call:
## lm(formula = y ~ g, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.850 -2.913 0.150 3.837 10.150
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.9333 0.8814 22.614 <2e-16 ***
## gB 2.9167 1.1660 2.501 0.0148 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.828 on 68 degrees of freedom
## Multiple R-squared: 0.08426, Adjusted R-squared: 0.07079
## F-statistic: 6.257 on 1 and 68 DF, p-value: 0.01479
------------------------------------------------------
水準の効果はP = 0.0148となった。下のF値による結果もP = 0.01479と確かに一致している。これは今回のデータでたまたま起こったことでなく、すべて一致する。自分でデータを生成して確認してみるとよいだろう。
このようなことが起こるということは、各解析間で互換性があることを意味している。それを以下に証明する。
Studentのt検定と線形回帰
以下のように証明ができる。
以上のように、最小二乗法によるパラメータ推定とパラメータの分散の計算から、各パラメータのt値を求めることができた。ところで、b1のt値は、水準Aの1標本のt検定の検定統計量となっている。一方で、b2のt値は水準Aと水準Bの2標本のt検定の検定統計量となっている。それゆえに、t検定と線形回帰の結果は一致するのである。というか、内部的には同じ計算をしている、ということである。
Studentのt検定と分散分析
以下のように証明ができる。
以上のように、分散分析の時に求める分散比のF値は、なんとt値の2乗と一致することがわかる。異なる検定統計量を計算しているように見えて、実は同じものだったというわけである。念のために、t値の2乗がF値と一致していることを、具体的計算で確認をしておこう。
------------------------------------------------------
b1e <- coef(summary(fitL))[1,1]
b2e <- coef(summary(fitL))[2,1]
ye <- c(rep(b1e, nrow(d[d$g == "A",])), rep(b1e + b2e, nrow(d[d$g == "B",])))
se <- sum((d$y - ye)^2)/(nrow(d) - 2)
sr <- sum((ye - mean(ye))^2)/1
Fvalue <- sr/se #F-value
Fvalue
## [1] 6.256704
b1e_h <- mean(d[d$g == "A","y"])#b1の推定値を手計算
b2e_h <- mean(d[d$g == "B","y"]) - b1e_h #b2の推定値を手計算
b1e_h
## [1] 19.93333
b2e_h
## [1] 2.916667
b1s_h <- sqrt(var(d[d$g == "A","y"])/nrow(d[d$g == "A",]))#b1の分散を手計算
b2s_h <- sqrt(se*(1/nrow(d[d$g == "A",]) + 1/nrow(d[d$g == "B",])))#b2の分散を手計算
b1s_h
## [1] 0.8454521
b2s_h
## [1] 1.166041
b1t <- b1e_h/b1s_h #b1のt値を手計算
b2t <- b2e_h/b2s_h #b2のt値を手計算
b1t
## [1] 23.57713
b2t
## [1] 2.50134
(b2t)^2 #b2のt値の2乗がFvalueと一致
## [1] 6.256704
------------------------------------------------------
以上のように正規性・等分散を仮定できる2標本の検定においては、Studentのt検定、一元配置分散分析、線形回帰のいずれを用いてもよいことがわかる。ただし、分散が異なるときにはt検定にはWelchのt検定があるし、分散分析や線形回帰は複数の説明変数を扱うことができる。通常はこれらを混同して扱うことはないだろう。しかし、以上の理論的背景は各解析の理解を深めてくれるものに違いない。