ポリコリック相関係数のお勉強

ポリコリック相関係数とは, 2つの順序付きのカテゴリー変数の間に2変量正規分布を仮定して母相関を最尤推定する方法である。Olsson(1979)の方法が有名である。Rのパッケージでも計算できるがお勉強ついでに自作してみる。解説はこの論文を参考にした。 https://www.ncbi.nlm.nih.gov/pubmed/15598100

データの準備

とりあえず順序つきのカテゴリー変数を作ってみる。2変量正規分布で乱数を生成してそれを離散化する。カテゴリー数は変数1は3, 変数2は4あたりにしてみる。

library(mvtnorm)
set.seed(1)
sigma <- matrix(c(1,0.5,0.5,1),2)
dat <- rmvnorm(100, mean=c(0,0), sigma=sigma)

cat1 <- cut(dat[,1], breaks=c(-Inf,-1,1,Inf), label=F)
cat2 <- cut(dat[,2], breaks=c(-Inf,-1,0,1,Inf), label=F)
cont.table <- table(cat1,cat2)

こんな感じの分割表が手に入る。

f:id:nekomosyakushimo:20190504011544p:plain

cor(dat) #0.5053
cor(cat1,cat2) #0.4049

カテゴリー変数になることにより相関係数が低くなっていることが分かる。

ポリコリック相関のやることはこの分割表から離散化前の相関係数を推定することである。

ステップ1 閾値の推定

Olsson(1979)の方法ではポリコリック相関係数を2つのステップで推定する。最初のステップは閾値の推定である。周辺確率を使って次のように定める。

$\displaystyle a_i = \Phi^{-1}_1(P_{i \cdot})$

ここで $a_i, 0...s$ は閾値である。同様に,

$\displaystyle b_j = \Phi^{-1}_1(P_{\cdot j})$

$b_j, 0...r$ はもう片方の変数の閾値である。

rではこんな感じでやると閾値が出る。後に使うので前後に-InfとInfを加えている。

# 閾値推定
marg.freq1 <- prop.table(rowSums(cont.table))
marg.freq2 <- prop.table(colSums(cont.table))

tau1 <- c(-Inf,qnorm(cumsum(marg.freq1)[-length(marg.freq1)]),Inf)
tau2 <- c(-Inf,qnorm(cumsum(marg.freq2)[-length(marg.freq2)]),Inf)

f:id:nekomosyakushimo:20190504012039p:plain

こんな感じで閾値が手に入る。

ステップ2 母相関の推定

さて, いま手に入れた閾値と分割表を使って母相関係数を最尤推定する。尤度は次のように求める。

$\displaystyle l = \ln K + \sum_1^{S} \sum_1^{J} n_{ij} \ln \pi_{ij}$

ここで, $K$ は定数, $n_{ij}$ は分割表のセルの度数であり, $pi_{ij}$ は変数がcell(i,j)に入る確率である。この確率というのは２変量正規分布を仮定すると,

$\displaystyle \pi_{ij} = \Phi_2(a_i, b_j) - \Phi_2(a_{i-1}, b_j) - \Phi_2(a_{i}, b_{j-1}) + \Phi_2(a_{i-1}, b_{j-1})$

で計算できる。ここで $\Phi_2$ は $\rho$ をパラメータに持つ２変量累積密度関数である。

Rでこれを計算する。mvtnormというパッケージのpmvnormはセルの範囲を指定すれば確率を計算してくれる。さっき作った閾値で範囲を指定しながら全セル分計算するようにしている。ここでぼは例えば母相関0.3の場合を計算してみる。

  pi <- numeric(length(cont.table))
  pi <- matrix(pi, length(marg.freq1), length(marg.freq2))
  for(ite.cat1 in 1:length(marg.freq1)){
    for(ite.cat2 in 1:length(marg.freq2)){
      pi[ite.cat1, ite.cat2] <- pmvnorm(mean=c(0,0),sigma=matrix(c(1,0.3,0.3,1),2),
                                    lower=c(tau1[ite.cat1], tau2[ite.cat2]),
                                    upper=c(tau1[ite.cat1+1],tau2[ite.cat2+1]))[1]
      }
  }

そうするとこんな感じで各セルの確率が出てくる。

f:id:nekomosyakushimo:20190504011617p:plain

で次は, 式にしたがってその確率の対数をとって度数をかけて合計する。

logL <- sum(cont.table * log(pi))
logL
#[1] -201.9102

これで母相関が0.3だとした場合の尤度が出てきた訳だ。定数はパラメータ変えようが関係ないから無視して大丈夫。あとは, 母相関を変えながら尤度が最も高くなる母相関を探せば良いので, 今までのものを次のようにループに組み込む。-0.99から0.99までを0.01刻みで動かしながら尤度を計算していく。

rho <- seq(-0.99,0.99, by=0.01)
logL <- numeric(length(rho))

for(ite.rho in 1:length(rho)){

  pi <- numeric(length(cont.table))
  pi <- matrix(pi, length(marg.freq1), length(marg.freq2))
  for(ite.cat1 in 1:length(marg.freq1)){
    for(ite.cat2 in 1:length(marg.freq2)){
      pi[ite.cat1, ite.cat2] <- pmvnorm(mean=c(0,0),sigma=matrix(c(1,rho[ite.rho],rho[ite.rho],1),2),
                                    lower=c(tau1[ite.cat1], tau2[ite.cat2]),
                                    upper=c(tau1[ite.cat1+1],tau2[ite.cat2+1]))[1]
      }
  }

  logL[ite.rho]<- sum(cont.table * log(pi))
}

#もっとも尤度が高い母相関は
rho[which.max(logL)]
#[1] 0.52

計算の結果, 母相関は0.52が最も尤もらしい値だということになった。

試しに、polycorパッケージのpolychorで計算してみる。

polychor(cat1,cat2)
#[1] 0.5165662

ほぼ同じ値を推定していることが分かる。母相関の区切れを細かくすればこの値に近づくだろう。という訳で、ポリコリック相関係数を自作するお話でした。

猫も杓子も構造化

発達障害、特別支援などについて書いています。最近は心理学関係の内容が多めです。

ポリコリック相関係数のお勉強

データの準備

ステップ1 閾値の推定

ステップ2 母相関の推定