ここでは、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検定

それでは、実際に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つの制約

ただし、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の検定を使う必要がある。
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"))