今回はR言語で、母平均についての検定であるt検定を行う方法を解説していきます。
t検定を実行する関数やその実行例を紹介します。
実行例では、実際にデータセットを用いて母平均の検定を行うだけでなく、検定結果を取得しcsvファイルへ出力する例もみていきます。
今回扱うプログラミングコードは以下よりダウンロードできます。
テンプレート用として、ぜひぜひダウンロードしてください。
2標本検定が見たい方は、次の記事をお勧めします。
-
【R言語】関数t.testを用いた2標本t検定・母平均の差の検定
母平均の差についての検定を行う方法を解説していきます。 R言語を用いることで、母平均の差の検定である2標本t検定が簡単に実行できます。 1標本の時と同様に、検定結果などは自分で計算する必要がなく非常に ...
続きを見る
また、t検定のより詳しい内容は次の記事を参照してください。
-
【統計学】t検定 母平均の検定・母平均の差の検定
ここでは、統計学の仮説検定において重要なStudentのt検定について解説する。 母分布である正規分布のパラメータによって様々なt検定の手法が提案されている。 その中でもよく使われるt検定についてみて ...
続きを見る
ポチップ
ポチップt検定
1標本t検定の概要とt検定を実行する関数を紹介します。
t検定の概要
t検定の検定統計量や検定手順は以下の通りです。
この後紹介する関数で以下の1標本t検定を実行することが可能です。
1標本t検定
\(x_1, \ldots, x_n\)をそれぞれ\(N(\mu, \sigma^2)\)からの独立同一な標本であるとし、\(X_1, \ldots, X_n\)をその確率変数とする。このとき、次の「ある母集団の平均\(\mu\)がある特定の値\(\mu_0\)であるか」の仮説検定をしたい場合、1標本検定を適用する。
検定統計量として次を用いる。
また、有意水準\(\alpha\)の棄却域は次で与えらえる。
ここに、\(\bar{X}\)と\(S^2\)は次で定義され、\(t_{\alpha/2, n-1}\)は自由度\(n-1\)のt分布の上側\(\alpha/2\)の確率点である。
関数t.test
1標本の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, …)
t.testの引数は以下の通りです。引数を見てもらうとわかる通り、1標本だけでなく2標本や対応のある検定にも行えます。
| 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を用います。まだインストールしてない場合は、install.packagesでインストールしておいてください。
1 2 3 4 | install.packages("ggplot2") install.packages("gridExtra") library(ggplot2) library(gridExtra) |
t検定を行うために、まず次のデータセットtreesを用意します。
1 | dataset <- trees |
treesは次のように、Girth（外周の長さ）、Height（高さ）、Volume（重さ）の3つの変数から構成されています。
1 2 3 4 5 6 7 8 | head(dataset) Girth Height Volume 1 8.3 70 10.3 2 8.6 65 10.3 3 8.8 63 10.2 4 10.5 72 16.4 5 10.7 81 18.8 6 10.8 83 19.7 |
今、Heightに関心があり、「Heightの平均は80であるか」という仮説検定を実行してみたいと思います。
検定を実行する前に、Heightがどのように分布しているのか確認します。次を実行すると、下の画像のHeightのヒストグラムと箱ひげ図をプロットすることができます。
1 2 3 4 5 6 7 8 9 | plot_hist <- ggplot(dataset, aes(x = Height, y = ..density..)) + geom_histogram(position = "identity", bins = 7, alpha = 0.8) + geom_density(stat = "density", position = "identity", fill = "black", alpha = 0.5) plot_box <- ggplot(dataset, aes(x = Height, y = 1)) + geom_point(size = 1.5) + geom_boxplot(alpha = 0.5) + labs(y = "") grid.arrange(plot_hist, plot_box) |
Heightは左右対称に分布しており、正規性をもつことが確認できました。
正規性をもつのでt検定を用いた母平均の検定を行うことができます。
「Heightの平均は80であるか」という仮説検定は、t.testを用いて次のように実行できます。
1 2 | mu0 <- 80 testResult <- t.test(dataset$Height, mu = mu0) |
testResultを参照すると、検定統計量やp値などが出力されます。
1 2 3 4 5 6 7 8 9 10 11 12 | > testResult One Sample t-test data: dataset$Height t = -3.4952, df = 30, p-value = 0.001496 alternative hypothesis: true mean is not equal to 80 95 percent confidence interval: 73.6628 78.3372 sample estimates: mean of x 76 |
ポイント
今、「p-value = 0.001496」より、帰無仮説は棄却されてHeightの平均は80であるとはいえないという結論を得ました。Heightの平均は80ではないことが分かりました。
この検定の検定結果を取得するには、次のように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 #母平均の信頼区間 estimate <- testResult$estimate #標本平均 nullValue <- testResult$null.value #帰無仮説の下での平均値の値 |
また、これらの検定結果をcsvファイル等に保存したい場合は、次のようにデータフレームに格納すると便利です。
write.csvなどでデータフレームを出力することで、csvファイル等の外部のファイルに保存することができます。
1 2 3 4 5 | resultTable <- data.frame(stat = stat, d.f. = df, p.value = pValue, estimate = estimate, c.i.lower = CI[1], c.i.upper = CI[2], null.value = nullValue, row.names = NULL) write.csv(resultTable, "1標本t検定結果.csv", row.names = FALSE) #csvファイルへの出力 |
片側検定の実行例も紹介します。
関数t.testで紹介したように、引数alternativeを"greater"もしくは"less"に指定することで右片側検定、左片側検定を実行することができます。
「Heightの平均は80より大きくないか」という仮説検定を実行したい場合は、次のように引数alternativeを"greater"に設定するだけです。
1 2 3 4 5 6 7 8 9 10 11 12 | t.test(dataset$Height, mu = mu0, alternative = "greater") #右側仮説検定 One Sample t-test data: dataset$Height t = -3.4952, df = 30, p-value = 0.9993 alternative hypothesis: true mean is greater than 80 95 percent confidence interval: 74.05764 Inf sample estimates: mean of x 76 |
ポイント
「p-value = 0.9993」より、帰無仮説は採択されて、Heightの平均は80より大きいとはいえないことが分かりました。
ポチップ検定結果の解釈
t検定の概要をプロットし、検定の棄却域や信頼区間がどのようになっているのか簡単にみていきます。
次を実行すると、下の図のt分布とt検定の棄却域、検定統計量の値がプロットされます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | #t分布 x <- seq(-4, 4, 0.001) dataset_plot <- data.frame(x, density = dt(x, nrow(dataset) - 1)) ggplot() + geom_line(data = dataset_plot, aes(x = x, y = density), size = 1) + geom_segment(data = data.frame(x = c(qt(c(0.025, 0.975), nrow(dataset) - 1), stat), y = c(dt(qt(c(0.025, 0.975), nrow(dataset) - 1), nrow(dataset) - 1), Inf), name = as.factor(c("reject rejon", "reject rejon", "statistics"))), aes(x = x, y = y, color = name), xend = c(qt(c(0.025, 0.975), nrow(dataset) - 1), stat), yend = rep(0, 3), size = 1) + geom_text(data = data.frame(x = c(qt(c(0.025, 0.975), nrow(dataset) - 1), stat), name = as.factor(c("reject rejon", "reject rejon", "statistics"))), aes(x = x, y = -0.01, color = name, label = round(x, 2))) + geom_ribbon(data = dataset_plot[x < qt(0.025, nrow(dataset) - 1), ], aes(x = x, y = density, ymin = 0, ymax = density), fill = "#F8766D", alpha = 0.5) + geom_ribbon(data = dataset_plot[x > qt(0.975, nrow(dataset) - 1), ], aes(x = x, y = density, ymin = 0, ymax = density), fill = "#F8766D", alpha = 0.5) + scale_y_continuous(limits = c(-0.01, 0.4)) + labs(color = "") |
図に示すように、検定統計量が赤い領域である棄却域に含まれているため、帰無仮説が棄却されて「Heightの平均は80であるとはいえない」という結論を得るのが分かります。
また、次を実行すると下の画像の母平均の信頼区間を描くことができます。
1 2 3 4 5 6 7 8 9 10 | #母平均の信頼区間 ggplot(dataset, aes(x = 1, y = Height)) + geom_dotplot(binaxis = "y", stackdir = "center") + stat_summary(aes(x = 1, y = Height), fun.y = mean, fun.ymin = function(x) mean(x), fun.ymax = function(x) mean(x), geom = "crossbar", width = 0.3, color = "#F8766D") + stat_summary(aes(x = 1, y = Height), fun.ymin = function(x) CI[1], fun.ymax = function(x) CI[2], geom = 'errorbar', width = 0.25, color = "#619CFF", size = 1) + geom_text(data = data.frame(x = c(estimate, CI), name = factor(c("sample mean", "c.i.", "c.i."), levels = c("sample mean", "c.i."))), aes(x = 0.7, y = x, color = name, label = paste0(c("sample mean", "c.i.lower", "c.i.upper"), ": ", round(x, 3)))) |
この図から、母平均の95%信頼区間は[73, 78]であり、母平均はこの標本の下で73から78の範囲にあることが確認できます。
帰無仮説の下での母平均の値80は、この信頼区間に含まれていないことが、先ほどの検定統計量が棄却域に含まれていることに対応しています。
このように、検定統計量や信頼区間に関するグラフを描くことで、検定結果が視覚的に理解しやすくなります。
まとめ
R言語で1標本t検定を実行する関数やその実行例を見ていきました。
関数t.testを用いることで簡単に母平均の検定を行うことができます。
