Rで検定する

はじめに †

ここでは、Rを使って検定を行います。

準備 †

Rのインストールについては、次のページを見てください。

• MacでRを使う
• WindowsでRを使う

ここでは、標準で使用できるirisデータセットを使います。

data(iris)

このデータセットは、アヤメの種類（Species）をがく片の長さ（Sepal.Length）、幅（Sepal.Width）、花びらの長さ（Petal.Length）、幅（Petal.Width）によって分類する問題です。 長さと幅は連続値、種類はsetosa, versicolor, virginicaのいずれかをとる離散値です。

このデータセットには、setosa, versicolor, virginicaという3種類のアヤメについて、それぞれ50個ずつ、合計150個のデータが含まれています。 ランダムに10個のデータを選択して、見てみましょう。

iris[sort(sample(1:150,10)),]

> iris[sort(sample(1:150,10)),]
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
4            4.6         3.1          1.5         0.2     setosa
22           5.1         3.7          1.5         0.4     setosa
65           5.6         2.9          3.6         1.3 versicolor
97           5.7         2.9          4.2         1.3 versicolor
100          5.7         2.8          4.1         1.3 versicolor
108          7.3         2.9          6.3         1.8  virginica
116          6.4         3.2          5.3         2.3  virginica
122          5.6         2.8          4.9         2.0  virginica
136          7.7         3.0          6.1         2.3  virginica
146          6.7         3.0          5.2         2.3  virginica

論文に使用するグラフはカッコイイ方がいいので、グラフ作成ライブラリーのggplot2でもグラフを作成します。

install.packages("ggplot2")
library(ggplot2)

検定 †

統計的検定では、「AとBの間には差がある」という仮説に対して、「AとBの間には差がない」という反対の仮説を立て、後者の仮説が成り立たないことを示すことによって、前者の仮説が成り立つことを示します。

否定される仮説を帰無仮説、帰無仮説を否定することを棄却するといいます。 帰無仮説が棄却されることによって支持される、帰無仮説の反対の仮説を対立仮説といいます。

平均の差の検定（t検定） †

Petal.Lengthについて、Versicolorの平均は4.906 cm、Virginicaの平均は5.552 cmです。

この二つの平均にの差には、意味があるのでしょうか。 （Virginicaの方がVersicolorよりも大きいと言えるでしょうか。）

これを、t検定という方法によって調べることができます。

t検定では、対立仮説を「二つの平均の差には意味がある」、帰無仮説を「二つの平均の差には意味がない」として、帰無仮説を棄却することによって対立仮説を支持します。

t検定を行うには、t.test関数を使います。

t.test(versicolor$Petal.Length, virginica$Petal.Length)

> t.test(versicolor$Petal.Length, virginica$Petal.Length)

	Welch Two Sample t-test

data:  versicolor$Petal.Length and virginica$Petal.Length
t = -12.604, df = 95.57, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.49549 -1.08851
sample estimates:
mean of x mean of y 
    4.260     5.552 
 

これは、二つの群の分散が同じであることを仮定しない、Welchのt検定といいます。

二つの群の分散が同じであることを仮定する検定はStudentのt検定といい、オプション var.equal=T を付けます。

t.test(versicolor$Petal.Length, virginica$Petal.Length, var.equal=T)

> t.test(versicolor$Petal.Length, virginica$Petal.Length, var.equal=T)

	Two Sample t-test

data:  versicolor$Petal.Length and virginica$Petal.Length
t = -12.604, df = 98, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.495426 -1.088574
sample estimates:
mean of x mean of y 
    4.260     5.552 
 

まず、p値（p-value）を見ます。

p値が、0.1よりも小さいときは「10%有意水準で統計的有意な差がある」、0.05よりも小さいときは「5%有意水準で統計的有意な差がある」、0.01よりも小さいときは「1%有意水準で統計的有意な差がある」といいます。

つぎに、平均の差の95%信頼区間（95 percent confidence interval）を見ます。

この区間がゼロを含むときは、5%有意水準で有意な差があるとはいえません。 「5%有意水準で有意な差がない」ではないことに注意しましょう。

ここで、p値が 1.612e-06 で、0.01よりも小さいため、「VersicolorのPetal.Lengthの平均とVirginicaのPetal.Lengthの平均には1%有意水準で統計的有意な差がある」といえます。

次に、Versicolorを半分ずつに分けて、そのPetal.Lengthの平均を調べてみましょう。

versicolor1 <- versicolor[1:25,]
versicolor2 <- versicolor[26:50,]

mean(versicolor1$Petal.Length)
mean(versicolor2$Petal.Length)

> mean(versicolor1$Petal.Length)
[1] 4.312
> mean(versicolor2$Petal.Length)
[1] 4.208

前半グループの平均は 4.312 cm、後半グループの平均は 4.208 cmで、その差が 0.104 cmありますが、これは意味のある差でしょうか？

t検定を行うと、次のようになります。

t.test(versicolor1$Petal.Length, versicolor2$Petal.Length, var.equal=T)

> t.test(versicolor1$Petal.Length, versicolor2$Petal.Length, var.equal=T)

	Two Sample t-test

data:  versicolor1$Petal.Length and versicolor2$Petal.Length
t = 0.77934, df = 48, p-value = 0.4396
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1643124  0.3723124
sample estimates:
mean of x mean of y 
    4.312     4.208 
 

p値が 0.4396 あり、平均の差の95%信頼区間にも0が含まれています。 したがって、統計的に有意な差があるとはいえません。

比率の差の検定 †

VersicolorとVirginicaについて、Petal.Lengthの値がその平均よりも大きいものがいくつあるか、調べてみましょう。

sum(versicolor$Petal.Length > mean(versicolor$Petal.Length))
sum(virginica$Petal.Length > mean(virginica$Petal.Length))

> sum(versicolor$Petal.Length > mean(versicolor$Petal.Length))
[1] 27
> sum(virginica$Petal.Length > mean(virginica$Petal.Length))
[1] 25

Versicolorでは50個中27個（54%）、Virginicaでは50個中25個（50%）ありました。

この二つの比率の差には、意味があるのでしょうか。 （Versicolorの方がVirginicaよりも比率が高いと言えるでしょうか。）

これを、2標本の比率検定という方法によって調べることができます。

2標本の比率検定では、、対立仮説を「二つの比率の差には意味がある」、帰無仮説を「二つの比率の差には意味がない」として、帰無仮説を棄却することによって対立仮説を支持します。

2標本の比率検定を行うには、prop.test関数を使います。

> prop.test(c(27, 25), c(50, 50))

	2-sample test for equality of proportions with continuity correction

data:  c(27, 25) out of c(50, 50)
X-squared = 0.040064, df = 1, p-value = 0.8414
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.1756826  0.2556826
sample estimates:
prop 1 prop 2 
  0.54   0.50 
 

一つ目の引数は条件を満たす測定値の数のベクトル、二つ目の引数は測定値全体の数のベクトルです。

結果の見方はt検定と同じです。

ここではp値が0.841と0.05よりも大きく、（比率の差の）95%信頼区間がゼロを含みますので、5%有意水準で有意な差があるとはいえません。

独立性の検定 †

今度は、Setosaについて、Sepal.LengthとPetal.Lengthについて、それぞれの値が平均より大きいものがいくつあるか調べてみましょう。

sum(setosa$Petal.Length > mean(setosa$Petal.Length) &
      setosa$Sepal.Length > mean(setosa$Sepal.Length))
sum(setosa$Petal.Length > mean(setosa$Petal.Length) &
      setosa$Sepal.Length <= mean(setosa$Sepal.Length))
sum(setosa$Petal.Length <= mean(setosa$Petal.Length) &
      setosa$Sepal.Length > mean(setosa$Sepal.Length))
sum(setosa$Petal.Length <= mean(setosa$Petal.Length) &
      setosa$Sepal.Length <= mean(setosa$Sepal.Length))

> sum(setosa$Petal.Length > mean(setosa$Petal.Length) &
+       setosa$Sepal.Length > mean(setosa$Sepal.Length))
[1] 15
> sum(setosa$Petal.Length > mean(setosa$Petal.Length) &
+       setosa$Sepal.Length <= mean(setosa$Sepal.Length))
[1] 11
> sum(setosa$Petal.Length <= mean(setosa$Petal.Length) &
+       setosa$Sepal.Length > mean(setosa$Sepal.Length))
[1] 7
> sum(setosa$Petal.Length <= mean(setosa$Petal.Length) &
+       setosa$Sepal.Length <= mean(setosa$Sepal.Length))
[1] 17

これを表にまとめると、次の表のようになります。 このような表を分割表といいます。

そこで、分割表を表す行列を作っておきます。

> m <- matrix(c(15, 11, 7, 17), nrow=2)
> m
     [,1] [,2]
[1,]   15    7
[2,]   11   17

この二つの事象は、独立に生じているでしょうか。 （Sepal.Lengthが大きいこととPetal.Lengthが平均より大きいことは無関係でしょうか。）

これをカイ二乗検定またはFisherの正確検定という方法によって調べることができます。

カイ二乗検定とFisherの正確検定では、対立仮説を「二つの事象は独立でない（関係がある）」、帰無仮説を「二つの事象は独立である（無関係である）」として、帰無仮説を棄却することによって対立仮説を支持します。

基本的には、カイ二乗検定を行いますが、分割表が次のようなときは、Fisherの正確検定を行います。

• 度数が0のセルが存在する、または
• 期待度数が5未満のセルが20%以上存在する

期待度数というのは、それぞれの事象が生じる確率から求められる度数です。 Petal.Lengthが平均より大きい確率は [math]26/50 = 0.52[/math]、Sepal.Lenghが平均より大きいものは22個ありますから、Petal.Lengthが平均より大きく、かつ、Sepal.Lengthが平均より大きいものの期待度数は [math]22 \times 0.52 = 11.44[/math] となります。

上の分割表は、この条件に当てはまらないので、カイ二乗検定を用います。

カイ二乗検定を行うには、chisq.test関数を使います。

> chisq.test(m)

	Pearson's Chi-squared test with Yates' continuity correction

data:  m
X-squared = 3.045, df = 1, p-value = 0.08099
 

分割表を表す行列を作って引数とします。

ここでは、p値が0.081と0.05より大きいので、5%有意水準で独立でない（関係がある）とは言えません。 （帰無仮説が棄却できなかったときは、対立仮説についてはなんとも言えません。「5%有意水準で独立である（関係がない）」とは言えないことに注意しましょう。）

Fisherの正確検定を行うには、fisher.test関数を使います。

> fisher.test(m)

	Fisher's Exact Test for Count Data

data:  m
p-value = 0.05176
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  0.8902322 12.7546580
sample estimates:
odds ratio 
  3.229179 

ここでも、p値が0.051と0.05より大きいので、5%有意水準で独立でない（関係がある）とは言えません。

演習 †

VersicolorのSepal.Lengthの平均とVirginicaのSepal.Lengthの平均の間に統計的有意な差があるかどうかを調べてみよう。

参考文献 †

• 樋口千洋: Rによるバイオインフォマティクスデータ解析 第2版－Bioconductorを用いたゲノムスケールのデータマイニング－, 共立出版 (2011)