多分相関係数／四分相関係数の簡単なベイズ推定
Pythonスクリプト
岡本安晴、2018
カテゴリ変数の潜在変数間の相関係数は、多分相関係数（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によるサンプリングによって推定する。事前分布は、有界な領域での一様分布を採用している。ただし、以下のプログラムでは、パラメータ
の推定とパラメータ
の推定の２段に分けて行っている。２段に分ける方法については、D. F. Alwin (2007) “Margins of Error” A John Wiley &
Sons, Inc., p. 45で論じられている。
第１段
パラメータの推定を式（１）によって行う。
第２段
第１段でのの推定値に基づいて、事後分布（２）をのみの関数とみての事後分布の推定をMCMCによって行う。
の事後分布（２）をMCMCにより推定する。プログラムでは、Metropolis-Hastings
algorithmと呼ばれている方法を用いている。このアルゴリズムについては、＜岡本安晴「心理学データ分析と測定」勁草書房＞の付録C「ベイズ的分析法」で説明している。Pythonスクリプトを以下に示す。
den_t = None
ck_den = 0
for t in range(0, NSimu):
if
(t % 100 == 0):
print(MsgStr, ' t = ', t,
'/', NSimu)
cnvs.create_rectangle(0, 0, 500, 150, outline = '#FFFFFF', fill =
'#FFFFFF')
cnvs.create_text(250, 75,
text = '{0} t =
{1}/{2}'.format(MsgStr, t, NSimu),
font = ('', 15), fill = '#000000')
cnvs.update()
num
= None
den
= None
y =
rn.normalMS(rhos[t], gen_rho)
a =
None
if
(math.fabs(y) >= 1.0):
a = -1.0
else:
num = LogL(Nij, Kx, Ky, Cx, Cy, y)
if ck_den == 0:
den = LogL(Nij, Kx, Ky, Cx, Cy, rhos[t])
else:
den = den_t
a = myExp(num - den)
if
rn.uni() < a:
rhos.append(y)
n_acpt += 1
ck_den += 1
den_t = num
else:
rhos.append(rhos[t])
上のスクリプトは、モジュールmcmc.pyの一部を示したものであるが、このモジュールを含めた全てのファイルとサンプルデータファイルは圧縮ファイルPolychoric.zipにまとめた。ダウンロードして自由に使用されたい。
主スクリプトファイルはMain.pyである。実行すると入力データファイル名と出力ファイル名を聞いてくる（図１）。
図１
入力データファイは、図２に示すように、スラッシュ / で始まる行で区切られている。最初のスラッシュ / で始まる行の次の行に、カテゴリ変数のカテゴリ数が書かれている。次のスラッシュで始まる行の次の行に変数のラベルが書かれている。変数のラベルは、各データのID用変数（例えば、通し番号）とカテゴリ変数に対するものである。図２のファイルでは、ID、X、およびYがラベルとして書かれている。３つ目のスラッシュで始まる行に続いて、データが「ケースID 第1カテゴリ変数 第2カテゴリ変数」の形式で１行に１組ずつのデータが書かれている。ケースIDは文字列として扱われるので、通し番号でも名前でもよい。
図２
データの最後は、スラッシュで始まる行を置いて示される（図３）。
図３
入力ファイル名と出力ファイル名の設定が終わると、計算が始まる（図４）。
図４
MCMCをデータに合わせて調整した後、主MCMCがスタートする。主MCMCが終了すると、相関係数
の事後分布が変数rhoに対するものとしてヒストグラムの形で表示される（図５）。
図５
ヒストグラムのウィンドウの右上角のX印をクリックして閉じると、プログラムの実行終了となる。
プログラムの実行終了後、出力ファイルを開くと図６のように出力されている。
図６
２つのカテゴリ変量に対するカテゴリ境界がリストCxとCyの要素として書き出されている。相関係数の事後分布の四分位数がQ1（第１四分位数）、Med（中央値）、Q3（第３四分位数）の値として書き出されている。
カテゴリ変量が
、・・・、
のM個の場合についてポリコリック相関係数（カテゴリ数が２個の場合は、テトラコリック相関係数）を求める。このとき、カテゴリ変量の同時分布はM次元正規分布で与えられるが、この場合の多重積分は実用的ではない。２個ずつの対に対してポリコリック相関係数を求める。上の２段推定を用いる。すなわち、M個の変数から２値ずつを順番に選び、上の２変量の方法を適用した。プログラムとサンプルデータファイルは圧縮ファイルPolychoricMV.zipとしてまとめた。ダウンロードして自由に利用することができる。
プログラムの使用法は、２変量の場合とほぼ同じである。主スクリプトファイルMain.pyを実行すると、入力データファイル名と出力ファイル名の設定が求められる。
図７
入力データファイルは、図８のような形式で用意する。
図８
行の先頭がスラッシュ / で始まる行を区切り行として用いる。１番目の行がスラッシュで始まる行の次の行に変数（変量）の数を書く（図８では、４が書かれている）。
２番目のスラッシュで始まる行の次の行に、各カテゴリ変量のカテゴリ数を書く。図８では、「3 3 5 5」が書かれている。カテゴリ変量のカテゴリ数がKのときは、カテゴリ変量の値は、１からKまでの整数値で表すものとする。データにおいて、用いられていないカテゴリがあれば、そのカテゴリは除く。例えば、K個のカテゴリのうち、用いられていない（回答のない）カテゴリが１つあれば、そのカテゴリを除いて考え、K−1個として扱う。このとき、カテゴリ値を表す整数値も、１からｋ−１になるように注意し、必要ならば変更する。
３番目のスラッシュで始まる行の次の行に、変数のラベル（名前）を並べる。最初に、各データ（ケース）の識別文字列を値とする変量のラベルを書き、その後に、カテゴリ変量のラベルを順番に並べる。図８では、IDはケース識別変量のラベルであり、カテゴリ変量のラベルはX1からX4である。
４番目のスラッシュで始まる行の次の行から、１行に１つのケースのデータを書いていく。各行の１番目の値は、ケース識別用文字列であり、図８では０から始まる通し番号が用いられている。ケース識別用文字列に続いて、各カテゴリ変量の値が書かれている。
最後のデータの次の行は、スラッシュで始まる行を置いて、データの終わりであることが示されている（図９）。
図９
入力データファイル名と出力ファイル名の設定が終わると、計算が始まる。MCMCによるサンプリングの経過は、２変量の場合と同じように表示用フォームが提示されて示される（図１0）。
図１０
M個のカテゴリ変数は、変数を０番からの整数で番号付けしたとき、変数―０と変数―１、変数―０と変数―２、という具合に２つの組み合わせが順に選ばれて計算され、変数―（M-2）と変数―（M-1）の組み合わせの計算が終了したところで全体の計算終了となる（図１１）。２変量の場合と異なり、各組合せにおけるの事後分布のヒストグラムの提示はない。
図１１
計算終了後、出力ファイルを開くと、ファイルの最後にポリコリック相関係数（テトラコリック相関係数）が行列の形で出力されている（図１２）。
図１２
引用文献
岡本安晴（2013）「２値データの相関係数のベイズ推定」日本心理学会第７７回大会発表論文集、ポスター、1AM-070
岡本安晴（2010）「大学生のための心理学VC++プログラミング入門」勁草書房