R言語でフィッシャーのZ変換を用いた相関係数の検定を行う方法を紹介します。
この記事では、正規近似であるフィッシャーのZ変換による相関係数の検定と相関係数の差の検定を解説していきます。
相関係数の検定では、無相関でなくある特定の値についての相関係数の検定方法について見ていきます。
この記事で紹介するプログラミングコードは以下からダウンロードできます。
t分布に基づいた無相関性の検定については、相関係数の検定の詳細については次の記事を参照してください。
-
【R言語】相関係数の検定・無相関性の検定 関数cor.testの使い方
R言語で相関係数の検定を行う方法を紹介します。 この記事では、t分布に基づく無相関性の検定を実行する関数やその使い方について見ていきます。 フィッシャーのz変換を用いた相関係数の検定、差の検定について ...
続きを見る
ポチップ
相関係数の検定
フィッシャーのZ変換を用いた相関係数の検定の関数をいくつか紹介します。
関数fisherZ
Fisherのz変換を行う関数であるFisherZやそれに関連する関数を紹介します。
FisherZInv(z)
CorCI(rho, n, conf.level = 0.95, alternative = c("two.sided", "less", "greater"))
上で紹介した関数の引数は次の通りです。
| rho | Pearsonの積率相関係数。 |
| z | Fisherのz変換後の値。 |
| n | 信頼区間を計算する際のサンプルサイズ。 |
| alternative | character型。"two.sided"で両側仮説、"greater"で右片側仮説。"less"で左片側仮説を指定する。 |
| conf.level | 信頼区間の信頼水準。 |
パッケージDescToolsをインストールすることで関数FisherZを使うことができます。
FisherZを用いることで、実数\(r,\ 0 < r <1\)を次のFisherのz変換によって変換することができます。
\begin{align}z = \cfrac{1}{2} \log \cfrac{1 + r}{1 - r}.\end{align}
このFisherのz変換の検定の関数はパッケージにありませんが、正規近似（デルタ法）によって次のように、相関係数の検定、相関係数の差の検定を行うことができます。
相関係数の検定（Fisherのz変換）
\((x_1, y_1), (x_2, y_2)\ldots, (x_n, y_n)\)の2組から成る標本が与えられているとし、ピアソンの標本相関係数を\(r\)とする。また、\(X\)と\(Y\)の相関係数を\(\rho\)とする。このとき次の仮説検定を考える。\begin{align}&H_0:\ \rho = \rho_0\\ &H_1:\ \rho \neq \rho_0\end{align}\(n\)が十分に大きいとき、この仮説検定の検定統計量は次で与えられる。\begin{align}\sqrt{n - 3}(Z - \zeta) \sim N(0, 1), \end{align}ここに\begin{align}Z &= \cfrac{1}{2} \log \cfrac{1 + R}{1 - R} ,\\ \zeta &= \cfrac{1}{2} \log \cfrac{1 +\rho_0}{1 - \rho_0}.\end{align}この検定の有意水準\(\alpha\)の棄却域は次で与えられる。\begin{align}(-\infty, - z(\alpha / 2)) \cup (z(\alpha / 2), \infty),\end{align}ここに、\(z(\alpha)\)は標準正規分布の上側\(\alpha\)点である。
また、相関係数の差の検定は以下の通りです。
相関係数の差の検定（Fisherのz変換）
\((x_{11}, y_{11}), \ldots, (x_{1n_1}, y_{1n_1})\)と\((x_{21}, y_{21}), \ldots, (x_{2n_2}, y_{2n_2})\)の2群から2組の標本が与えられているとし、それぞれの群のピアソンの標本相関係数を\(r_1\)、\(r_2\)とする。また、それぞれの群の\(X\)と\(Y\)の母相関係数を\(\rho_1\)、\(\rho_2\)とする。このとき次の仮説検定を考える。\begin{align}&H_0:\ \rho_1 = \rho_2\\ &H_1:\ \rho_1 \neq \rho_2\end{align}\(n_1\)、\(n_2\)が十分に大きいとき、この仮説検定の検定統計量は次で与えられる。\begin{align}(Z_1 - Z_2) \left/ \sqrt{\cfrac{1}{\sqrt{n_1 - 3}} + \cfrac{1}{\sqrt{n_2 - 3}}}\right. \sim N(0, 1), \end{align}ここに\begin{align}Z_i &= \cfrac{1}{2} \log \cfrac{1 + R_i}{1 - R_i},\quad i =1,2 .\end{align}この検定の有意水準\(\alpha\)の棄却域は次で与えられる。\begin{align}(-\infty, - z(\alpha / 2)) \cup (z(\alpha / 2), \infty),\end{align}ここに、\(z(\alpha)\)は標準正規分布の上側\(\alpha\)点である。
実行例
上で紹介したFisherZを用いて相関係数の検定を行う例をみていきます。
相関係数の検定
Fisherのz変換を用いた相関係数の検定方法について紹介します。
次のパッケージDescToolsをインストールすることでFisherのz変換に関する関数をインストールすることができますが、cor.testのような検定を行ってくれる関数は実装されていないため、自分で実装する必要があります。（検定方法はデルタ法に基づく正規近似を用いているので、FisherZのRDocumentationに載っているt分布による方法とは異なります。）
早速、パッケージDescToolsをインストールしましょう。
1 2 | install.packages("DescTools") library(DescTools) |
次のデータセットtreesを例に相関係数の検定を行いたいと思います。
1 2 3 4 5 6 7 8 9 | > dataset <- trees > 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 |
GirthとHeightの散布図は下の図のようになります。
1 2 3 4 | library(ggplot2) ggplot(dataset, aes(x = Girth, y = Height)) + geom_point(size = 2) + stat_ellipse() |
散布図
この散布図から、GirthとHeightに相関があるのが見て取れます。
「GirthとHeightの相関が高い、すなわちρ=0.8である」という仮説検定を行っていきます。
以下のように帰無仮説の下での相関係数の値と信頼水準を与え、標本数、標本相関係数を計算しておきます。
1 2 3 | conf.level <- 0.95 #信頼水準 rho0 <- 0.8 #帰無仮説の下での値 n <- nrow(dataset) #標本数<br />r <- cor(dataset$Girth, dataset$Height) #標本相関係数 |
Fisherのz変換により、この仮説検定のp値と信頼区間は次のように計算できます。
1 2 3 4 5 6 | z <- FisherZ(r) zeta <- FisherZ(rho0) stat <- sqrt(nrow(dataset) - 3) * (z - zeta) pValue <- 2 * (1 - pnorm(abs(stat))) ci <- CorCI(r, n, conf.level) |
検定結果は以下のようになります。
1 2 3 4 5 | > resultTable <- data.frame(stat = stat, p.value = pValue, + c.i.lower = ci["lwr.ci"], c.i.upper = ci["upr.ci"], row.names = NULL) > resultTable stat p.value c.i.lower c.i.upper 1 -2.768825 0.005625882 0.2021327 0.7378538 |
ポイント
p値=0.005625882 < 0.05より、優位水準0.05の下で「GirthとHeightの相関係数は0.8である（0.8ではないことはいえない）」ことが分かりました。
検定結果を保存するには次のように検定結果を保存したデータフレームをwrite.csv等で出力します。
1 | write.csv(resultTable, "相関係数の検定（フィッシャーのZ変換）.csv", row.names = FALSE) |
相関係数の検定結果
相関係数の差の検定
次に相関係数の差の検定の実行例を紹介します。
データセットとして次のSpeciesがsetosaとversicolorのirisを用います。
1 | dataset <- iris[iris$Species %in% c("setosa", "versicolor"), ] |
SpeciesがsetosaとversicolorのSepal.LengthとPetal.Lengthの散布図を描くと次のようになります。
1 2 3 | ggplot(dataset, aes(x = Sepal.Length, y = Petal.Width, color = Species)) + geom_point(size = 2) + stat_ellipse() |
2群の散布図
この散布図からsetosaとversicolorによってSepal.LengthとPetal.Lengthの相関が異なることがわかります。 setosaのときよりもversicolorのほうが相関が強いことが分かります。
関数FisherZで相関係数の差の検定を行うには、次のように検定統計量statと、p値pValueを計算します。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | sepalLength_setosa <- dataset$Sepal.Length[dataset$Species == "setosa"] sepalLength_versicolor <- dataset$Sepal.Length[dataset$Species == "versicolor"] petalLength_setosa <- dataset$Petal.Length[dataset$Species == "setosa"] petalLength_versicolor <- dataset$Petal.Length[dataset$Species == "versicolor"] n_setosa <- length(sepalLength_setosa) n_versicolor <- length(sepalLength_versicolor) r_setosa <- cor(sepalLength_setosa, petalLength_setosa) #Setosa種の標本相関係数 r_versicolor <- cor(sepalLength_versicolor, petalLength_versicolor) #Versicolor種の標本相関係数 z_setosa <- FisherZ(r_setosa) z_versicolor <- FisherZ(r_versicolor) stat <- (z_setosa - z_versicolor) / sqrt(1 / (n_setosa - 3) + 1 / (n_versicolor - 3)) pValue <- 2 * (1 - pnorm(abs(stat))) |
検定統計量statとp値pValueはそれぞれ次の通りです。
1 2 3 4 | > stat [1] -3.434362 > pValue [1] 0.0005939503 |
ポイント
「p値=0.0005939503 < 0.05」から、優位水準0.05の下でsetosaとversicolorによってSepal.LengthとPetal.Lengthの母相関係数は等しいとは言えないことが分かりました。
これらの検定結果をcsv等に保存するには、次のようにデータフレームに格納します。
関数write.csvを用いることで検定結果をcsvファイルに出力することが可能です。
1 2 | resultTable <- data.frame(stat = stat, p.value = pValue, row.names = NULL) write.csv(resultTable, "相関係数の差の検定（フィッシャーのZ変換）.csv", row.names = FALSE) |
相関係数の差の検定結果
まとめ
R言語でフィッシャーのZ変換を用いた相関係数の検定の実行例を紹介しました。
関数fisherZを用いることで相関係数の検定、相関係数の差の検定を行うことができます。
