テトラコリック相関の初期

テトラコリック相関係数関連の論文を読んだメモ。テトラコリック相関を計算するときに便利な表を作ったぞという話。

書誌情報

Everitt, P. F. (1910). Tables of the tetrachoric functions for fourfold correlation talbes. Biometrika, 7(4), 437–451. https://doi.org/10.1093/biomet/7.4.437

テーブルの説明

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の表を使う，

この値を用いて，論文巻末の表を見る。巻末にはこんな感じで表が載っている。 (p.443)

表を用いると次のように値が求まる。

これを(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で実験

データをまず入力。

dat <- matrix(c(1668,137,131,64),2,2)

周辺頻度からhとkおよびHとKを計算する。

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

τを求める計算式を作る。愚直に数式を入力する。

tau <- numeric(6)
v_bar <- numeric(6)

for(n in 1:6){
  if(n==1){
    v_bar[n] <- h
    tau[n] <- (H*1)/sqrt(factorial(n))
  } else if (n==2){
    v_bar[n] <- h^2 - 1
    tau[n] <- H * v_bar[n-1] /sqrt(factorial(n))
  } else{
    v_bar[n] <- h * v_bar[n-1] - (n-1) * v_bar[n-2]
    tau[n] <- (H * v_bar[n-1])/sqrt(factorial(n))  
  }
}
## tauの中身
[1]  0.172276537  0.157867340  0.047785511 -0.060181440 -0.066934070  0.008538122

同じ要領で，もう一つの方のtauも求めて(3)式を作る。Rでニュートン法を使う際にはpolyroot()の関数を用いる。この関数では，多項式の係数を並べたベクトルを引数として与えると，多項式=0となる解を計算してくれる。

# dat[2,2]/Nは(2)式の左辺を右辺に移行したもの
coef_vec <- c(b_plus_d_byN * c_plus_d_byN - dat[2,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$があり，テトラコリック相関が求まったことが分かる。

• ←Previous
• Next→