対応のある2標本t検定を行う方法を解説していきます。
相関のある2つのデータの母平均に差はあるのかという仮説検定の実行例を紹介します。2標本t検定では2つのグループは独立であることを仮定していましたが、対応のある検定によって独立ではないグループ間の平均に関する検定を行うことができます。
2標本t検定について詳しく見たい方は次の記事を参照。
今回扱うプログラミングコードは以下よりダウンロードできます。
テンプレート用として、ぜひぜひダウンロードしてください。
普通の2標本t検定が見たいという方は次のリンクをクリックしてください。
-
【R言語】関数t.testを用いた2標本t検定・母平均の差の検定
母平均の差についての検定を行う方法を解説していきます。 R言語を用いることで、母平均の差の検定である2標本t検定が簡単に実行できます。 1標本の時と同様に、検定結果などは自分で計算する必要がなく非常に ...
続きを見る
2標本t検定
検定の概要
対応のある2標本t検定の概要を以下にまとめました。以下、2つの要素から成る観測値の組が得られたときの母平均の差の検定の手順です。
対応のある2標本t検定
\((x_1, y_1), \ldots, (x_n, y_n)\)はそれぞれ2変量正規分布\(N(\boldsymbol{\mu}, \boldsymbol{\Sigma})\)からの無作為標本とし、\((X_1, Y_1), \ldots, (X_n, Y_n)\)を対応する確率変数の組とする。ここに
このとき、次の「2つの母集団の平均\(\mu_1\)と\(\mu_2\)は等しいか」の仮説検定を考える。
標本に対応がある場合、検定統計量として次を用いる。
また、有意水準\(\alpha\)の棄却域は次で与えらえる。
ここに、\(\bar{X}_D\)と\(S_D^2\)は次で定義される。
通常の2標本t検定との違いとして、2つの母集団は独立であるという仮定が要らないという点があります。
通常の2標本t検定ではそれぞれの標本分散の和から成るカイ二乗統計量を検定統計量に用いていますが、これは2つの母集団が独立であるときにしか成り立ちません。
また、対応のある検定によって、検定に必要なサンプル数も減らすことができるという利点があります。
対応のあるt検定を使う場面として、例えば、薬の投与前と投与後について薬の効果を調べたいときなどがあります。投与前と投与後の2つのグループに同じ被験者が存在するため、2つの母集団は独立であるとはいえないため、対応のある2標本t検定が適しています。
関数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を含む場合どうするか。 |
実行例
対応のある2標本t検定の実行例について見ていきます。
プロットとデータセットのために、次のパッケージggplot2、ggsignif、PairedDataを使います。また、インストールしていない方は、install.packagesを使ってインストールしてください。
1 2 3 4 5 6 | install.packages("PairedData") install.packages("ggsignif") library(ggplot2) library(ggsignif) library(gridExtra) library(PairedData) |
データセットとして次のパッケージPairedDataのAnorexiaを使います。
1 2 | data(Anorexia) dataset <- Anorexia |
次に示すように、Anorexiaは拒食症患者の治療前Priorと治療後Postの体重を変数として持つデータセットとなります。
1 2 3 4 5 6 7 8 | > head(dataset) Prior Post 1 83.8 95.2 2 83.3 94.3 3 86.0 91.5 4 82.5 91.9 5 86.7 100.3 6 79.6 76.7 |
このデータのヒストグラムと箱ひげ図を描くと次のようになります。治療前と治療後で値の変化が分かりやすいように、箱ひげ図には個々のデータの治療前と治療後の体重の遷移も入れました。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | dataset_plot <- data.frame(id = rep(paste("id:", seq_len(nrow(dataset))), 2), weight = c(dataset$Prior, dataset$Post), therapy = as.factor(rep(colnames(dataset), each = nrow(dataset)))) plot_hist <- ggplot(dataset_plot, aes(x = weight, y = ..density.., color = therapy, fill = therapy)) + geom_histogram(position = "identity", bins = 10, alpha = 0.8) + geom_density(stat = "density", position = "identity", alpha = 0.4) plot_box <- ggplot(dataset_plot, aes(x = weight, y = therapy, fill = therapy)) + geom_point(size = 1.5) + geom_boxplot() + geom_segment(data = dataset_plot[dataset_plot$therapy == "Prior", ], aes(x = weight, y = therapy), xend = dataset$Post, yend = rep("Post", nrow(dataset)), color = "#84919E") grid.arrange(plot_hist, plot_box) |
ヒストグラムと箱ひげ図
ヒストグラムと箱ひげ図から治療後の体重は治療前よりも大きいことが分かり、治療によって体重が増えることが推測できます。
また、このデータの散布図を描くと次のように、当たり前ですが治療前と治療後ではある程度相関があるのが確認できます。
1 2 3 4 5 6 7 8 | #変数の相関 r <- cor(dataset$Prior, dataset$Post) ggplot(dataset, aes(x = Prior, y = Post)) + geom_point(size = 3) + stat_ellipse() + annotate(geom = "text", x = -Inf, y = Inf, label = paste0("Sample correlation coefficient:", round(r, 4)), hjust = -0.1 , vjust = 5) |
散布図
同じ拒食症患者のデータが治療前と治療後で含まれているため、治療前と治療後のデータの母集団の分布は独立ではないことがいえます。
2変数間の相関を考慮する必要があるので、「治療前Priorと治療後Postで体重に変化があるか」という仮説検定を対応のある2標本t検定で行います。
対応のある2標本t検定は次のように、関数t.testの引数pairをTRUEにすることで実行することができます。
1 2 | #対応のある2標本t検定 testResult <- t.test(dataset$Prior, dataset$Post, data = dataset, pair = TRUE) |
testResultを参照すると、対応のある2標本t検定の結果がコンソールに出力されます。検定統計量やp値などの検定結果が出力されます。
1 2 3 4 5 6 7 8 9 10 11 12 | > testResult Paired t-test data: dataset$Prior and dataset$Post t = -4.1849, df = 16, p-value = 0.0007003 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -10.94471 -3.58470 sample estimates: mean of the differences -7.264706 |
ポイント
今、「p-value = 0.0007003」と書いてあることから、有意水準0.05の下で治療前Priorと治療後Postの体重に差が無いとは言えない、すなわち治療前Priorと治療後Postで体重に変化があるという結論が得られました。
これらの検定結果を取得したいときは、次のようにtestResultのうしろに$をつけます。検定統計量、自由度、p値、母平均の信頼区間、標本平均、帰無仮説の下での母平均の差の値などを取得できます。
1 2 3 4 5 6 | stat <- testResult$statistic #検定統計量 df <- testResult$parameter #自由度 pValue <- testResult$p.value #p値 ci <- testResult$conf #母平均の差の信頼区間 estimate <- testResult$estimate #母平均の差の推定値 nullValue <- testResult$null.value #帰無仮説の下での母平均の差の値 |
また、次のようにパッケージggsignifを使うことで検定結果を箱ひげ図にまとめることができます。p値を箱ひげ図に加えることができます。
1 2 3 4 5 6 7 8 | #検定結果箱ひげ図 ggplot(dataset_plot, aes(x = therapy, y = weight, fill = therapy)) + geom_boxplot() + geom_signif(comparisons = list(levels(dataset_plot$therapy)), test = "t.test" , test.args = list(alternative = "two.sided", paired = TRUE), map_signif_level = TRUE) + geom_text(x = 1.5, y = max(dataset_plot$weight) + 5, label = paste0("p-value: ", round(pValue, 6)), fontface = 1) + scale_y_continuous(limits = c(min(dataset_plot$weight), max(dataset_plot$weight) + 5)) |
箱ひげ図とp値
また、次のようにデータフレームに格納すると、csvファイル等に保存するときに便利です。下の画像のcsvファイルが保存されます。
1 2 3 4 5 6 | resultTable <- data.frame(stat = stat, d.f. = df, c.i.lower = ci[1], c.i.upper = ci[2], estimate = estimate, p.value = pValue, null.value = nullValue, row.names = NULL) write.csv(resultTable, "対応のある2標本t検定.csv", row.names = FALSE) |
対応のある2標本t検定結果
まとめ
R言語で対応のある2標本t検定を実行する関数やその実行例をみていきました。
関数t.testの引数pairをTRUEにすることで、対応のある2標本t検定が実行することができます。
通常のt検定と同様に、相関のある2変数のデータの母平均の差の検定が可能です。
