テトラコリック相関の初期
書誌情報
テーブルの説明
Pearson(1900)が発表したテトラコリック相関では以下のようなクロス表を想定して相関係数を求める。
(p.437)
左の図の実線はx軸とy軸である。点線の表す平面が二変量正規分布を切って4つの区画a,b,c,dに分けられてると考える。周辺確率を用いてhとkの値を求めることができる。
さて,この表から相関係数が求めるには以下の式が使われる。
$ \frac{d}{N} = \frac{b+d}{N} \cdot \frac{c+d}{N}+S _ {1}^{\infty} (\frac{r^n}{n!}HK\bar{v} _ {n-1}\bar{w} _ {n-1}) \tag{1}
ここでHとKは正規曲線の縦座標(確率密度)である(Sは今の$\Sigma$)。$\bar{v}と\bar{w}$ は以下の式で与えられる。
$ \bar{v} = h\bar{v} _ {n-1} - (n-1)\bar{v} _ {n-2}, \bar{v} _ 0 = 1, \bar{v} _ 1 = h
$ \bar{w} = k\bar{w} _ {n-1} - (n-1)\bar{w} _ {n-2}, \bar{w}_0 = 1, \bar{w} _ 1 = k
これらを用いて,
$ \tau_n = \frac{H\bar{v} _ s{n-1}}{\sqrt{n!}}
$ \tau_n' = \frac{K\bar{w} _ {n-1}}{\sqrt{n!}}
とすると最初の(1)式は以下のように書き換えることができる。
$ \frac{d}{N} = \frac{b+d}{N}\cdot \frac{c+d}{N} + S_{1}^{\infty}(\tau_n \tau_n'r^n) \tag{ 2 }
この論文では,$\tau_6$までの数表が乗っていて,相関表からテトラコリック相関rを簡単に求められるようになっている。
テーブル利用の実際
論文では例として次のようなテーブルが挙げられている
こんな感じの2×2の表を使う,
table:table
1668 (a) 131 (b)
137 (c ) 64 (d)
この値を用いて,論文巻末の表を見る。巻末にはこんな感じで表が載っている。
(p.443)
表を用いると次のように値が求まる。
table:table
$ \tau_1$ $\tau _ 2$ $\tau _ 3$ $\tau _ 4$ $\tau _ 5$ $\tau _ 6$
0.17228 0.15787 0.04779 -0.06018 -0.06693 0.00854
0.17614 0.159276 0.04567 -0.06275 -0.06652 0.01111
これを(1)式に代入すると,
[$ 0.032000 = 0.009799+0.030345r+0.025142r^2+0.002183r^3+0.003776r^4+0.004452r^5+0.000095r^6 \tag{3}
となる。
あとはこの方程式を解けば$r=0.501$のようにテトラコリック相関係数が求まるそうである(論文だとニュートン法で解いたと書いてある)。
Rで実験
データをまず入力。
code:R
dat <- matrix(c(1668,137,131,64),2,2)
周辺頻度からhとkおよびHとKを計算する。
code:R
N <- sum(dat)
b_plus_d_byN <- colSums(dat)2/N # 0.0975
c_plus_d_byN <- rowSums(dat)2/N # 0.1005
h <- qnorm(b_plus_d_byN, lower.tail = F)
k <- qnorm(c_plus_d_byN, lower.tail = F)
H <- dnorm(h)
# 0.1722765
K <- dnorm(k)
# 0.1761384
τを求める計算式を作る。愚直に数式を入力する。
code:R
tau <- numeric(6)
v_bar <- numeric(6)
for(n in 1:6){
if(n==1){
taun <- (H*1)/sqrt(factorial(n)) } else if (n==2){
taun <- H * v_barn-1 /sqrt(factorial(n)) } else{
v_barn <- h * v_barn-1 - (n-1) * v_barn-2 taun <- (H * v_barn-1)/sqrt(factorial(n))} }
# tauの中身
1 0.172276537 0.157867340 0.047785511 -0.060181440 -0.066934070 0.008538122 同じ要領で,もう一つの方のtauも求めて(3)式を作る。Rでニュートン法を使う際にはpolyroot()の関数を用いる。この関数では,多項式の係数を並べたベクトルを引数として与えると,多項式=0となる解を計算してくれる。
code:R
# dat2,2/Nは(2)式の左辺を右辺に移行したもの coef_vec <- c(b_plus_d_byN * c_plus_d_byN - dat2,2/N, tau * tau2 ) round(polyroot(coef_vector),3)
# 結果のベクトル
1 0.502+0.00i -1.544+0.67i -1.544-0.67i 0.866-1.68i 0.866+1.68i -46.007+0.00i
1番目の引数に論文とほぼ同じ数値であるr=0.502があり,テトラコリック相関が求まったことが分かる。