About connecting the dots.

data science related trivial things

Rでの平均値の差の検定のスニペット

仕事でぱぱっと平均値の差の検定のコードを書かないといけなくて,そういやRでどう書くんだっけとか思って若干調べたので,スニペットがわりにまとめておくことにします.

今回は,サンプルデータとして有名なirisを使います.irisのデータは以下の構成です.

変数名 概要
Sepal.Length センチメートル単位のがく片の長さ
Sepal.Width がく片の幅
Petal.Length 花弁の長さ
Petal.Width 花弁の幅
Species 品種(iris, setosa, versicolor, virginicaの4種類)

チェック項目

基本は以下の3項目をチェック

  • 比較したいのが2群か,3群以上か
  • 対応がある*1
  • 等分散仮定が成り立つか

2群

対応がある

この場合は,対応するサンプル同士の変数の差から,平均と分散を求めることができるので,等分散性の検定を行う必要はありません.単にt検定を行えばOK.下の例のように,花弁の長さと幅を比べた場合,花弁の方が有意に長いことがわかります.

> data(iris)
> t.test(iris$Petal.Length, iris$Petal.Width, paired=TRUE)

	Paired t-test

data:  iris$Petal.Length and iris$Petal.Width
t = 29.7968, df = 149, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 2.388985 2.728348
sample estimates:
mean of the differences 
               2.558667 

対応がない

この場合は等分散性の検定を行う必要があります.setosa種とvirginica種について,花弁の長さの分散が等しいかを調べたところ,(両者の分散の比率が1であるという)帰無仮説が棄却されたため,分散が異なるものと考えられます.

> # サマリの出力
> summary(subset(iris, Species=="setosa")$Petal.Length)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   1.400   1.500   1.462   1.575   1.900 
> summary(subset(iris, Species=="virginica")$Petal.Length)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.500   5.100   5.550   5.552   5.875   6.900 

> # 等分散性の検定
> var.test(subset(iris, Species=="setosa")$Petal.Length, subset(iris, Species=="virginica")$Petal.Length)

	F test to compare two variances

data:  subset(iris, Species == "setosa")$Petal.Length and subset(iris, Species == "virginica")$Petal.Length
F = 0.099, num df = 49, denom df = 49, p-value = 1.875e-13
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.05618945 0.17448557
sample estimates:
ratio of variances 
         0.0990164 
等分散である

上記のように,等分散性が保証されない場合には,等分散性を仮定しないt検定*2を行う必要があります.また等分散性が保証された場合には,普通のt検定を行うことができます.こちらは等分散の場合のt検定になります.

> t.test(subset(iris, Species=="setosa")$Petal.Length, subset(iris, Species=="virginica")$Petal.Length, var.equal=TRUE)

	Two Sample t-test

data:  subset(iris, Species == "setosa")$Petal.Length and subset(iris, Species == "virginica")$Petal.Length
t = -49.9862, df = 98, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -4.252374 -3.927626
sample estimates:
mean of x mean of y 
    1.462     5.552 
等分散ではない

この場合には,var.equal=FALSEにして検定を行えばOKです.

> t.test(subset(iris, Species=="setosa")$Petal.Length, subset(iris, Species=="virginica")$Petal.Length, var.equal=FALSE)

	Welch Two Sample t-test

data:  subset(iris, Species == "setosa")$Petal.Length and subset(iris, Species == "virginica")$Petal.Length
t = -49.9862, df = 58.609, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -4.253749 -3.926251
sample estimates:
mean of x mean of y 
    1.462     5.552 

3群以上

対応がある

必ず等分散性の検定をする必要があります.今度は全品種における花弁の長さの比較を行いたいと思います.対応のある分散分析を行うためには,データをスタック型*3にしておく必要があるので,reshape2パッケージのmelt関数を使います.

今回の例では,(全群の分散の比率が1であるという)帰無仮説が棄却されたため,分散が異なるものと考えられます.

> iris.stack <- melt(iris, id.var="Species")
> head(iris.stack)
  Species     variable value
1  setosa Sepal.Length   5.1
2  setosa Sepal.Length   4.9
3  setosa Sepal.Length   4.7
4  setosa Sepal.Length   4.6
5  setosa Sepal.Length   5.0
6  setosa Sepal.Length   5.4

> bartlett.test(iris.stack$value~iris.stack$variable)

	Bartlett test of homogeneity of variances

data:  iris.stack$value by iris.stack$variable
Bartlett's K-squared = 294.1945, df = 3, p-value < 2.2e-16

通常の分散分析および事後検定をすればOKです.今回の例では,すべての群間で違いがみられましたね.

> summary(aov(iris.stack$value~iris.stack$variable+iris.stack$Species))
                     Df Sum Sq Mean Sq F value Pr(>F)    
iris.stack$variable   3 1656.3   552.1   882.1 <2e-16 ***
iris.stack$Species    2  309.6   154.8   247.3 <2e-16 ***
Residuals           594  371.8     0.6                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> pairwise.t.test(iris.stack$value, iris.stack$variable, p.adjust.method="holm", paired=TRUE)

	Pairwise comparisons using paired t tests 

data:  iris.stack$value and iris.stack$variable 

             Sepal.Length Sepal.Width Petal.Length
Sepal.Width  <2e-16       -           -           
Petal.Length <2e-16       3e-05       -           
Petal.Width  <2e-16       <2e-16      <2e-16      

P value adjustment method: holm 
等分散でない

うーん,ここでは何を使うのかが適切なのかが実は良くわかっていません.いろいろがんばって特殊なノンパラ手法を使うのがいいみたいなんですが,どうなんでしょう...

対応がない

等分散性の検定を行います.

> bartlett.test(iris$Petal.Length~iris$Species)

	Bartlett test of homogeneity of variances

data:  iris$Petal.Length by iris$Species
Bartlett's K-squared = 55.4225, df = 2, p-value = 9.229e-13
等分散である

通常の分散分析および事後検定をすればOKです.

> summary(aov(iris$Petal.Length~iris$Species))
              Df Sum Sq Mean Sq F value Pr(>F)    
iris$Species   2  437.1  218.55    1180 <2e-16 ***
Residuals    147   27.2    0.19                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> pairwise.t.test(iris$Petal.Length, iris$Species, p.adjust.method="bonferroni")

	Pairwise comparisons using t tests with pooled SD 

data:  iris$Petal.Length and iris$Species 

           setosa versicolor
versicolor <2e-16 -         
virginica  <2e-16 <2e-16    

P value adjustment method: bonferroni 
等分散でない

ウェルチの修正分散分析を用いて検定を行います.

> oneway.test(iris.stack$value~iris.stack$variable, var.equal=FALSE)

	One-way analysis of means (not assuming equal variances)

data:  iris.stack$value and iris.stack$variable
F = 860.426, num df = 3.000, denom df = 306.797, p-value < 2.2e-16

> pairwise.t.test(iris$Petal.Length, iris$Species, p.adjust.method="bonferroni")

	Pairwise comparisons using t tests with pooled SD 

data:  iris$Petal.Length and iris$Species 

           setosa versicolor
versicolor <2e-16 -         
virginica  <2e-16 <2e-16    

P value adjustment method: bonferroni 

*1:同一サンプルに対して,複数変数の平均値の比較を行う場合に,対応があるという言い方をします

*2:ウェルチのt検定と呼ばれます

*3:縦にサンプル,横に変数が並んでいるものをアンスタック型,サンプル,変数名,値の3カラムにしてすべてを縦並べにするものをスタック型と呼びます