Brief Introduction to R
どんぐり研究所 孫 建強
6 統計検定
R には多くの統計処理用の関数が実装されている。この資料では R の使い方を紹介する目的で、よく使われる t 検定、分散分析（ANOVA）、フィッシャーの正確確率検定およびカイ二乗検定を簡単に紹介する。統計検定を詳しく勉強されたい方は、下記の資料を参考にしてください。
- 生物統計学入門 （成山堂書店）
- 統計学入門 − 基礎統計学 Ⅰ（東京大学出版会）
- 基礎統計（biopapyrus）
6.1 t 検定
t 検定は、母集団から抽出された標本を利用して、その母集団の平均（母平均）を区間推定し、検定を行う方法である。ある母集団の平均が 1 つの値に等しいかどうかを検定する 1 標本 t 検定と 2 つの母集団の平均が等しいかどうかを検定する 2 標本 t 検定がある。3 つ以上の母集団の平均を比べるとき、t 検定を使うことができないため、分散分析やその他の多重比較検定の方法が使われる。なお、t 検定を利用するには、母集団が正規分布に従う必要がある。母集団が正規分布でない場合は、正しい検定結果が得られない。
1 標本 t 検定
1 標本 t 検定は、ある母集団から抽出された標本を利用して、その母集団の平均がある値 k に等しいかどうかを検定する方法である。このとき、帰無仮説として「母平均は k である」を用いる。例えば、ある母集団から 8 標本（3.2, 3.6, 2.9, 2.5, 3.1, 2.7, 3.0, 3.2）を抽出したとき、危険率 5% において、その母平均が2.0 かどうか検定する場合は、次のようにデータ x と比較したい値 2.0 を t.test 関数に代入して実行する。
x <- c(3.2, 3.6, 2.9, 2.5, 3.1, 2.7, 3.0, 3.2)
z <- t.test(x, mu = 2.0)検定結果の変数をそのまま出力させると、t-value や p-value などの値が出力される。特定の値だけを取得したければ、$ を利用する。
z##
## One Sample t-test
##
## data: x
## t = 8.6027, df = 7, p-value = 5.716e-05
## alternative hypothesis: true mean is not equal to 2
## 95 percent confidence interval:
## 2.743258 3.306742
## sample estimates:
## mean of x
## 3.025z$conf.int## [1] 2.743258 3.306742
## attr(,"conf.level")
## [1] 0.95z$p.value## [1] 5.716248e-05この検定結果をみると、標本から推定される母集団は 95% の信頼区間は 2.743258 から 3.306742 までの間にある。この区間に 2.0 が含まれていないため、この母集団の平均は 2.0 ではないと考えられる。また、p-value に着目しても、同様に考えることができる。p-value は検定前に設定した危険率 5% よりも小さいため、帰無仮説（母平均が 2.0 という仮説）を棄却できる。したがって、この母集団の平均は 2.0 ではないと考えられる。
対応なし 2 標本 t 検定
2 つの母集団それぞれから標本を抽出し、それら標本を利用して 2 つの母集団の平均（母平均）が同じかどうか検定する方法として 2 標本 t 検定が使われる。このとき、帰無仮説として「2 つの母平均は同じである」を用いる。実際の検定では、両者の母平均の差の信頼区間にゼロが含まれているかどうかを調べる。
2 つの母集団が互いに独立している場合、両者の母分散が同じかどうかを考慮する必要がある。母分散が同じであれば、t.test 関数に 2 つの母集団から抽出された標本とともに、両者の母分散が同じであることを明示するために var.equal = TRUE オプションを与える。 例えば、x グループと y グループの標本を利用して、危険率 5% において両グループの母平均に差があるかどうかを検定するには、次のように実行する。
x <- c(20.5, 5.3, 12.4, 2.9, 12.3, 6.7, 2.1, 13.1)
y <- c(2.4, 16.1, 21.0, 10.9, 20.6, 25.7, 24.2, 18.0)
z <- t.test(x, y, var.equal = TRUE)
z##
## Two Sample t-test
##
## data: x and y
## t = -2.2788, df = 14, p-value = 0.03888
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -15.4324621 -0.4675379
## sample estimates:
## mean of x mean of y
## 9.4125 17.3625この検定結果を確認すると、両者の母平均の差の 95% 信頼区間は -15.4324621 から -0.4675379 までである。この区間にゼロが含まれていないため、両者の母平均が同じではないと考えられる。危険率 5% のもとで p-value に着目しても、帰無仮説を棄却することができ、両者の母平均に差が認められるとわかる。
2 つのグループの母分散が同じでない場合は、t.test を実行する時に var.equal = FALSE オプションを付ける。このときの t 検定をウェルチ（Welch）の t 検定などとよぶことがある。
x <- c(20.5, 5.3, 12.4, 2.9, 12.3, 6.7, 2.1, 13.1)
y <- c(1.4, 16.1, 31.0, 10.9, 20.6, 15.7, 24.2, 28.0)
z <- t.test(x, y, var.equal = FALSE)
z##
## Welch Two Sample t-test
##
## data: x and y
## t = -2.2376, df = 12.026, p-value = 0.04495
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -17.909469 -0.240531
## sample estimates:
## mean of x mean of y
## 9.4125 18.48752 つのグループの分散が同じではないため、両者の母平均の差の信頼区間が広く見積もられる。
母分散が未知の場合、2 つのグループの母分散が同じかどうかを判断できない。このとき、F 検定を用いて両者の母分散が同じかどうかを検定してから t 検定を行うことがある。また、F 検定を用いずに、初めから両者の母分散が異なるものとしてウェルチの t 検定を行うこともある。
対応あり 2 標本 t 検定
同じ母集団を時系列的に観測していったとき、2 つの時点で同じ標本を抽出し、2 時点における母平均に差があるかどうかを検定する場合、対応あり 2 標本 t 検定を使う。例えば、薬剤投与前と投与後の体重に変化があったかどうか（同じかどうか）を検定する場合に使うことができる。R で対応ありの 2 標本 t 検定を行う場合は t.test 関数に paired = TRUE オプションを付け加える。
before <- c(591, 615, 602, 618, 596)
after <- c(585, 590, 583, 594, 589)
t.test(before, after, paired = TRUE)##
## Paired t-test
##
## data: before and after
## t = 3.9595, df = 4, p-value = 0.01669
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 4.840301 27.559699
## sample estimates:
## mean of the differences
## 16.2なお、この対応あり 2 標本 t 検定は、各標本の投与前と投与後の体重の差がゼロかどうかを検定していることに等しいため、次のようにしても同様な結果が得られる。
diff <- after - before
t.test(diff, mu = 0)##
## One Sample t-test
##
## data: diff
## t = -3.9595, df = 4, p-value = 0.01669
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -27.559699 -4.840301
## sample estimates:
## mean of x
## -16.26.2 多重検定補正
検定を繰り返すと、帰無仮説を間違って棄却される確率が始めに設定した危険率を超える。例えば、独立な検定を 3 回行ったときに、3 つ帰無仮説のうち少なくとも 1 つ以上が間違って棄却される確率は 1 - (1 - 0.05)3 = 0.142625 となる。これが本来設定していた危険率 0.05 よりも約 3 倍の値となっている。このように、危険率 0.05 のもとで検定していたはずなのに、いつの間にか危険率が 0.05 を超えてしまう。この誤りを防ぐために、複数の検定を同時に行ったときに、多重検定補正を行う必要がある。
Bonferroni 補正
n 回の検定を行ったときに、検定全体で危険率を α% に抑えたい場合、それぞれの検定において危険率を (α/n)% として設定する。例えば、危険率 5% のもとで 3 回の検定を行ったとき、Bonferroni 補正を行うと、各検定の危険率を 5/3% に抑えれば良い。このとき、検定全体の危険率は 1 - (1 - 0.05/3)3 = 0.0491713 < 0.05 となる。
Bonferroni 補正の考え方や計算方法が簡単だが、その補正が厳しいことが知られている。例えば、分子生物学の分野で使われるマイクロアレイや RNA-Seq のデータに対して検定を行うとき、一度に数千から数万回の検定を行う必要がある。このとき、Bonferroni 補正を行うと、例えば各回の検定で危険率を 0.05/10000 = 0.000005 未満に抑える必要がある。危険率が低すぎるために、ほとんどの検定で有意差が見られなくなり、重要な遺伝子を拾えなくなる。そのため、検定数が多いとき、Bonferroni 補正を行うのではなく、false discovery rate (FDR) を調整する方法が使われる。
false discovery rate (FDR)
危険率は、帰無仮説が正しいにもかかわらずそれを誤って棄却した過誤を抑えるときに使う。これに対して、検定の過誤を抑える指標として、棄却した帰無仮説の中に真の帰無仮説が含まれる割合を使うこともできる。この割合のことを false discovery rate (FDR) とよぶ。検定の FDR を抑える補正方法として様々に提唱されているが、初期に提唱されている Benjamini と Hochberg による補正方法は、次のように p.adjust 関数にオプションとして method = 'BH' を与えて実行することができる。
p.values <- runif(100)
q.values <- p.adjust(p.values, method = "BH")
head(cbind(p.values, q.values))## p.values q.values
## [1,] 0.02696253 0.4493755
## [2,] 0.60212363 0.9839870
## [3,] 0.87574254 0.9839870
## [4,] 0.92669686 0.9839870
## [5,] 0.53641651 0.9839870
## [6,] 0.78231390 0.9839870FDR を 0.10 未満に抑えたければ、上で計算した q-value が 0.10 未満の検定を有意として扱えばよい。
6.3 分散分析（ANOVA）
3 群間以上の母平均を比較して、互いに差があるかどうかを検定するには分散分析（ANOVA）を用いる。なお、分散分析を用いるとき、データが正規分布に従う必要がある。データが正規分布に従わない場合、正しい結果が得られない。
R で分散分析を行うとき aov 関数または anova 関数を使う。aov 関数は、変数を徐々に追加していったときに平方和がどれぐらい増えるのかを計算する、というタイプ I の平行和に基づいて分散分析を行っている。この aov 関数について見ていく。
ここで R のデータセット PlantGrowth を使って分散分析する例を示す。PlantGrowth データセットは 2 列からなるデータフレームで、3 つの処理群（group 列）における植物の乾燥重量（weight 列）を測定したデータが記録されている。
data(PlantGrowth)
head(PlantGrowth)## weight group
## 1 4.17 ctrl
## 2 5.58 ctrl
## 3 5.18 ctrl
## 4 6.11 ctrl
## 5 4.50 ctrl
## 6 4.61 ctrlこのデータセットを使って、各処理で育てられた植物の乾燥重量に差があるかどうかを分散分析で検定してみる。各処理によって乾燥重量に変化があったかどうかを検定したいので、因子が処理（group）で、値が乾燥重量（weight）となる。このとき、aov で分散分析を行うとき、因子と値の間の関係を weight ~ group のように記述する。ただし、因子は数値ではないので、明示的に factor(group) とした方が間違いが起きにくい。例えば、群名が文字ではなく 1, 2, 3 などのように記載されたときに factor(group) で記述しなかった場合、回帰分析のような計算が行われる。間違って回帰分析が行われないようにするためには、factor(group) を記述するか、予めデータフレームの group 列の値を因子型に変換しなければならない。
z <- aov(weight ~ factor(group), data = PlantGrowth)
summary(z)## Df Sum Sq Mean Sq F value Pr(>F)
## factor(group) 2 3.766 1.8832 4.846 0.0159 *
## Residuals 27 10.492 0.3886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1aov 関数の出力結果を確認すると、factor(group) の p-value が 0.159 であることがわかる。危険率 5% のもとで、乾燥重量の平均値が処理間で有意に異なっていると考えられる。ただし、分散分析の場合は、どの群の平均値が大きくて、どの群の平均値が小さいのかまでを検定することはできない。母平均の大小関係も検定したい場合は、Tukey 検定、Dunnett 検定や Williams 検定など利用する。
次に anova 関数の使い方を見ていく。anova 関数では、与えられた回帰モデルを基準として、そのモデルから説明変数（因子）を取り除いたときに、平方和が減るかどうかを検定している。タイプ II の分散分析とも呼ばれている。そのため、anova 関数を使用するときは、すべての説明変数を使って構築した回帰モデルを与えることで、anova 関数の内部で説明変数を一つずつ取り除きながら平方和の減りを計算し、取り除かれた説明変数が、目的変数を説明するために有意であるかどうかを検定している。PlantGrowth データセットの場合、説明変数が処理（group）だけであるので、group を取り除いたときの変化を見ることになる。
anova(lm(weight ~ factor(group), data = PlantGrowth))## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(group) 2 3.7663 1.8832 4.8461 0.01591 *
## Residuals 27 10.4921 0.3886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1次に、複数の因子を持つデータセットを分散分析で解析する例を示す。npk データセットを使う。このデータセットは、各圃場実験区（block）において、肥料 N、P、K の有無によって収量がどのぐれいになるのかを記録したデータである。
data(npk)
head(npk)## block N P K yield
## 1 1 0 1 1 49.5
## 2 1 1 1 0 62.8
## 3 1 0 0 0 46.8
## 4 1 1 0 1 57.0
## 5 2 1 0 0 59.8
## 6 2 1 1 1 58.5N、P、K の有無によって収量が変化するかどうかを aov 関数と anova 関数で検定していく。
z <- aov(yield ~ factor(N) + factor(P) + factor(K), data = npk)
summary(z)## Df Sum Sq Mean Sq F value Pr(>F)
## factor(N) 1 189.3 189.28 6.488 0.0192 *
## factor(P) 1 8.4 8.40 0.288 0.5974
## factor(K) 1 95.2 95.20 3.263 0.0859 .
## Residuals 20 583.5 29.17
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(yield ~ factor(N) + factor(P) + factor(K), data = npk))## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(N) 1 189.28 189.282 6.4880 0.01919 *
## factor(P) 1 8.40 8.402 0.2880 0.59743
## factor(K) 1 95.20 95.202 3.2632 0.08592 .
## Residuals 20 583.48 29.174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1これらの関数の出力結果を確認すると、危険率 5% のもとで、N の有無が収量に影響を与えていることがわかる。これに対して、P および K は収量に影響を与えていないと考えられる。
6.4 マン・ホイットニーの U 検定
マン・ホイットニーの U 検定（Mann–Whitney U test）は、分布の形が同じ 2 つの母集団の値に順序関係があるときに、 それらの分布の位置が同じかどうかを検定するノンパラメトリック方法の一つである。対応のない 2 標本 t 検定に対応したノンパラメトリックな検定方法である。なお、マン・ホイットニーの U 検定と実質的に同じ検定としてウィルコクソンの順位和検定と呼ばれるものがある。R では、ウィルコクソンの順位和検定にちなんで wilcox.test 関数が用意されている。
次は、分布の形は同じだが、一方の分布の値が他方よりも大きい母集団から抽出した標本に対して検定を行う例である。危険率 5% のもとでマン・ホイットニーの U 検定の結果を確認すると、確かに両者の分布が異なっていることがわかる。
x <- c(rnorm(10, 60, 10), rnorm(20, 60, 20))
y <- c(rnorm(10, 90, 10), rnorm(20, 90, 20))
par(mfrow = c(1, 2))
hist(x)
hist(y)
wilcox.test(x, y)##
## Wilcoxon rank sum exact test
##
## data: x and y
## W = 97, p-value = 1.537e-08
## alternative hypothesis: true location shift is not equal to 0マン・ホイットニーの U 検定（ウィルコクソンの順位和検定）と似たノンパラメトリックな検定方法に、ウィルコクソンの符号付順位検定とよばれる検定方法がある。ウィルコクソンの符号付順位検定は、対応のある 2 標本 t 検定に対応したノンパラメトリックな検定方法である。R でウィルコクソンの符号付順位検定を行うには、coin パッケージに実装されている wilcox_test 関数を使用する。
library(coin)
x <- c(rnorm(10, 60, 10), rnorm(20, 60, 20))
y <- x + rnorm(length(x), 10, 6)
wilcoxsign_test(x ~ y, distribution = 'exact')##
## Exact Wilcoxon-Pratt Signed-Rank Test
##
## data: y by x (pos, neg)
## stratified by block
## Z = -4.741, p-value = 5.588e-09
## alternative hypothesis: true mu is not equal to 06.5 フィッシャーの正確確率検定
フィッシャー（Fisher）の正確確率検定は割合に関する検定方法の一つである。例えば、あるシーズンにおいて、インフルエンザ予防接種を受けなかった 532 人を追跡調査した際に、238 人が感染し、294 人が感染しなかった。これに対して、予防接種を受けた 386 人のうち 121 人が感染し、265 人が感染しなかった。このデータを利用して、インフルエンザの予防接種に効果があるかどうかを検定したい。
この場合、予防接種を受けたグループとそうでないグループの感染者・非感染者の割合が同じであることを帰無仮説としてフィッシャーの正確確率検定を用いる。
dat <- matrix(c(238, 121, 294, 265), ncol = 2)
rownames(dat) <- c("uninoculated", "inoculation")
colnames(dat) <- c("infection", "uninfected")
dat## infection uninfected
## uninoculated 238 294
## inoculation 121 265fisher.test(dat)##
## Fisher's Exact Test for Count Data
##
## data: dat
## p-value = 3.96e-05
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.335179 2.357439
## sample estimates:
## odds ratio
## 1.771845予防接種を受けたグループの非感染者の割合は、そうでないグループの非感染者の割合に比べて大きく、両者の比率（オッズ比）の 95% 信頼区間は 1.335179 から 2.357439 までである。この 95% 信頼区間に 1.0 が含まれていないため、予防接種を受けたグループとそうでないグループの割合が同じではないと考えられる。
6.6 カイ二乗検定
カイ二乗検定は適合度検定および独立性検定の目的で使われている。適合度検定では、観測された度数分布が理論分布（または母集団の分布）と同じかどうかを検定することになる。例えば、あるサイコロが公平であるかどうかを判断するときに、カイ二乗検定（適合度検定）が使われる。
サイコロを 60 回振ったとき、1 から 6 の目が出る回数がそれぞれ 8, 11, 10, 9, 12, 10 であったとする。このとき、もしこのサイコロが公平であれば、1 から 6 の目が出る確率は、1/6, 1/6, 1/6, 1/6, 1/6, 1/6 となるはずである。このことを利用して、観測数が実際の理論値（確率）に従っているかどうかを検定する。
x <- c(8, 11, 10, 9, 12, 10)
chisq.test(x, p = rep(1/6, length = 6))##
## Chi-squared test for given probabilities
##
## data: x
## X-squared = 1, df = 5, p-value = 0.9626次に、独立性検定について考える。フィッシャーの正確確率検定の例と同じ設定で考えてみる。あるシーズンにおいて、インフルエンザ予防接種を受けなかった 532 人を追跡調査した際に、238 人が感染し、294 人が感染しなかった。これに対して、予防接種を受けた 386 人のうち 121 人が感染し、265 人が感染しなかった。このデータを利用して、インフルエンザの予防接種に効果があるかどうかを検定してみる。
このデータでは、調査対象者数は 532 + 386 = 918 人であり、そのうち 238 + 121 = 359 人がインフルエンザに感染した。感染率は 359/918 = 0.391 である。もし、予防接種を受けても効果がないと仮定した場合、予防接種を受けたグループも、そうでないグループも、感染率は理論上 0.391 になるはずである。カイ二乗検定は、両グループの感染率が理論値（0.391）に等しいかどうかを検定する。
dat <- matrix(c(238, 121, 294, 265), ncol = 2)
rownames(dat) <- c("uninoculated", "inoculation")
colnames(dat) <- c("infection", "uninfected")
dat## infection uninfected
## uninoculated 238 294
## inoculation 121 265chisq.test(dat)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: dat
## X-squared = 16.284, df = 1, p-value = 5.453e-05