【R言語】対応のある2標本t検定　相関のある母平均の差の検定

対応のある2標本t検定を行う方法を解説していきます。

相関のある2つのデータの母平均に差はあるのかという仮説検定の実行例を紹介します。2標本t検定では2つのグループは独立であることを仮定していましたが、対応のある検定によって独立ではないグループ間の平均に関する検定を行うことができます。

2標本t検定について詳しく見たい方は次の記事を参照。

今回扱うプログラミングコードは以下よりダウンロードできます。

R言語　対応のある2標本t検定

1 ファイル 1.51 KB

ダウンロード

テンプレート用として、ぜひぜひダウンロードしてください。

普通の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)$を対応する確率変数の組とする。ここに

\begin{align} \boldsymbol{\mu} &= \begin{pmatrix} \mu_1 \\ \mu_2 \end{pmatrix}, \\ \boldsymbol{\Sigma} &= \begin{pmatrix} \sigma_1^2 & \sigma_{12}\\ \sigma_{21} & \sigma_2^2\end{pmatrix}.\end{align}

このとき、次の「2つの母集団の平均$\mu_1$と$\mu_2$は等しいか」の仮説検定を考える。

\begin{align}&H_0: \delta = \mu_1 -\mu_2 = 0\\&H_1: \delta\neq 0\end{align}

標本に対応がある場合、検定統計量として次を用いる。

\begin{align}\label{eq5}t = \cfrac{\bar{X}_D}{s_D/\sqrt{n}} \sim t_{n- 1}.\tag{5}\end{align}

また、有意水準$\alpha$の棄却域は次で与えらえる。

\begin{align}(-\infty, -t_{\alpha/2, n-1})\cup ( t_{\alpha/2, n-1}, \infty),\end{align}

ここに、$\bar{X}_D$と$S_D^2$は次で定義される。

\begin{align}X_{Di} &= X_{1i}- X_{2i},\\ \bar{X}_D &= \cfrac{1}{n}\sum_{i=1}^n X_{Di}, \\ S_D^2 &= \cfrac{1}{n-1} \sum_{i=1}^n (X_{Di}- \bar{X}_D)^2.\end{align}

通常の2標本t検定との違いとして、2つの母集団は独立であるという仮定が要らないという点があります。

通常の2標本t検定ではそれぞれの標本分散の和から成るカイ二乗統計量を検定統計量に用いていますが、これは2つの母集団が独立であるときにしか成り立ちません。

また、対応のある検定によって、検定に必要なサンプル数も減らすことができるという利点があります。

対応のあるt検定を使う場面として、例えば、薬の投与前と投与後について薬の効果を調べたいときなどがあります。投与前と投与後の2つのグループに同じ被験者が存在するため、2つの母集団は独立であるとはいえないため、対応のある2標本t検定が適しています。

関数t.test

関数t.testでt検定を実行することができます。以下、t.testとその引数です。

# S3 method for default
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の引数
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を使ってインストールしてください。

install.packages("PairedData")
install.packages("ggsignif")
library(ggplot2)
library(ggsignif)
library(gridExtra)
library(PairedData)

install.packages("PairedData")

install.packages("ggsignif")

library(ggplot2)

library(ggsignif)

library(gridExtra)

library(PairedData)

データセットとして次のパッケージPairedDataのAnorexiaを使います。

data(Anorexia)
dataset <- Anorexia

1 2	data(Anorexia) dataset <- Anorexia

次に示すように、Anorexiaは拒食症患者の治療前Priorと治療後Postの体重を変数として持つデータセットとなります。

> 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

> 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

このデータのヒストグラムと箱ひげ図を描くと次のようになります。治療前と治療後で値の変化が分かりやすいように、箱ひげ図には個々のデータの治療前と治療後の体重の遷移も入れました。

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)

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)

ヒストグラムと箱ひげ図

ヒストグラムと箱ひげ図から治療後の体重は治療前よりも大きいことが分かり、治療によって体重が増えることが推測できます。

また、このデータの散布図を描くと次のように、当たり前ですが治療前と治療後ではある程度相関があるのが確認できます。

#変数の相関
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)

#変数の相関

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にすることで実行することができます。

#対応のある2標本t検定
testResult <- t.test(dataset$Prior, dataset$Post, data = dataset, pair = TRUE)

1 2	#対応のある2標本t検定 testResult <- t.test(dataset$Prior, dataset$Post, data = dataset, pair = TRUE)

testResultを参照すると、対応のある2標本t検定の結果がコンソールに出力されます。検定統計量やp値などの検定結果が出力されます。

> 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

> 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値、母平均の信頼区間、標本平均、帰無仮説の下での母平均の差の値などを取得できます。

stat <- testResult$statistic       #検定統計量
df <- testResult$parameter         #自由度
pValue <- testResult$p.value       #p値
ci <- testResult$conf              #母平均の差の信頼区間
estimate <- testResult$estimate    #母平均の差の推定値
nullValue <- testResult$null.value #帰無仮説の下での母平均の差の値

stat <- testResult$statistic #検定統計量

df <- testResult$parameter #自由度

pValue <- testResult$p.value #p値

ci <- testResult$conf #母平均の差の信頼区間

estimate <- testResult$estimate #母平均の差の推定値

nullValue <- testResult$null.value #帰無仮説の下での母平均の差の値

また、次のようにパッケージggsignifを使うことで検定結果を箱ひげ図にまとめることができます。p値を箱ひげ図に加えることができます。

#検定結果箱ひげ図
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))

#検定結果箱ひげ図

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ファイルが保存されます。

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)

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変数のデータの母平均の差の検定が可能です。

usagi-san

統計学とゲームとかをメインに解説していくよ。数式とかプログラミングコードにミスがあったり質問があったりする場合はコメントで受け付けます。すぐに対応します。

2024/04/28

【R言語】F検定テンプレートスクリプト

2024/04/28

【R言語】経験密度関数・経験分布関数のプロット　関数densityとecdfの使い方

2024/04/10

【R言語】ベータ関数とガンマ関数　関数beta, gammaの使い方

usagi-sanの記事をもっと見る

-R言語
-R言語, プログラミング, 統計学, 統計解析

comment コメントをキャンセル

: R言語
【R言語】ベクトルやリストの長さの計算方法　関数lengthの使い方
R言語のベクトルやリストの長さを計算する関数とその使い方について解説していきます。この記事では、関数lengthを用いたベクトルやリストの長さ（要素の数）の計算を見ていきます。また、ベクトルやリス ...

: R言語
R言語入門　データフレーム練習【初心者向け】
こんにちはusagi-sanです。この記事では、R言語のデータフレームに関して重要な知識を紹介します。データフレームの操作はデータ解析や統計解析をするうえで、絶対に必要になってくるので参考になった ...

: R言語
R言語　データ整形練習【初心者向け】
こんにちは、usagi-sanです。今回は、統計解析に用いるデータを整形する方法について解説します。統計解析をする際に、約8割から9割の作業がこのデータ整形に費やされます。マジで面倒くさいです。 ...

: R言語
【R言語】F検定テンプレートスクリプト
R言語のF検定のテンプレートファイルを紹介します。 F検定のスクリプトファイルは以下のgithubまたは、ダウンロードリンクより入手できます。 F検定の実行したいときのテンプレートファイルとして使って ...

: R言語
【R言語】パッケージのインストール　CRAN, GitHub, Bioconductor
R言語のパッケージのインストールの仕方についてまとめました。この記事では、以下のようにCRANだけでなくや他のgit等のリモートリポジトリからパッケージをインストールする手順を解説しています。 CR ...

【R言語】繰り返しの処理　for文・while文・repeat文

【計算フォーム】ガンマ分布の分布関数　下側確率の計算

【R言語】対応のある2標本t検定 相関のある母平均の差の検定

【R言語】関数t.testを用いた2標本t検定・母平均の差の検定

2標本t検定

検定の概要

関数t.test

実行例

まとめ

【R言語】対応のある2標本t検定　相関のある母平均の差の検定