はじめに
最初に、本講義に必要なパッケージをロードする。初めてMASSパッケージを使う場合はインストールが必要となる。
require(MASS)## Loading required package: MASS相関
様々な場面において、2つの変数間の関係の強さを評価し、そこから何らかの推論を行いたい場合がある。本講義では、2変量間の関連の強さを表す統計量である相関係数を計算し、得られた相関係数をもとに推論を行うための統計的手法について説明する。
ピアソンの相関係数
\(X\)と\(Y\)が2つの正規分布\(N(0,\sigma_X^2)\)と\(N(0, \sigma_Y^2)\)に従っており、\(X\)と\(Y\)の間の相関係数が\(\rho_{XY}\)（以降、単に\(\rho\)と表す）であるとする。複数の標本\(1, 2, ..., n\)について計測された\(X\)と\(Y\)の観察値の組\((X_1, Y_1),(X_2, Y_2),...,(X_n,Y_n)\)が得られている場合に、その相関係数\(\rho\)について、推定・検定を行う方法について説明する。
観察値の組\((X_1, Y_1),(X_2, Y_2),...,(X_n,Y_n)\)から、以下の平方和を計算することができる。 \[ S_{xx} = \sum_{i = 1}^n (X_i - \bar X)^2 \\ S_{yy} = \sum_{i = 1}^n (Y_i - \bar Y)^2 \\ S_{xy} = \sum_{i = 1}^n (X_i - \bar X)(Y_i - \bar Y) \] このとき、ピアソンの相関係数\(\rho\)の推定値は、 \[ r = \frac {S_{xy}}{\sqrt{S_{xx}S_{yy}}} \] として表せる。
この推定値の別の表現の仕方は、 \[ r = \frac {\sum_i X_iY_i - n\bar X \bar Y}{\sqrt{\left(\sum_i X_i^2 - n \bar X^2\right)\left(\sum_i Y_i^2 - n \bar Y^2\right)}} \] である。これは、例えば、 \[ S_{xy} = \sum_{i = 1}^n (X_i - \bar X)(Y_i - \bar Y) \\ = \sum_{i = 1}^n (X_i Y_i - X_i \bar Y - \bar X Y_i + \bar X \bar Y) \\ = \sum_{i = 1}^n X_i Y_i - \bar X \sum_{i = 1}^n Y_i - \bar Y \sum_{i = 1}^n X_i + n \bar X \bar Y \\ = \sum_{i = 1}^n X_i Y_i - n \bar X \bar Y \] というかたちで簡単に導くことができる。
様々な相関をもつ二変量正規分布からの乱数を発生して、図を描いてみる。多変量正規分布の乱数を発生するための関数mvrnormは、パッケージMASSに含まれている。
## 11.2
op <- par(mfrow = c(2,2))
r <- 0.1
xy <- mvrnorm(100, c(0,0), Sigma = matrix(c(1,r,r,1), nrow = 2))
plot(xy, main = r, col = "green")
r <- 0.6
xy <- mvrnorm(100, c(0,0), Sigma = matrix(c(1,r,r,1), nrow = 2))
plot(xy, main = r, col = "green")
r <- -0.6
xy <- mvrnorm(100, c(0,0), Sigma = matrix(c(1,r,r,1), nrow = 2))
plot(xy, main = r, col = "green")
r <- 0.9
xy <- mvrnorm(100, c(0,0), Sigma = matrix(c(1,r,r,1), nrow = 2))
plot(xy, main = r, col = "green")
par(op)なお、2つの独立した正規乱数\(X\)と\(Y\)から、\(X\)と相関\(\rho\)をもつ正規乱数\(Z\)を \[ Z = \rho X+\sqrt{1-\rho^2}Y \] として生成することもできる。
実際に生成して、図を描いてみる。
r <- 0.9
x <- rnorm(100)
y <- rnorm(100)
z <- r * x + sqrt(1 - r^2) * y
plot(x, z, main = r, col = "green")
偏相関係数
2変数\(X\),\(Y\)の相関関係が、それら2変数とは別の変数\(Z\)との関係から生じている場合がある。このような場合は、\(Z\)によって生じている相関を除いた場合の相関を計算することができる。上述した乱数発生の方法をもちいて、正規乱数\(Z\)と、\(Z\)と相関\(\rho\)をもつ正規乱数\(X\)と\(Y\)を独立に生成してみる。
ここでは、\(X\)と\(Z\)が相関0.8、\(Y\)と\(Z\)が相関0.8をもつとして、乱数\(X\), \(Y\), \(Z\)を生成する。
r <- 0.8
z <- rnorm(100)
x <- r * z + sqrt(1 - r^2) * rnorm(100)
y <- r * z + sqrt(1 - r^2) * rnorm(100) # Xと同じ計算法で独立して生成する
df <- data.frame(z, x, y)
pairs(df, col = "green")
すると、\(X\)と\(Y\)間にも相関関係が生じる。
このとき、次のようにすると、\(Z\)をとの関係を考慮に入れた（\(Z\)の効果を除いた）\(X\)と\(Y\)の相関係数（偏相関係数 partial correlation）を計算できる。
今、\(r_{xy}, r_{xz}, r_{yz}\)を変数\(X, Y, Z\)間の相関とする。このとき、\(Z\)を考慮に入れた偏相関係数は、 \[ r_{xy.z} = \frac {r_{xy} - r_{xz}r_{yz}} {\sqrt{1-r_{xz}^2}\sqrt{1-r_{yz}^2}} \] として計算される。
上で生成した\(X, Y, Z\)のデータで計算をしてみる。
(r_xy <- cor(x, y))## [1] 0.6000321(r_xz <- cor(x, z))## [1] 0.7381909(r_yz <- cor(y, z))## [1] 0.7445033(r_xy.z <- (r_xy - r_xz * r_yz) / (sqrt(1 - r_xz^2) * sqrt(1 - r_yz^2)))## [1] 0.1120114\(Z\)による影響を考慮した場合には、\(X\)の\(Y\)の相関はほとんど0となることが分かる。すなわち、\(X\)と\(Y\)間の関係は、もっぱら\(Z\)との関係によって生じていたと解釈できる。
なお、\(X\)による影響を考慮して、\(Y\)と\(Z\)の相関を計算した場合は、
(r_yz.x <- (r_yz - r_xy * r_xz) / (sqrt(1 - r_xy^2) * sqrt(1 - r_xz^2)))## [1] 0.5588085相関係数の値は小さくなるが、\(X\)と\(Y\)間の偏相関係数のように0にはならない。
相関係数の検定
帰無仮説\(H_0: \rho = 0\)を、どれか一つの対立仮説\(\rho >, \ne, < 0\)に対して検定するために、\(t\)統計量 \[ t = r \sqrt{\frac {n-2}{1-r^2}} \] が用いられる。この統計量は、帰無仮説のもとで自由度\(df = n-2\)の\(t\)分布に従う。
ここで、各対立仮説とそれに対応する\(p\)値の計算方法は、以下の通り。
| 対立仮説 | \(\alpha\)水準における棄却域 | \(p\)値（Rでの計算法） |
|---|---|---|
| \(H_1: \rho > 0\) | \([t_{df, 1-\alpha}, \infty)\) | 1 - pt(t, df) |
| \(H_1: \rho \ne 0\) | \((-\infty, t_{df, \alpha/2}] \cup [t_{df, 1-\alpha/2}, \infty)\) | 2 * pt(-abs(t), df) |
| \(H_1: \rho < 0\) | \((-\infty, t_{df, \alpha}]\) | pt(t, df) |
例1 ラットは迷路を抜けだすための方法を学習できるのか
以下のデータは、ラットに迷路を体験させた回数\(X\)と、最後の回に迷路を抜けるまでにかかった時間\(Y\)の関係を表したものである。
- 相関係数を計算せよ。
- 迷路の学習は有意か？ 母相関係数\(\rho\)が0であるという帰無仮説について、\(\rho\)は負であるという対立仮説に対して検定せよ。
まずは、\(X\)と\(Y\)を図示して確認してみる。
X <- c(8, 9, 6, 5, 3, 6, 3, 2)
Y <- c(10.9, 8.6, 11.4, 13.6, 10.3, 11.7, 10.7, 14.8)
plot(X, Y)
体験回数\(X\)が増えるに従って、迷路を抜け出すまでの時間\(Y\)は減少しているようにみえる。
次に、相関係数を計算する。まず、1つめの計算法に従って計算する。
sxx <- sum((X - mean(X))^2)
syy <- sum((Y - mean(Y))^2)
sxy <- sum((X - mean(X)) * (Y - mean(Y)))
(r <- sxy / sqrt(sxx * syy))## [1] -0.5687298もちろん、2つめの計算法でも同じ値が得られる。
sxx <- sum(X^2)
syy <- sum(Y^2)
sxy <- sum(X*Y)
mx <- mean(X)
my <- mean(Y)
n <- length(X)
(r <- (sxy - n * mx * my) / sqrt((sxx - n * mx^2) * (syy - n * my^2)))## [1] -0.5687298次に、検定を行う。まずは、\(t\)値を計算し、\(H_1: \rho < 0\)に対して\(H_0\)を検定する。
t <- r * sqrt((n - 2) / (1 - r^2))
t## [1] -1.693685pt(t, n - 2)## [1] 0.07063285%水準で有意ではない。すなわち、このデータからは、ラットが迷路を体験して抜け出すための方法を学習しているとは言えない。
Rの関数corを使っても計算できる。
cor(X,Y)## [1] -0.5687298cor.test(X,Y, alternative = "less") # 下側確率で計算##
## Pearson's product-moment correlation
##
## data: X and Y
## t = -1.6937, df = 6, p-value = 0.07063
## alternative hypothesis: true correlation is less than 0
## 95 percent confidence interval:
## -1.00000000 0.08971573
## sample estimates:
## cor
## -0.5687298なお、偏相関の検定には、検定統計量 \[ t = r \sqrt{\frac {n - q - 2}{1 - r^2}} \] を用いる。自由度は\(df = n - q - 2\)である。ここで、\(q\)は、偏相関の計算で考慮した他の変数の数である。
例2 母親の教育水準と子供の死亡率の関係
標本数\(n=28\)のデータにおいて、母親の教育水準\(X\)と子供の死亡率\(Y\)に負の相関\(r_{xy} = -0.6\)がみられた。ただし、\(X\)と\(Y\)の関係は、その家族の社会経済的状況\(Z\)を考慮に入れて剣とうされるべきだと考えられる。なお、\(Z\)は\(X\)、\(Y\)と、\(r_{xz} = 0.65\)および\(r_{yz} = -0.7\)の相関があった。
- \(r_{xy.z}\)を求めなさい
- 母集団における\(\rho_{xy.z}\)の有意性を検定しなさい。
まずは、\(Z\)を考慮にいれて、\(X\)と\(Y\)間の偏相関係数を計算する。
rxy <- -0.6
rxz <- 0.65
ryz <- -0.7
# 偏相関の計算
rxy_z <- (rxy - rxz * ryz) / sqrt((1 - rxz^2)*(1 - ryz^2))
rxy_z## [1] -0.2671818\(Z\)を考慮に入れた\(X\), \(Y\)間の偏相関係数（\(r = -0.27\)）は、\(X\), \(Y\)間の相関係数（\(r = -0.60\)）に比べて0に近い値をとる。
次に、検定を行う。まずは、相関係数の検定を行ってみる。
# 検定
n <- 28
# 相関係数の検定
t <- rxy * sqrt((n - 2) / (1 - rxy^2))
t## [1] -3.8242652 * pt(- abs(t), n - 2)## [1] 0.00073809610.1%水準で有意である。
つぎに、偏相関係数の検定を行う。
# 偏相関係数の検定
t <- rxy_z * sqrt((n - 1 - 2)/ (1 - rxy_z^2))
t## [1] -1.3863072 * pt(- abs(t), n - 1 - 2)## [1] 0.17789395%水準で有意ではない。
社会経済的状況\(Z\)を考慮しない場合は、母親の教育水準\(X\)と子供の死亡率\(Y\)の相関係数は高度に有意になる。しかし、偏相関係数は、通常の相関係数に比べて0に近い値\(r_{xy.z} = -0.27\)になり、両側検定の結果も有意ではない。この結果から、母親の教育水準\(X\)と子供の死亡率\(Y\)の関係は、相関係数でみると関連があるように見えるが、社会経済的状況\(Z\)を考慮に入れると両者の関係は有意ではない、ということが分かる。
相関係数の信頼区間
母相関係数\(\rho\)の信頼区間は、直接計測される値から求められるのではなく、変換された値に対して得られる。この変換は、Fisherの\(z\)変換として知られているものである。
\((X_{11},X_{21}),...,(X_{1n},X_{2n})\)を二変量正規分布\(N(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2,\rho)\)からの標本とするまた、\(\bar{X_i} = \sum_{j=1}^n X_{ij} / n \;(i = 1,2)\)とします。このとき、ピアソンの相関係数 \[ r = \frac {\sum_{i=1}^n (X_{1i} - \bar X_1)(X_{2i} - \bar X_2)}{\biggl[\sum_{i=1}^n(X_{1i}-\bar X_1)^2 \sum_{i=2}^n(X_{2i}-\bar X_2)^2\biggr]^{1/2}} \]
は複雑な分布をもつ。しかし、\(r\)の漸近分布は正規分布\(N(\rho, (1-\rho^2)^2 / n)\)であることが知られている。ただし、\(r\)の漸近分布は、分散が平均\(\rho\)の関数\(\sigma^2(\rho) = \frac {(1-\rho^2)^2}{n}\)であることから、以下の変換を行って分散の平均依存性を除く。
\(r\)と\(\rho\)を \[ w = \frac 1 2 \log \left(\frac{1 + r}{1-r}\right) = \text{atanh} \;r \\ \zeta = \frac 1 2 \log \left(\frac{1 + \rho}{1-\rho}\right) = \text{atanh} \;\rho \] と変換すると、\(w\)の分布は、正規分布\(N\left(\zeta, \frac 1 {n-3}\right)\)で近似できる。これが、相関係数のFisherの\(z\)変換である。これは、逆双曲線正接変換とよばれる変換であり、Rでは、atanh関数を用いて変換できる。
逆変換は、 \[ r = \tanh w = \frac {\exp{(2w)}-1}{\exp{(2w)}+ 1} \] となる。これは、双曲線正接変換である。Rでは、tanh関数を用いて計算できる。
上述した変換を例示するために、\(n = 30\)の値の組について相関\(r = 0.6\)をもつ正規分布からの乱数を生成する。これを 10,000回繰り返し、得られた10,000標本のヒストグラムを描いみる。また、\(z\)変換をした\(r\)のヒストグラムも描いてみる。後者には、その正規近似\[ N\left(\text{atanh}\: 0.6, \frac{1}{30 - 3}\right) \]を重ね書きしてみる。
まずは、10000回のシミュレーションをおこない、\(r\)とその逆双曲線正接\(w\)を計算する。
n <- 10000 # シミュレーションの回数
rho <- 0.6 # 真の相関係数
r <- rep(NA, n) # 入れ物の準備
for(i in 1:n) { # for文で、シミュレーションをn回回す。
x <- rnorm(30) # 30個の乱数を発生
y <- rho * x + sqrt(1 - rho^2) * rnorm(30) # この計算によりxとyは相関rho（真値）をもつ乱数となる
r[i] <- cor(x, y)
}
w <- atanh(r) # 逆双曲線正接変換まずは、\(r\)のヒストグラムの描画を行う。
hist(r, freq = F)
右側に偏った分布となっていて、正規分布とみなすのは難しそうである。
次に、\(w\)のヒストグラムを描画し、近似される正規分布\(N\left(\zeta, \frac 1 {n-3}\right)\)を描く。
hist(w, freq = F) # wのヒストグラムを書く
x <- seq(0, 2, 0.01) # 理論的に期待される正規分布を描く準備（xを準備）
y <- dnorm(x, atanh(rho), sqrt(1/(30 - 3))) # xに対応するyを計算。
# ばらつきは標準偏差で指定することに注意
lines(x, y, col = "green") # 正規分布の曲線を描く
このように、\(w\)の標本分布は、平均\(\zeta = \text{atanh}\;\rho\)および分散\(1/(n-3)\)の正規分布で精度良く近似される。この近似は、標本数が20を超えている場合には十分な精度を与えることが知られている。 \(\zeta\)の\((1-\alpha)100\)%信頼区間は、以下のように計算できる。 \[ [w_L, w_U] = [w - \frac {z_{1-\alpha/2}} {\sqrt{n-3}}, w + \frac {z_{1-\alpha/2}} {\sqrt{n-3}}] \] ここで、\(w= \text{atanh}\;r\)である。
逆変換を行うと、\(\rho\)についての信頼区間を求めることができる。すなわち、 \[ [r_L, r_U] = [\tanh w_L, \tanh w_U] = \left[\frac {\exp{(2w_L)}-1}{\exp{(2w_L)}+ 1}, \frac {\exp{(2w_U)}-1}{\exp{(2w_U)}+ 1}\right] \]
例3 Fisherの\(z\)変換
\(r = -0.5687\)、\(n = 8\)のときに\(\zeta = \frac 1 2 \log \frac {1 + \rho}{1 - \rho}\)の95%信頼区間を求めよ（これは、例1のデータに基づいている）。また、\(\rho\)の信頼区間も逆変換を用いて計算しなさい。
まずは、Fisherの\(z\)変換を行う。
r <- -0.5687
n <- 8
(w <- atanh(r)) # Fisherのz変換## [1] -0.6455993次に、\(w\)の信頼区間を求める。
alpha <- 0.05
wL <- w - qnorm(1 - alpha/2) / sqrt(n - 3)
wU <- w + qnorm(1 - alpha/2) / sqrt(n - 3)
# zetaの信頼区間
c(wL, wU)## [1] -1.5221219 0.2309232\(\rho\)の信頼区間を逆変換を用いて求める
# rhoの信頼区間
tanh(c(wL, wU)) ## [1] -0.9090667 0.2269042cor.testを使って信頼区間を求める。先ほどの例1では、下側の信頼区間が計算されていた。ここでは両側検定として信頼区間を求める。
X <- c(8, 9, 6, 5, 3, 6, 3, 2)
Y <- c(10.9, 8.6, 11.4, 13.6, 10.3, 11.7, 10.7, 14.8)
cor.test(X, Y)##
## Pearson's product-moment correlation
##
## data: X and Y
## t = -1.6937, df = 6, p-value = 0.1413
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9090744 0.2268625
## sample estimates:
## cor
## -0.5687298同じ信頼区間が計算される。
\(\rho = \rho_0\)の検定
帰無仮説\(H_0\)が\(\rho_0 = 0\)のときは、先に述べたように統計量\(t = r \sqrt {\frac {n-2}{1-r^2}}\)を用いた\(t\)検定により検定でる。しかし、この検定を、\(\rho_0 \ne 0\)の場合にも拡張することはできない。そこで、\(\rho_0 \ne 0\)の場合は、正規近似に基づく検定を行うことになる。 \(H_0: \rho = \rho_0\)のもとで、検定統計量 \[ z = \sqrt{n - 3} \; (\text{atanh}\;r-\text{atanh}\;\rho_0)=\frac {\sqrt{n - 3}}{2} \left[\log \left(\frac{1 + r}{1 - r}\right) - \log \left(\frac{1 + \rho_0}{1 - \rho_0}\right)\right] \] は、近似的に標準正規分布に従う。これにより、以下のルールのもとで検定を行うことができる。
| 対立仮説 | \(\alpha\)水準における棄却域 | \(p\)値（Rでの計算法） |
|---|---|---|
| \(H_1: \rho > \rho_0\) | \([z_{1-\alpha}, \infty)\) | 1 - pnorm(z) |
| \(H_1: \rho \ne \rho_0\) | \((-\infty, z_{\alpha/2}] \cup [z_{1-\alpha/2}, \infty)\) | 2 * pnorm(-abs(z)) |
| \(H_1: \rho < \rho_0\) | \((-\infty, z_{\alpha}]\) | pnorm(z) |
例4 \(H_0: \rho=0.5\)に対する検定
2つの変数間の母相関係数が\(\rho=0.5\)となることが期待される実験を行い、\(r=0.23\)という推定値が30組の標本が得られとする。期待される相関係数に比べて小さな値が得られたが、これが偶然によるものかどうかを検討したい。そこで、\(\rho=0.5\)であるという仮説を\(\rho \ne0.5\)に対して検定する。
まずは、正規分布にしたがう検定統計量\(z\)を計算する。
r <- 0.23
n <- 30
rho0 <- 0.5
(z <- sqrt(n - 3) * (atanh(r) - atanh(rho0)))## [1] -1.637394つぎに、検定をおこなう。
2 * pnorm(-abs(z))## [1] 0.1015481検定の結果、5％水準で有意ではない。すなわち、30組の標本では、\(r=0.23\)という推定値が\(\rho=0.5\)のもとで10回に1回程度は偶然に得られる値であることが分かる。
2つの相関係数の同一性の検定
2つの相関係数\(\rho_1\)、\(\rho_2\)が等しいかどうかに興味がある場合がある。今、1つめの母集団から、\(X_i^{(1)},Y_i^{(1)}\;(i=1,...,n_1)\)が標本として得られ、それとは独立して、2つめの母集団から\(X_i^{(2)},Y_i^{(2)}\:(i=1,...,n_2)\)が標本として得られたとする。ここで、それぞれの標本から相関係数の推定値\(r_1, r_2\)が計算されたとする。このとき、これら2つの独立した母集団間で、母相間係数\(\rho_1\)、\(\rho_2\)が等しいかどうか検定する。
\(H_0: \rho_1 = \rho_2\)を\(H_1: \rho_1 >,\ne,< \rho_2\)に対して検定をおこなう検定統計量は、標本相関\(r_1, r_2\)のFisherの\(z\)変換である\(w_1, w_2\)をもとに表すことができる。 \[ z = \frac{w_1 - w_2}{\sqrt{\frac 1 {n_1-3} + \frac 1 {n_2 - 3}}} \] ここで、
\[ w_i = \text{atanh}\;r_i = \frac 1 2 \log \frac {1 + r_i}{1-r_i}\;\; (i = 1,2) \] で、\(n_1,n_2\)は2つの集団における標本（\(X_i, Y_i\)のペア）数である。
| 対立仮説 | \(\alpha\)水準における棄却域 | \(p\)値 （Rでの計算法） |
|---|---|---|
| \(H_1: \rho_1 > \rho_2\) | \([z_{1-\alpha}, \infty)\) | 1 - pnorm(z) |
| \(H_1: \rho_1 \ne \rho_2\) | \((-\infty, z_{\alpha/2}] \cup [z_{1-\alpha/2}, \infty)\) | 2 * pnorm(-abs(z)) |
| \(H_1: \rho_1 < \rho_2\) | \((-\infty, z_{\alpha}]\) | pnorm(z) |
例5 2つの計測法の比較
手法2に対して、新たな計測法である手法1を開発した。手法1を用いて\(n_1=60\)回の計測を行った結果、真値と計測値の相関係数の推定値が\(r_1 = 0.862\)であった。それとは独立して、手法2を用いて\(n_2=49\)回の計測を行った結果、真値と計測値の相関係数の推定値が\(r_2 = 0.720\)であった。この結果から、新たな手法1は、旧手法2に比べて精度が高いといえるか。有意水準\(\alpha = 0.05\)で検定せよ。
まずは、\(r_1, r_2\)のFisherの\(z\)変換\(w_1, w_2\)を計算する。
r1 <- 0.862
r2 <- 0.720
# Fisherのz変換
(w1 <- atanh(r1))## [1] 1.301076(w2 <- atanh(r2))## [1] 0.907645次に、\(\rho_1\)と\(\rho_2\)が等しいかどうか検定を行う。ここでは、新手法1が旧手法2に対して精度が高いことに興味があるので、上側の片側検定を行う。
n1 <- 60
n2 <- 49
# 2つの母相関係数 rho_1とrho_2が等しいかどうかを検定
z <- (w1 - w2) / sqrt(1 / (n1 - 3) + 1 / (n2 - 3))
# 上側の片側検定を行う
1 - pnorm(z)## [1] 0.023570655%水準で有意である。したがって、新手法1は旧手法2より精度が高いといえる。
複数の相関係数の同一性の検定
ここでは、3つ以上の母集団から得られた相関係数をもとに、仮説\(H_0: \rho_1 = \rho_2 = \cdots = \rho_k\)を検定することを考える。
\(r_i\)を\(n_i\)個の標本（ペア）から計算された相関係数として、\(w_i = \text{arctanh}\;r_i = \frac 1 2 \log \frac{1 + r_i}{1 - r_i}\)をFisherの\(z\)変換とする。ここで、 \[ N = (n_1 - 3) + (n_2 - 3) + \cdots + (n_k - 3) = \sum_{i=1}^k n_i - 3k \] \[ \bar w = \frac {(n_1-3)w_1 + (n_2-3)w_2 + \cdots + (n_k - 3)w_k}{N} \] とすると、統計量 \[ \chi^2 = (n_1-3)w_1^2+(n_2-3)w_2^2+\cdots+(n_k-3)w_k^2-N\bar w^2 \] は、帰無仮説\(H_0: \rho_1 = \rho_2 = \cdots = \rho_k\)のもとで自由度\(k-1\)のカイ二乗分布に従う。
例6 Fisherのアヤメデータ
Fisherのアヤメデータでは、アヤメ3種（setosa, versicolor, virginica）の萼片の長さ(Sepal.Length)と花弁の長さ(Petal.Length)間に以下の相関がある。これら3種の相関に違いがあるかどうか検定せよ。
| 品種 | setosa | versicolor | virginica |
|---|---|---|---|
| 標本数 | 50 | 50 | 50 |
| \(r\) | 0.267 | 0.754 | 0.864 |
# 標本数等のセット
n <- c(50, 50, 50)
k <- 3 # k <- length(n)としてもよい
N <- sum(n) - 3 * k
# 相関のセット
r <- c(0.267, 0.754, 0.864)
w <- atanh(r)
wbar <- sum((n-3) * w) / N # sum((n-3) * w)を(n-3)%*%wとしてもよい
(chi2 <- sum((n-3) * w^2) - N * wbar^2) # sum((n-3) * w^2)を(n-3)%*%w^2としてもよい## [1] 26.32925(pval <- 1 - pchisq(chi2, k - 1))## [1] 1.917238e-06計算される統計量の値は26.3となり、自由度2（\(k=3\)）の\(\chi^2\)分布として検定すると3種間の相関係数の違いは、高度に有意である。
スピアマン（Spearman）の順位相関
\(n\)個の標本において、\(R_1, ...,R_n\)を\(X\)の昇順の順位、\(S_1,...,S_n\)を\(Y\)の昇順の順位とする。例えば、\((R_1,S_1)=(3,5)\)とすると、1番目の標本の\(X_1\)と\(Y_1\)が、それぞれ、\(X\)、\(Y\)の中で3番目、5番目に小さい値となる。このとき、スピアマンの相関係数は以下のように定義される。 \[ \hat \rho = \frac{\sum_{i=1}^n (R_i - \bar R)(S_i - \bar S)}{\sqrt{\sum_{i=1}^n(R_i-\bar R)^2\cdot \sum_{i=1}^n(S_i-\bar S)^2}} \]
この式は、より簡単な式に変形できる。それには、 \[ \bar R = \bar S = \frac{n+1} 2 \\ \sum(R_i - \bar R)^2=\sum(S_i - \bar S)^2 = nVar(R_i) = \frac{n(n^2-1)} {12} \] という性質を利用する。
\(D\)を順位の違いを、\(D_i = R_i - S_i\)と定義すると、\(\bar R = \bar S\)より、以下のようになる。 \[ D_i = (R_i-\bar R)-(S_i-\bar S) \] このとき、 \[ \sum_{i=1}^n D_i^2 = \sum_{i=1}^n(R_i - \bar R)^2 + \sum_{i=1}^n(S_i - \bar S)^2 - 2\sum_{i=1}^n(R_i - \bar R)(S_i - \bar S) \] となるため、上式の最後の項は、以下のように表せる。 \[ \sum_{i=1}^n (R_i - \bar R)(S_i - \bar S) = \frac{n(n^2-1)}{12}-\frac1 2 \sum_{i=1}^n D_i^2 \]
さらに、上式の両側を\(\sqrt{\sum_{i=1}^n(R_i - \bar R)^2\cdot \sum_{i=1}^n (S_i -\bar S)^2}=\sum_{i=1}^n(R_i-\bar R)^2=n(n^2-1)/12\)で割ると、最終的に、 \[ \hat \rho = 1 - \frac{6\sum_{i=1}^n D_i^2}{n(n^2-1)} \] が得られる。
ピアソンの相関係数と同じように、スピアマンの順位相関も\(-1\)から\(1\)の値をとる。順位が完全に一致している場合は、順位の違いは0となり、\(\hat \rho =1\)となります。順位が完全に逆転している場合は、\(R_i = n - S_i + 1\)となり、\(\sum D_i^2\)が最大化される。
標本が十分に大きい場合は、スピアマンの統計量は、正規分布で近似できる。なお、\(n>10\)の場合、 \[ Z = (\hat \rho - \rho) \sqrt{n-1} \sim N(0,1) \] が良い近似となることが知られています。
ケンドール(Kendall)のタウ(\(\tau\))
ケンドール（Kendall）は、19世紀に提案された依存性（dependence）について別の尺度を公式化した。そこでは、二変量のペアのうちどれだけのものが合致（concordant）しているか、すなわち、\(X\)と\(Y\)の間で符号が一致するか、を見つけ出す。一方の符号が正で、もう一方が負であれば、合致しない（discordant）とする。 \(X_i, Y_i\:(i=1,\dots,n)\)からは、\(\left(\begin{matrix}n\\2\end{matrix}\right)\)の異なるペアを選ぶことができる。ペア、\((X_i,Y_i),(X_j,Y_j)\)において、もし\(X_i\le X_j\)かつ\(Y_i \le Y_j\)、あるいは、\(X_i\ge X_j\)かつ\(Y_i \ge Y_j\)であれば、合致する（concordant）とする。逆に、\(X_i\le X_j\)かつ\(Y_i \ge Y_j\)あるいは、\(X_i\ge X_j\)かつ\(Y_i \le Y_j\)であれば、合致しない（discordant）とする。例えば、\((2,4)\)と\((1,-1)\)のペアは合致(\(2 \ge 1\)かつ\(4 \ge -1\)) し、\((-2,4)\)と\((1,-1)\)は合致しない(\(-2 \le 1\)だが\(4 \ge -1\))。
ケンドールの\(\tau\)統計量は、以下のように定義される。 \[ \hat \tau = \frac {2 S_r}{n(n-1)}\\ S_r = \sum_{i=1}^{n-1} \sum_{j=i+1}^n \text{sign}(r_j-r_i) \] ここで、\(r_i\)は、1つめの変数における順位にしたがって {1,2,…,n}の順に並べたときの、2つめの変数における順位である。すなわち、 \[ \left( \begin{matrix} 1 & 2 & ... & n \\ r_1 & r_2 & ... & r_n \end{matrix}\right) \] なお、合致の組数を\(n_C\)、非合致の組数を\(n_D=n-n_C\)とすると、 \[ \hat \tau = \frac {2(n_C-n_D)}{n(n-1)} \] と表せる。
例7
コムギの収量（\(X\)）とタンパク質含量（\(Y\)）のデータ（Snedecor and Cochran 1989より改変）について、スピアマンの順位相関、ケンドールの相関係数、ピアソンの相関係数を計算せよ。
X <- c(5, 8, 10, 11, 14, 16, 17, 18, 19, 20)
Y <- c(16.2, 14.2, 14.6, 18.3, 13.2, 13.0, 12.9, 13.4, 10.6, 12.8)
plot(X, Y)
相関の計算、検定を行う。
# speaman
cor(X, Y, method = "spearman")## [1] -0.830303# kendall
cor(X, Y, method = "kendall")## [1] -0.6444444# pearson
cor(X, Y)## [1] -0.7113396# testing H_0: rho = 0
cor.test(X, Y, method = "spearman")##
## Spearman's rank correlation rho
##
## data: X and Y
## S = 302, p-value = 0.005557
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.830303cor.test(X, Y, method = "kendall")##
## Kendall's rank correlation tau
##
## data: X and Y
## T = 8, p-value = 0.009148
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## -0.6444444cor.test(X, Y)##
## Pearson's product-moment correlation
##
## data: X and Y
## t = -2.8626, df = 8, p-value = 0.02107
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9261594 -0.1479987
## sample estimates:
## cor
## -0.7113396参考文献
本テキストを作成するにあたり参考にした書籍
Vidakovic B (2011) Statistics for Bioengineering Sciences With MATLAB and WinBUGS Support. Springer.
Snedecor GW and Cochran WG (1989) Statistical Methods. Eighth Edition. Iowa State University Press.