ここでは、2標本の間に平均値の差があるのかを検定することを目的とします。
小難しいことはひとまず置いておいて、Rで実行して見ましょう。
データはRに標準装備されているアヤメ3種の萼長・幅と花弁長・幅のデータ(iris
)を使います。
data(iris)
まずは、どんな風にデータが格納されているかをhead()
関数で見ます。
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
次に、データの特徴を掴むためにsummary()
関数を使って主要な統計量を見てみましょう。
summary(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
##
各種50花分のデータがあり、平均値(Mean)や中央値(Median)などがどの程度かがわかります。
それでは、実際にt検定を行なってみましょう。
今回は、setosaとvirginicaの萼長(Sepal.Length
)の間に差があるのかを検定して行きます。
# それぞれのデータを別々のデータフレームに分けて格納
setosa <- iris %>% filter(Species == "setosa")
virginica <- iris %>% filter(Species == "virginica")
# 平均値を確認しておく
mean(setosa$Sepal.Length); mean(virginica$Sepal.Length)
## [1] 5.006
## [1] 6.588
# t.testでt検定の実行
t.test(setosa$Sepal.Length, virginica$Sepal.Length, var.equal = TRUE)
##
## Two Sample t-test
##
## data: setosa$Sepal.Length and virginica$Sepal.Length
## t = -15.386, df = 98, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.786042 -1.377958
## sample estimates:
## mean of x mean of y
## 5.006 6.588
結果は、p値が0.001以下であるため、
有意水準5%のもとで、setosaとvirginicaの萼長の平均は有意に異なることがわかった。
ただし、t検定は2つの重要な制約がある。
1つ目は、データが正規分布にしたがっているかどうか
2つ目は、両標本の分散が等しいかどうか である。
1つ目のデータが正規分にしたがっているかどうかは、正規性の検定という方法で確かめることができる。
(帰無仮説:本データは正規分布に従う、対立仮説:本データは正規分布に従わない)
それでは、先ほどのデータが正規性を満たすかをシャピロ・ウィルク検定(shapiro.test()
)で確かめる。
shapiro.test(setosa$Sepal.Length)
##
## Shapiro-Wilk normality test
##
## data: setosa$Sepal.Length
## W = 0.9777, p-value = 0.4595
shapiro.test(virginica$Sepal.Length)
##
## Shapiro-Wilk normality test
##
## data: virginica$Sepal.Length
## W = 0.97118, p-value = 0.2583
両標本ともにp値が0.05以上であったため、有意水準5%のもとで帰無仮説は棄却できない。
つまり、正規性を満たすと考えて良いことになる。
2つ目の両標本の分散が等しいかどうかは、等分散性の検定という方法で確かめることができる。
(帰無仮説:両標本の分散は等しい、対立仮説:両標本の分散は異なる)
それでは、先ほどのデータが等分散性を満たすかをF検定(var.test()
)で確かめる。
var.test(setosa$Sepal.Length, virginica$Sepal.Length)
##
## F test to compare two variances
##
## data: setosa$Sepal.Length and virginica$Sepal.Length
## F = 0.30729, num df = 49, denom df = 49, p-value = 6.366e-05
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1743776 0.5414962
## sample estimates:
## ratio of variances
## 0.3072862
p値が0.001以下であったため、有意水準5%のもとで帰無仮説を棄却する。
つまり、両標本の等分散性は満たされない(分散は異なる)ということになる。
ここで、問題となるのが正規生は満たすが、等分散性を満たさないデータにおいて、平均の差を検定したい時である。
(正規性を満たさない場合は、データを変形(変数変換)するか他の方法で検定するか)
そういう場合は、Welchの検定を使う必要がある。
Welchの検定は、等分散性を仮定しないt検定と認識してもらえれば良い。
R上でWelchの検定を実行するのは非常に簡単で、t.test()
のオプションで、var.equal = FALSE
と指定すれば良い。
それでは、先ほどのデータをWelchの検定で検定しみる。
t.test(setosa$Sepal.Length, virginica$Sepal.Length, var.equal = FALSE)
##
## Welch Two Sample t-test
##
## data: setosa$Sepal.Length and virginica$Sepal.Length
## t = -15.386, df = 76.516, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.78676 -1.37724
## sample estimates:
## mean of x mean of y
## 5.006 6.588
結果はほとんど変わらなが、Two Sample t-test
の文言がWelch Two Sample t-test
に、自由度(df)が変わっているのがわかる。
最後に可視化して終わり。
エラーバーは標準偏差(SD)にしておく。
setosa %>%
bind_rows(virginica) %>%
group_by(Species) %>%
summarise(mean = mean(Sepal.Length), sd = sd(Sepal.Length)) %>%
ggplot(aes(x = Species, y = mean)) +
geom_bar(stat = "identity", aes(fill = Species), color = "black", size = 0.8) +
geom_errorbar(aes(ymax = mean + sd, ymin = mean - sd), width = 0.5, size = 0.8) +
scale_y_continuous(breaks = seq(0, 8, 2), limits = c(0, 8)) +
geom_hline(aes(yintercept = 0), size = 0.8) +
theme_bw(base_family = "HiraKakuPro-W3") +
labs(x = "", y = "萼長の平均±SD") +
theme(legend.position = "none",
panel.background = element_rect(fill = "transparent", colour = "black", size = 1.2),
panel.grid = element_blank(),
strip.background = element_blank(),
strip.text = element_text(size = 15, face = "bold"),
axis.title = element_text(size = 15, face = "bold"),
axis.text.y = element_text(size = 15, face = "bold", color = "black"),
axis.text.x = element_text(size = 15, face = "bold", color = "black"))