母平均の差についての検定を行う方法を解説していきます。
R言語を用いることで、母平均の差の検定である2標本t検定が簡単に実行できます。
1標本の時と同様に、検定結果などは自分で計算する必要がなく非常に楽です。
また統計学に詳しくなくても、検定結果の各項目について詳しく解説していきます。
今回扱うプログラミングコードは以下よりダウンロードできます。
R言語 2標本t検定 Rスクリプト
テンプレート用として、ぜひぜひダウンロードしてください。
1標本t検定が見たいという方は次のリンクをクリックしてください。
【R言語】関数t.testを用いた母平均の検定 1標本t検定
今回はR言語で、母平均についての検定であるt検定を行う方法を解説していきます。 t検定を実行する関数やその実行例を紹介します。 実行例では、実際にデータセットを用いて母平均の検定を行うだけでなく、検定 ...
続きを見る
また、2標本t検定のより詳しい内容は次の記事を参照してください。
2標本t検定
検定の概要
2標本t検定の概要を以下にまとめました。以下、母分散が等しい場合と等しくない場合の2通りの母平均の差の検定の手順です。
分散が等しいとき
この仮説検定の検定統計量は次で与えられる。
ここに、\(\bar{X}_1\)、\(\bar{X}_2\)、\(S^2\)はそれぞれ
であり
また、有意水準\(\alpha\)の棄却域は次で与えられる。
ここに、\(t_{n, \alpha}\)は自由度\(n\)のt分布の上側\(\alpha\)点である。
分散が異なるとき
この仮説検定の検定統計量は次で与えられる。
ここに、\(\bar{X}_1\)、\(\bar{X}_2\)、\(S_1^2\)、\(S_2^2\)はそれぞれ
であり、\(\nu\)は次で与えられる。
また、有意水準\(\alpha\)の棄却域は次で与えられる。
関数t.test
関数t.testでt検定を実行することができます。以下、t.testとその引数です。
t.test(x, y = NULL, alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, …)
# S3 method for formula
t.test(formula, data, subset, na.action, …)
x | numeric型のベクトル。1つ目のグループの標本(観測値)。 |
y | numeric型のベクトル。2つ目のグループの標本(観測値)。 |
alternative | character型。"two.sided"で両側検定、"greater"で右片側検定、"less"で左片側検定を指定できる。 |
mu | numeric型。帰無仮説の下での母平均、または母平均の差の値。 |
paired | logical型。TRUEのとき対応のある2標本検定を行う。 |
var.equal | logical型。TRUEのとき等分散性を仮定する。 |
conf.level | numeric型。信頼区間を計算する際の信頼度数。 |
formula | formula型。lhs ~ rhsと書き、lhsはnumeric型の変数でrhsは2つのlevelsから成るfactor型の変数。 |
data | データフレーム。検定に用いるデータセット。 |
subset | どの観測値を検定に用いるか指定するためのベクトル。 |
na.action | dataがNAを含む場合どうするか。 |
実行例
次のパッケージggplot2とgridExtraをインストールしましょう。ヒストグラムや箱ひげ図のプロットの際に用います。
1 2 3 4 5 | #2標本検定 install.packages("ggplot2") install.packages("gridExtra") library(ggplot2) library(gridExtra) |
2標本のt検定は、関数t.testを用いることで、行うことができます。データセットとして次のToothLengthを使います。
1 2 3 4 5 6 7 8 9 10 | > #2群の分散が等しいとき > dataset <- ToothGrowth > head(dataset) len supp dose 1 4.2 VC 0.5 2 11.5 VC 0.5 3 7.3 VC 0.5 4 5.8 VC 0.5 5 6.4 VC 0.5 6 10.0 VC 0.5 |
ToothLengthは歯の長さlenや薬の種類suppなどの変数から構成されています。
supp別にlenのヒストグラムや経験密度、箱ひげ図をプロットすると下の画像のようになります。
1 2 3 4 5 6 7 8 9 | histPlot <- ggplot(dataset, aes(x = len, y = ..density.., colour = supp, fill = supp)) + geom_histogram(position = "identity", bins = 25, alpha = 0.75) + geom_density(stat = "density", position = "identity", alpha = 0.75) + xlim(min(dataset$len), max(dataset$len)) boxPlot <- ggplot(dataset, aes(x = len, y = supp, fill = supp)) + geom_boxplot(alpha = 0.75, colour = "gray") grid.arrange(histPlot, boxPlot) |
また、supp別のlenの散布図は以下の通りです。
1 2 3 4 5 6 7 8 9 | ggplot(dataset, aes(x = supp, y = len, color = supp)) + geom_dotplot(binaxis = "y", stackdir = "center") + stat_summary(fun.y = mean, fun.ymin = function(x) mean(x), fun.ymax = function(x) mean(x), geom = "crossbar", width = 0.4) + stat_summary(fun.ymin = function(x) mean(x) - sd(x), fun.ymax = function(x) mean(x) + sd(x), geom = 'errorbar', width = 0.3) |
suppの種類によってlenの異なることが分かります。VCよりもOJの方がlenが大きいことが図から読み取れます。
実際に、t検定で統計的に2つのグループに差があるのかを確かめてみようと思います。
2群の分散が等しいとき
まずは、2つのグループの母分散が等しいと仮定した場合の検定について解説します。
「suppがVCのグループとOJのグループでlenが異なるか」という仮説検定を実行したいと思います。
この仮説検定を実行するには、関数t.testの引数を次のように設定します。
引数var.equal=TRUEにすることで、2グループの等分散性を仮定したt検定を実行することができます。
また、1行目のようにt.testの引数formulaを指定する方法と5行目のように引数xとyを指定する方法の2通りがあります。
1 2 3 4 5 | testResult <- t.test(len ~ supp, data = dataset, var.equal = TRUE) len_VC <- dataset$len[dataset$supp == "VC"] len_OJ <- dataset$len[dataset$supp == "OJ"] t.test(len_VC, len_OJ, var.equal = TRUE) #formulaを用いない場合 |
testResult を参照すると、コンソールにt検定の結果が表示されます。
1 2 3 4 5 6 7 8 9 10 11 12 | > testResult Two Sample t-test data: len by supp t = 1.9153, df = 58, p-value = 0.06039 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.1670064 7.5670064 sample estimates: mean in group OJ mean in group VC 20.66333 16.96333 |
検定統計量やp値、信頼区間などの値が書いてあるのが分かります。
ポイント
これらの検定結果を取得したいときは、次のようにtestResultのうしろに$をつけます。検定統計量、自由度、p値、母平均の信頼区間、標本平均、帰無仮説の下での母平均の差の値などを取得できます。
1 2 3 4 5 6 | stat <- testResult$statistic #検定統計量t値 df <- testResult$parameter #t分布の自由度 pValue <- testResult$p.value #p値 CI <- testResult$conf.int #母平均の差の信頼区間 sampleMean <- testResult$estimate #2群の標本平均 diffMu <- testResult$null.value #帰無仮説の下での母平均の差:μ1-μ2=0 |
また、次のようにデータフレームに格納すると、csvファイル等に保存するときに便利です。
1 2 3 4 | #検定結果をデータフレームにまとめる resultTable_equalVar <- data.frame(stat = stat, d.f. = df, c.i.lower = CI[1], c.i.upper = CI[2], p.value = pValue, row.names = NULL) |
2群の分散が等しくないとき
次に、2群の分散が等しくないときのt検定、すなわちウェルチのt検定の実行方法について紹介します。
ウェルチのt検定を行うには、次のようにt.testの引数var.equalをFLASEにします。関数t.testの引数var.equalはデフォルトでFALSEであるため、var.equal = FALSEは書かなくても大丈夫です。
1 2 | #2群の分散が等しくないとき testResult <- t.test(len ~ supp, data = dataset, var.equal = FALSE) |
testResultを参照すると、次のように検定結果が格納されていることが分かります。
1 2 3 4 5 6 7 8 9 10 11 12 | > testResult Welch Two Sample t-test data: len by supp t = 1.9153, df = 55.309, p-value = 0.06063 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.1710156 7.5710156 sample estimates: mean in group OJ mean in group VC 20.66333 16.96333 |
ポイント
分散が等しい場合と同様に次のようにすることで、このウェルチのt検定の結果を取得することが可能です。これらの検定結果をデータフレームに格納する例を以下に載せます。
1 2 3 4 5 6 7 8 9 10 | stat <- testResult$statistic #検定統計量t値 df <- testResult$parameter #t分布の自由度 pValue <- testResult$p.value #p値 CI <- testResult$conf.int #母平均の差の信頼区間 sampleMean <- testResult$estimate #2群の標本平均 diffMu <- testResult$null.value #帰無仮説の下での母平均の差:μ1-μ2=0 resultTable_notEqualVar <- data.frame(stat = stat, d.f. = df, c.i.lower = CI[1], c.i.upper = CI[2], p.value = pValue, row.names = NULL) |
最後に、検定結果をcsvファイルへ保存する例を紹介します。
分散が等しい場合と異なる場合の2種類のt検定の結果をcsvファイルへ出力したいと思います。
まず、2種類の検定結果が格納されたデータフレームを次のように縦に結合させて、名前を各検定名に変更します。あとは結合させたデータフレームを関数write.csvでcsvファイルへ出力すればおkです。
1 2 3 4 | resultTable <- rbind(resultTable_equalVar, resultTable_notEqualVar) rownames(resultTable) <- c("Student_t", "Welch_t") #検定結果表の行名を検定名にする write.csv(resultTable, "母平均の差の検定.csv") #csvファイルへの出力 |
上記を実行すると次の画像のcsvファイルを作ることができます。
片側仮説検定
補足として、片側検定の実行方法についても紹介します。
片側検定は次のように、t.testの引数alternativeを指定することで行うことができます。
alternative="greater"で右片側検定、alternative = "less"で左片側検定が実行できます。
1 2 3 | #片側仮説検定 testResult <- t.test(len ~ supp, data = dataset, alternative = "greater") t.test(len ~ supp, data = dataset, alternative = "less") |
上のように引数formulaで実行した場合、suppのlevelsの順番がグループの番号に対応します。つまりsuppのlevelsはOJ、VCの順なので、2行目の左片側検定の対立仮説は「OJのグループのlenの母平均<VCのグループのlenの母平均」となることに注意が必要です。
testResultを参照すると次のようになります。
1 2 3 4 5 6 7 8 9 10 11 12 | > testResult Welch Two Sample t-test data: len by supp t = 1.9153, df = 55.309, p-value = 0.03032 alternative hypothesis: true difference in means is greater than 0 95 percent confidence interval: 0.4682687 Inf sample estimates: mean in group OJ mean in group VC 20.66333 16.96333 |
ポイント
OJのグループよりもVCのグループの方がlenは大きいということが分かりました。
以下、検定結果をデータフレームに格納する例です。
1 2 3 4 5 6 7 8 9 10 | stat <- testResult$statistic #検定統計量t値 df <- testResult$parameter #t分布の自由度 pValue <- testResult$p.value #p値 CI <- testResult$conf.int #母平均の差の信頼区間 sampleMean <- testResult$estimate #2群の標本平均 diffMu <- testResult$null.value #帰無仮説の下での母平均の差:μ1-μ2=0 resultTable <- data.frame(stat = stat, d.f. = df, c.i.lower = CI[1], c.i.upper = CI[2], p.value = pValue, row.names = NULL) |
まとめ
R言語で2標本t検定を実行する関数やその実行例をみていきました。
関数t.testを用いることで2標本の母平均の差を実行することができます。
引数var.equalによってStudentのt検定、Welchのt検定を選択することができます。