Rでの平均値の差の検定のスニペット
仕事でぱぱっと平均値の差の検定のコードを書かないといけなくて,そういやRでどう書くんだっけとか思って若干調べたので,スニペットがわりにまとめておくことにします.
今回は,サンプルデータとして有名なirisを使います.irisのデータは以下の構成です.
変数名 | 概要 |
---|---|
Sepal.Length | センチメートル単位のがく片の長さ |
Sepal.Width | がく片の幅 |
Petal.Length | 花弁の長さ |
Petal.Width | 花弁の幅 |
Species | 品種(iris, setosa, versicolor, virginicaの4種類) |
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