ポリコリック／テトラコリック相関係数のベイズ推定

多分相関係数の簡単なベイズ推定

Simple Bayesian Estimation of Polychoric Correlation Coefficients

岡本安晴

2016.03

　カテゴリ変数の潜在変数間の相関係数は、多分相関係数（polychoric correlation coefficient）と呼ばれているが、特にカテゴリ数が２個同士の場合は四分相関係数（tetrachoric correlation coefficient）と呼ばれている。これらを簡明なモデルとアルゴリズムによって求めるプログラムを作成してみた。まず、基本モデルについて説明し、その後、一括推定、２段推定、多変量推定に分けてアルゴリズムとプログラムの使用法について説明する。

基本モデル

順序カテゴリ変量の値が正規分布に従う潜在変数によって決まるとするモデルを採用する。これは、Thurstone’s law of categorical judgment (Torgerson, 1958)、Thurstone’s method of successive intervals (Guilford, 1954) あるいは ordinal probit regression (Kruschke, 2011)と呼ばれている。２つのカテゴリ変量X と Y に対してそれぞれの潜在変数 U と V を考える。相関係数について考えるので、U と V は平均０、分散１の標準正規分布に従うとしてよい。カテゴリ変量X と Y のとるカテゴリ値は、それぞれ１からK および１から H までの整数値であるとし、潜在変数に対するカテゴリ境界値をそれぞれ

および

とおく。すなわち、

　のとき　

である。したがって、

（１）　、　

となる。は、累積標準正規分関数である。

　２変量が相関係数の２変量正規分布に従うとする。

である。ここで、

、　

である。このとき、

である。は、の確率密度関数である。上の積分は、ガウスの積分公式を適応的方法によって用いて、計算精度と計算速度の両立を計っている（岡本安晴「大学生のための心理学VC++プログラミング入門」勁草書房、pp. 114-117 参照）。

　いま、データとしての組が組与えられたとする。その内、である組が組であるとすると

となる。したがって、事後分布は

（２）　　

で与えられる。ここで、は事前分布である。

事後分布の推定を、以下のMCMCによるサンプリングによって推定する。事前分布は、有界な領域での一様分布を採用している。

一括推定

　の事後分布（２）をMCMCにより推定する。プログラムでは、Metropolis-within-Gibbsあるいはcomponentwise version of the Metropolis-Hastings algorithmと呼ばれている方法を用いている。このアルゴリズムについては、＜岡本安晴「心理学データ分析と測定」勁草書房＞の付録C「ベイズ的分析法」で説明している。プログラムのソースコードファイルなどは圧縮ファイルEstPolychoric.zipにまとめたので、クリック・ダウンロード・解凍（展開）すればVisual Studio 2015（無料版が公開されている）で開くことができる。

プログラムを開いた時のフォーム・デザイナーのエラーは無視できる。エラーの出ているフォームのタブのX印をクリックして一旦消した後、ソリューションエクスプローラ内のフォームのファイル名をクリックすればよいが、このことは次のウェブサイトで説明している。

http://y-okamoto-psy1949.la.coocan.jp/VCpp/FormError/

Visual Studioでプロジェクトを開いたときに下段にエラーメッセージが表示されることがあるが、これも無視できる。

実行形式のファイルEstPolychoric.exeはVisual Studio 2015（無料版が公開されている）のインストールされているPC上でクリック・ダウンロードして実行することができる。

Windows 8/8.1 での実行形式ファイル*.exeのダウンロード・実行については

http://y-okamoto-psy1949.la.coocan.jp/VCpp/Win8EXEFile/

で説明している。

プログラムを実行すると図1.1のフォームが提示される。

図1.1

「Go」ボタンのクリックで、図1.2の入力ファイル名設定のダイアログウィンドウが表示される。

図1.2

入力データファイルは、図1.3に示す形式で用意する。

図1.3

第１文字がスラッシュ / で始まる最初の行の次の行に2つのカテゴリ変量のカテゴリ数を設定する。図1.3の場合は、カテゴリ数として２つの変量に対してそれぞれ４が設定されている。カテゴリ数としてKが設定された変量は、１からKまでの整数値をとるものとする。

第１文字がスラッシュ / で始まる２番目の行の次の行には、変量名を設定する。図1.3では、「ID」、「X」および「Y」が設定されている。スラッシュ / で始まる３番目の行の次の行から、１行に１組のデータを書いていく。図1.3のデータでは、ID、X、Yの３つの値で１組のデータとなっている。変数IDは、データを区別するためのもので、プログラムでは文字列として読み込まれている。IDは、空白などを含まない連続した文字列である。図1.3では、「数値：」が文字列の値として使われている。このIDの値に続けて、Xのカテゴリ値とYのカテゴリ値が並べられる。

データの最後は、スラッシュ / で始まる行で終わりを示す（図1.4）。

図1.4

図1.2において、入力データファイル名の設定後、「開く」ボタンをクリックすると、図1.5の出力ファイル名の設定を求めるダイアログウィンドウが表示される。

図1.5

適当なテキストファイル名を設定して、「保存」ボタンをクリックすると計算が始まる。

最初は、図1.6に示すような「Adaptive stage」の文字列のあるフォームが表示される。

図1.6

これは、MCMCにおける提案分布や初期値の調整が行われていることを示すものである。

調整後、メインのMCMCがスタートするが、これは「Main」の文字列のあるフォームの表示によって示される。

図1.7

メインのMCMCでは４本のマルコフ鎖が同時並行で実行される。

図1.8

MCMCによる計算が終了すると、の事後分布のグラフが表示される（図1.8）。図1.8はサンプリングのスムージングによって描かれている。図1.8のように滑らかさに問題がある場合は、スムージングの幅を調節する。

図1.8では、スムージングの幅は図1.9の「kw＝」に示されている0.005である。

図1.9

これを、図1.10に示すように0.007に設定して「Draw」ボタンをクリックすると、図1.11に示すように滑らかなグラフが得られる。

図1.10

図1.11

「Save」ボタンをクリックすると、描かれている描画画面を保存することが出来る。

「Close」ボタンのクリックで、描画画面は閉じられて、図1.12のメッセージボックスが表示される。

図1.12

出力ファイル名のファイルが保存されたことが示されている。「OK」ボタンのクリックでプログラムの終了である。

保存された出力ファイルを開くと、図1.12のようになっている。

図1.13

図1.13のファイルの内容は、以下のようである。

Input Data File = D:\yasuharu\www\LaCoocan\bayesian\polychoric\EstPolychoric\EstPolychoric\Data.txt

Kx = 4

Ky = 4

Variable Names... X, Y

Data...

0: 2 3

1: 1 2

2: 1 2

3: 2 2

4: 2 2

5: 2 2

6: 3 3

7: 2 2

8: 2 2

9: 2 2

291: 4 4

292: 3 2

293: 3 2

294: 1 1

295: 4 4

296: 3 3

297: 3 4

298: 2 1

299: 3 3

Correlation coefficient = 0.61068

rho...

mean = 0.67504

Q1 = 0.65010 Med = 0.67593 Q3 = 0.70141

95% CI = [0.59240, 0.74528]

Cx[0]...

mean = -1.10732

Q1 = -1.16881 Med = -1.10600 Q3 = -1.04509

95% CI = [-1.28329, -0.93039]

Cx[1]...

mean = -0.02605

Q1 = -0.07481 Med = -0.02494 Q3 = 0.02046

95% CI = [-0.16455, 0.11427]

Cx[2]...

mean = 1.12226

Q1 = 1.06123 Med = 1.12125 Q3 = 1.18488

95% CI = [0.94253, 1.30882]

Cy[0]...

mean = -1.00482

Q1 = -1.06499 Med = -1.00288 Q3 = -0.94328

95% CI = [-1.17620, -0.83400]

Cy[1]...

mean = 0.01618

Q1 = -0.03324 Med = 0.01660 Q3 = 0.06394

95% CI = [-0.12702, 0.16108]

Cy[2]...

mean = 1.08846

Q1 = 1.02624 Med = 1.08822 Q3 = 1.14910

95% CI = [0.91725, 1.26485]

GraphRho was called with kw = 0.005

Rho_MAP = 0.6685

GraphRho was called with kw = 0.007

Rho_MAP = 0.671

「Correlation coefficient = 0.61068」は、カテゴリ値に対してそのまま相関係数を算出したときの値である。

rho...

mean = 0.67504

Q1 = 0.65010 Med = 0.67593 Q3 = 0.70141

95% CI = [0.59240, 0.74528]

は、の事後分布に対して、平均値、四分位数、および95%確信区間を算出したものである。

Cx[k] は、に対するもの、Cy[h] は、に対する事後分布から算出された値を示している。

「GraphRho」に続く

GraphRho was called with kw = 0.005

Rho_MAP = 0.6685

GraphRho was called with kw = 0.007

Rho_MAP = 0.671

は、k = 0.005 および k=0.007 に対して求められたMAPの値である。図1.8あるいは図1.11に示さているグラフから適切なものを選べばよい。

２段推定

　パラメータの推定とパラメータの推定の２段に分けて行う。

第１段

パラメータの推定を式（１）によって行う。

第２段

第１段でのの推定値に基づいて、事後分布（２）をのみの関数とみての事後分布の推定をMCMCによって行う。

２段推定法のプログラムも作成した。ソースコードファイルなどは圧縮ファイルEstPolychoricTwoStep.zipとしてまとめたので、クリック・ダウンロード・展開すればVisual Studio 2015で開くことができる。実行形式のファイルEstPolychoricTwoStep.exeは、Visual Studio 2015のインストールされているPC上で、クリック・ダウンロードして実行することができる。

プログラムの実行は、ほぼ上の一括推定のプログラムの場合と同様である。

出力ファイルの例を以下に示す。

Input Data File = D:\yasuharu\www\LaCoocan\bayesian\polychoric\EstPolychoricTwoStep\EstPolychoricTwoStep\Data.txt

Kx = 4

Ky = 4

Variable Names... X, Y

Data...

0: 2 3

1: 1 2

2: 1 2

3: 2 2

4: 2 2

5: 2 2

6: 3 3

7: 2 2

8: 2 2

9: 2 2

291: 4 4

292: 3 2

293: 3 2

294: 1 1

295: 4 4

296: 3 3

297: 3 4

298: 2 1

299: 3 3

Correlation coefficient = 0.61068

Cumx[0] = 0.13667

Cumx[1] = 0.48667

Cumx[2] = 0.86667

Cumy[0] = 0.16000

Cumy[1] = 0.50333

Cumy[2] = 0.86000

Cx[0] = -1.09542

Cx[1] = -0.03343

Cx[2] = 1.11077

Cy[0] = -0.99446

Cy[1] = 0.00836

Cy[2] = 1.08032

rho...

mean = 0.67939

Q1 = 0.65654 Med = 0.68180 Q3 = 0.70475

95% CI = [0.60534, 0.74301]

GraphRho was called with kw = 0.005

Rho_MAP = 0.68495

「Cumx[k]」あるいは「Cumy[h]」は、「」あるいは「」の累積比率である。

「Cx[k]」あるいは「Cy[h]」は、対応するカテゴリ境界値を示している。

以下の、

rho...

mean = 0.67939

Q1 = 0.65654 Med = 0.68180 Q3 = 0.70475

95% CI = [0.60534, 0.74301]

GraphRho was called with kw = 0.005

Rho_MAP = 0.68495

は、上の一括推定の場合と同じである。

多変量推定

　カテゴリ変量が、・・・、のM個の場合についてポリコリック相関係数を求める。このとき、カテゴリ変量の同時分布はM次元正規分布で与えられるが、この場合の多重積分は実用的ではない。２個ずつの対に対してポリコリック相関係数を求める。上の２段推定を用いる。

第１段

　カテゴリ変量のカテゴリ境界を、・・・、とおく。各に対して式（１）によってを推定する。

第２段

　各カテゴリ変量の対に対して式（２）によるベイズ推定をMCMCによって行い、事後分布を推定する。

　上の方法による推定を行うプログラムを作成した。ソースコードファイルなどは圧縮ファイルEstPolychoricMV.zipとしてまとめた。クリックしてダウンロード・展開すれば、Visual Studio 2015で開くことができる。実行形式のファイルEstPolychoricMV.exeは、Visual Studio 2015のインストールされているPC上で、クリック・ダウンロード・実行することができる。

プログラムを実行すると、図3.1のフォームが表示される。

図3.1

「Go」ボタンのクリックで、図3.2の入力ファイル名設定のダイアログウィンドウが表示される。

図3.2

入力データファイルは、図3.2に示す形式で用意する。

図3.3

第１文字目がスラッシュ / で始まる最初の行の次の行に、カテゴリ変数の数を設定する。

２番目のスラッシュ / で始まる行の次の行に、各カテゴリ変数のカテゴリ数を書く。カテゴリ数がKであれば、そのカテゴリ変数は、１～Kの範囲の整数値をとるものとする。

３番目のスラッシュ / で始まる行の次の行には、変数のラベルを並べる。先頭の変数は、データの識別変数である。

４番目のスラッシュ / で始まる行に続けて、各ケースのデータを１組ずつ書く。第１番目の変数はIDであるが、これは文字列（ただし、空白などの区切り文字を含まないもの）の値として読み込まれる。図3.2では、「数値：」がIDの文字列として使われている。

データの最後は、スラッシュ / で始まる行をおいてデータの終わりであることを示す（図3.4）。

図3.4

入力データファイル名を設定して、「開く」ボタンをクリックすると、出力ファイル名の設定を求めるダイアログウィンドウが表示される（図3.5）。

図3.5

出力ファイル名の設定後、「保存」ボタンのクリックで計算が始まる。Adaptive stageのフォームの提示後、メインMCMCのフォームが表示されることは、上の一括推定の場合と同様であるが、図3.6の場合は、どの組み合わせのカテゴリ変数に対して推定が行われているかの情報も示されている

図3.6

カテゴリ変数はM個あれば、０からまでの番号を付けて表されている。すべての組み合わせについて計算が終了すると、図3.7のメッセージボックスが表示される。

図3.7

出力ファイルの名前を確認して、「OK」ボタンをクリックするとプログラムは終了する。

出力ファイルを開くと、図3.8のようになっている。

図3.8

図3.8のファイルは以下のような内容である。

Input Data File = D:\yasuharu\www\LaCoocan\bayesian\polychoric\EstPolychoricMV\EstPolychoricMV\DataHeteroN10000.txt

Number of variables = 4

Variable Names...

X1 X2 X3 X4

Number of categories...

3 3 5 5

Data...

0: 1 2 2 3

1: 2 2 2 2

2: 2 2 4 4

3: 2 2 3 3

4: 2 2 4 4

5: 2 2 4 3

6: 1 2 3 2

7: 2 3 4 4

8: 1 2 3 2

9: 2 2 3 4

9991: 2 2 3 3

9992: 1 2 1 1

9993: 1 2 1 2

9994: 2 2 3 3

9995: 2 2 4 2

9996: 2 2 4 3

9997: 2 3 4 4

9998: 2 2 2 5

9999: 2 2 3 3

r =

1.00000 0.52124 0.57227 0.56587

0.52124 1.00000 0.56819 0.57009

0.57227 0.56819 1.00000 0.63428

0.56587 0.57009 0.63428 1.00000

Cumx[0][0] = 0.15720

Cumx[0][1] = 0.84280

Cumx[1][0] = 0.14760

Cumx[1][1] = 0.84560

Cumx[2][0] = 0.06370

Cumx[2][1] = 0.31120

Cumx[2][2] = 0.68970

Cumx[2][3] = 0.93290

Cumx[3][0] = 0.06380

Cumx[3][1] = 0.30350

Cumx[3][2] = 0.69550

Cumx[3][3] = 0.93290

Cx[0][0] = -1.00603

Cx[0][1] = 1.00603

Cx[1][0] = -1.04678

Cx[1][1] = 1.01774

Cx[2][0] = -1.52444

Cx[2][1] = -0.49245

Cx[2][2] = 0.49500

Cx[2][3] = 1.49774

Cx[3][0] = -1.52363

Cx[3][1] = -0.51436

Cx[3][2] = 0.51150

Cx[3][3] = 1.49774

Estimation for X1 and X2 started.

rho...

mean = 0.69508

Q1 = 0.68964 Med = 0.69495 Q3 = 0.70066

95% CI = [0.67872, 0.71104]

Estimation for X1 and X3 started.

rho...

mean = 0.69624

Q1 = 0.69175 Med = 0.69642 Q3 = 0.70086

95% CI = [0.68235, 0.70964]

Estimation for X1 and X4 started.

rho...

mean = 0.68747

Q1 = 0.68293 Med = 0.68759 Q3 = 0.69206

95% CI = [0.67393, 0.70086]

Estimation for X2 and X3 started.

rho...

mean = 0.69798

Q1 = 0.69361 Med = 0.69812 Q3 = 0.70259

95% CI = [0.68452, 0.71067]

Estimation for X2 and X4 started.

rho...

mean = 0.69963

Q1 = 0.69516 Med = 0.69978 Q3 = 0.70436

95% CI = [0.68533, 0.71272]

Estimation for X3 and X4 started.

rho...

mean = 0.69372

Q1 = 0.69016 Med = 0.69383 Q3 = 0.69745

95% CI = [0.68247, 0.70412]

r (means) ...

X1 X2 X3 X4

X1 1.00000 0.69508 0.69624 0.68747

X2 0.69508 1.00000 0.69798 0.69963

X3 0.69624 0.69798 1.00000 0.69372

X4 0.68747 0.69963 0.69372 1.00000

r (medians) ...

X1 X2 X3 X4

X1 1.00000 0.69495 0.69642 0.68759

X2 0.69495 1.00000 0.69812 0.69978

X3 0.69642 0.69812 1.00000 0.69383

X4 0.68759 0.69978 0.69383 1.00000

r =

1.00000 0.52124 0.57227 0.56587

0.52124 1.00000 0.56819 0.57009

0.57227 0.56819 1.00000 0.63428

0.56587 0.57009 0.63428 1.00000

は、カテゴリ変量の値から直接相関係数を計算したものである。

「Cumx[j][k]」は、ｊ番目のカテゴリ変数に対するカテゴリの累積比率であり、「Cx[j][k]」はそのカテゴリ境界値である。

Estimation for X1 and X2 started.

rho...

mean = 0.69508

Q1 = 0.68964 Med = 0.69495 Q3 = 0.70066

95% CI = [0.67872, 0.71104]

は、カテゴリ変数X1とX2に対するの事後分布の推定であることを示す。

r (means) ...

X1 X2 X3 X4

X1 1.00000 0.69508 0.69624 0.68747

X2 0.69508 1.00000 0.69798 0.69963

X3 0.69624 0.69798 1.00000 0.69372

X4 0.68747 0.69963 0.69372 1.00000

は、事後分布の平均値を並べたものである。

r (medians) ...

X1 X2 X3 X4

X1 1.00000 0.69495 0.69642 0.68759

X2 0.69495 1.00000 0.69812 0.69978

X3 0.69642 0.69812 1.00000 0.69383

X4 0.68759 0.69978 0.69383 1.00000

は、事後分布の中央値を並べたものである。

なお、相関行列（相関係数を行列の形に並べたもの）に基づいて因子分析を行うプログラムについては別にウェブサイトを用意して扱っている。

引用文献

岡本安晴（2014）「心理学データ分析と測定」勁草書房

岡本安晴（2013）「２値データの相関係数のベイズ推定」日本心理学会第７７回大会発表論文集、ポスター、1AM-070

岡本安晴（2010）「大学生のための心理学VC++プログラミング入門」勁草書房