潜在クラス項目反応理論
発達段階理論による尺度構成モデル
現代テスト理論の項目反応理論（Item Response Model: IRT）では、正答確率を特性値の関数として
（１）
と表す。ここで、特性値
は連続変数である。しかし、Piagetの発達段階理論に見られるように、心理学的概念あるいは特性は心理特性の段階を表すカテゴリ変数の場合もある。このとき、正答確率は心理特性の段階の関数として
（２）
と表される。生澤はモデル（２）を用いた研究を報告しており、学童の興味発達の潜在構造分析（生澤、1965）、知能発達の基本構造（生澤、1972）、潜在クラス分析におけるGreen解の研究（生澤、1974）を上げることができる。
モデル（２）による確率モデルのベイズ分析を岡本（2008、2009、2013）は報告しているが、Metropolis-Hasting within GibbsアルゴリズムをC++で実装したものが用いられている。
本ウェブサイトでは、ベイズ分析を行うスクリプトをPyMC3で開発してみたところ、以下に示すように簡単なスクリプトになることが分かった。本ウェブサイトのスクリプトはUbuntuで開発したが、Windowsパソコンの場合はJupyter Notebookで用いることができる（2020.03現在）。
スクリプトを示す前に、モデル（２）について少し説明を加えておく。
まず、モデル（２）で表される正答確率は、段階が進むにつれて増加する条件Stageの単調増加関数である。
（３） Stage k -> Stage (k+1) のとき、 P(Correct|Stage
k) ≤ P(Correct|Stage k+1)
条件（３）の単調性の条件を置かないとき、すなわち正答確率を単に潜在クラスのみの関数とするときは、MCMCによるサンプリングにおいてラベルスウィッチングが起きないように適当な制約等を置く必要があることに注意。
いま、Stageの数がStage 1からStage KまでのK段階であるとする。各段階ごとの正答確率の増加分をp_stepとおく。すなわち、
p_step[1]=P(Correct|Stage 1)
p_step[k]=P(Correct|Stage
k)-P(Correct|Stage k-1)
p_step[K+1]=1-P(Correct|Stage K)
とおく。このとき、
P(Correct|Stage k)=p_step[1]+…+p_step[k]
p_step[k] ≥ 0
p_step[1]+…+p_step[K+1]=1
である。確率P(Correct|Stage k)とp_step[k]の関係は、行列を用いて以下のように表せる（例として、K=3の場合）。
これは、以下に示すスクリプトでは次のようにprob_reponseとして与えられている。
prob_response_
= pm.math.matrix_dot(coeff, prob_steps.T)
prob_stepsはp_stepを表すもので、Dirichlet型のオブジェクトとして
prob_steps
= pm.Dirichlet('prob_steps', a = a_Diri, shape = (M, NStg + 1))
と宣言されている。NStgはカテゴリ数K、Mは項目数である。
段階理論あるいは潜在クラスモデルでは、各回答者の属するクラスはカテゴリ変数である。以下に示すスクリプトでは、Categorical型のオブジェクトとして
p_stages
= pm.Categorical('p_stages', p = a_p, shape = N)
と宣言され、p_stagesで表されている。Nは回答者数である。
PyMC3はカテゴリ変数も扱うことができるが、Stanではパラメータとしてカテゴリ変数を選定することはできない。これに対応するため、パラメータの確率の周辺化を行う。Stanによるスクリプト例は別のウェブサイトで説明している。
上で説明したモデルのサンプルスクリプトをリスト１に示すが、スクリプトファイルは圧縮ファイルIRTLClassPyMC3.zipとして用意した。また、データをファイルから入力する一般的な使用を想定したサンプルスクリプトは、このウェブサイトの後半にリスト２として用意し、圧縮ファイルIRTLClassPyMC3RD.zipとしても用意した。PyMC3の説明は＜岡本安晴「いまさら聞けないPythonでデータ分析――多変量解析、ベイズ分析（PyStan、PyMC）――」丸善出版＞で行っている。
なお、PyMC3をWindowsで利用する場合の説明を、このウェブサイトで行っている。
リスト１ 発達段階理論用尺度構成モデルのデモ用サンプルスクリプト
"""
Yasuharu Okamoto, 2019.02
"""
import
matplotlib.pyplot as plt
import numpy
as np
import pymc3
as pm
import arviz
as az
X = ([[0]*5
+ [0]*5] * 7) + ([[1]*5 + [0]*5] * 7) +\
([[1]*5
+ [1]*5] * 7)
print('X =')
for v in X:
print(v)
N = len(X)
print('N = ',
N)
M =
len(X[0])
print('M =
', M)
NStg = 3 # Number of
stages
print('NStg
= ', NStg)
a_Diri =
np.ones((M, NStg + 1))
print('a_Diri
=\n', a_Diri)
a_coeff =
np.zeros((NStg, NStg + 1))
for i in
range(NStg):
for j
in range(i + 1):
a_coeff[i][j]
= 1
print('a_coeff
=\n', a_coeff)
a_p =
np.ones(NStg)
print('a_p
=\n', a_p)
with
pm.Model() as stage_model:
prob_steps
= pm.Dirichlet('prob_steps', a = a_Diri, shape = (M, NStg + 1))
coeff
= pm.math.constant(a_coeff, 'coeff')
prob_response_
= pm.math.matrix_dot(coeff, prob_steps.T)
prob_response
= pm.Deterministic('prob_response', prob_response_)
p_stages
= pm.Categorical('p_stages', p = a_p, shape = N)
p_at_stage
= pm.Deterministic('p_at_stage', prob_response[p_stages])
y =
pm.Bernoulli('y', p = prob_response[p_stages], observed = X)
trace
= pm.sample()
summary =
pm.summary(trace)
print('summary...\n',
summary)
ary_prob_response
= trace['prob_response']
print('\nP(correct|item,stage)\n')
print('{0:
>10s}'.format(' '), end = '')
for i in range(NStg):
print("{0:>10s}".format('Stage-{}'.format(i)),
end = '')
print('
')
for j in
range(M):
print('Item-{0:
<5d}'.format(j), end = '')
for s
in range(NStg):
print('{0:
>10.5f}'.format(ary_prob_response[:,s,j].mean()), end = '')
print('
')
ary_p_stages
= trace['p_stages']
ary_p_stages
= ary_p_stages.T
print("\nMean
stages and their SD's\n")
print('{0:
>11s}{1:>10s}{2:>10s}'.format(' ', 'Mean', 'SD'))
for i in
range(len(ary_p_stages)):
print('person-{0:
<4d}{1:>10.3f}{2:>10.3f}'.format(i, ary_p_stages[i].mean(),
ary_p_stages[i].std()))
v_waic =
az.waic(trace, scale = 'deviance')
print(v_waic)
print('\nWAIC
= {0:.3f} for the number of stages
= {1}'.format(v_waic['waic'], NStg))
リスト１のスクリプトを実行すると、各回答者の反応データが次のように出力される。
X =
[0, 0, 0, 0,
0, 0, 0, 0, 0, 0]
[0, 0, 0, 0,
0, 0, 0, 0, 0, 0]
[0, 0, 0, 0,
0, 0, 0, 0, 0, 0]
[0, 0, 0, 0,
0, 0, 0, 0, 0, 0]
[0, 0, 0, 0,
0, 0, 0, 0, 0, 0]
[0, 0, 0, 0,
0, 0, 0, 0, 0, 0]
[0, 0, 0, 0,
0, 0, 0, 0, 0, 0]
[1, 1, 1, 1,
1, 0, 0, 0, 0, 0]
[1, 1, 1, 1,
1, 0, 0, 0, 0, 0]
[1, 1, 1, 1,
1, 0, 0, 0, 0, 0]
[1, 1, 1, 1,
1, 0, 0, 0, 0, 0]
[1, 1, 1, 1,
1, 0, 0, 0, 0, 0]
[1, 1, 1, 1,
1, 0, 0, 0, 0, 0]
[1, 1, 1, 1,
1, 0, 0, 0, 0, 0]
[1, 1, 1, 1,
1, 1, 1, 1, 1, 1]
[1, 1, 1, 1,
1, 1, 1, 1, 1, 1]
[1, 1, 1, 1,
1, 1, 1, 1, 1, 1]
[1, 1, 1, 1,
1, 1, 1, 1, 1, 1]
[1, 1, 1, 1,
1, 1, 1, 1, 1, 1]
[1, 1, 1, 1,
1, 1, 1, 1, 1, 1]
[1, 1, 1, 1,
1, 1, 1, 1, 1, 1]
１０項目、２１人分のデータであり、各人のデータが各１行に表示される。正答は1で、誤答は0で表されている。
PyMC3によるサンプリングが終了すると、パラメータのサンプリング値から以下の数値が出力される。
まず、
の平均値が、項目ごとに各ステージにおける値として以下のように表示される。
P(correct|item,stage)
Stage-0 Stage-1 Stage-2
Item-0
0.11299 0.83643 0.94046
Item-1
0.10882 0.83649 0.94223
Item-2
0.11266 0.83717 0.94057
Item-3 0.12034 0.83957 0.94359
Item-4
0.10975 0.83419 0.93951
Item-5
0.06003 0.16121 0.89432
Item-6
0.05926 0.16450 0.88910
Item-7
0.05962 0.16370 0.88480
Item-8
0.05887 0.16182 0.88767
Item-9 0.05980 0.16246 0.88645
ステージ０のときは、どの項目も正答確率は0.1前後と低く、ステージ１では、Item-0からItem-4までは正答確率は0.8に近いが、Item-5からItem-9までは約0.2で低い。ステージ２になるとどの項目の正答確率も約0.9と高い値になる。
次に各回答者がどのステージにいるかを表すp_stagesの平均値がステージごとに以下のように表示される。
Mean stages
and their SD's
Mean SD
person-0
0.000
0.000
person-1
0.000
0.000
person-2
0.000
0.000
person-3
0.000
0.000
person-4
0.000
0.000
person-5
0.000
0.000
person-6
0.003
0.055
person-7
1.000 0.000
person-8
1.000
0.000
person-9
1.000
0.000
person-10 1.001 0.032
person-11 1.000 0.000
person-12 1.000 0.000
person-13 1.000 0.000
person-14 2.000 0.000
person-15 2.000 0.000
person-16 2.000 0.000
person-17 2.000 0.000
person-18 2.000 0.000
person-19 2.000 0.000
person-20 2.000 0.000
Person-0からPerson-6まではステージ０であり、Person-7からPerson-13まではステージ１であり、Person-14からPerson-20までは、ステージ２であることがわかる。標準偏差SDが実質的に０であることは、サンプリングにおけるトレースがほとんど変化していない、すなわち分散が０であることを意味する。このとき、n_effあるいはRhatは計算不能でNaNとなる。サンプリングの終了後トレースのsummaryを表示すると、以下のようになっている。
Multiprocess
sampling (2 chains in 2 jobs)
CompoundStep
>NUTS:
[prob_steps]
>CategoricalGibbsMetropolis:
[p_stages]
Sampling 2
chains, 0 divergences: 100%|█| 2000/2000 [00:06<00:00,
316.30draws/s]
summary...
mean sd hpd_3% ...
ess_bulk ess_tail r_hat
p_stages[0] 0.000 0.000 0.000 ... 1000.0 1000.0 NaN
p_stages[1] 0.000 0.000 0.000 ... 1000.0 1000.0 NaN
p_stages[2] 0.000 0.000 0.000 ... 1000.0 1000.0 NaN
p_stages[3] 0.000 0.000 0.000 ... 1000.0 1000.0 NaN
p_stages[4] 0.000 0.000 0.000 ... 1000.0 1000.0 NaN
...
... ... ... ... ... ... ...
p_at_stage[20,5] 0.887 0.102 0.707 ... 625.0 430.0 1.0
p_at_stage[20,6] 0.888 0.096 0.719 ... 1307.0 575.0 1.0
p_at_stage[20,7] 0.892 0.095 0.720 ... 1020.0 494.0 1.0
p_at_stage[20,8] 0.888 0.099 0.709 ... 1061.0 417.0 1.0
p_at_stage[20,9] 0.892 0.092 0.723 ... 1029.0 693.0 1.0
[301 rows x
11 columns]
リスト１のデモスクリプトでは、データXはプログラム内で生成されているが、実データを分析する場合は、データをファイルとして用意してプログラムで読み込めるようにしておくと使い易い。CSV形式で用意されたデータを読み込むスクリプト例を作成した（リスト２）。CSV形式のファイルはExcelなどで作成できる。図１は、上のリスト１におけるデータXをLibreOffice Calcに入力したものである。１行目に変数名を設定して、２行目以降に１行に1人のデータを書いていく。第１欄に回答者のIDとしての文字列を設定し（数値は実行時にエラーとなる）、第２欄目以降に各項目の解答の正誤を記入する。正答は１、誤答は０である。
図１
ファイルの保存のときに、「Save As...」を選ぶと図２のダイアログボックスが表示される。
図２
右下の赤丸で囲まれた下三角印をクリックすると様々な形式が提示される（図３）。提示されたものからText CSV(.csv)を選ぶと、図２のようになる。
図３
リスト２にデータをファイルから読み込むスクリプト例を示す。スクリプトファイルとサンプルデータファイルは圧縮ファイルIRTLClassPyMC3RD.zipとして用意した。
リスト2 データをファイルから読み込むスクリプト例
"""
Yasuharu Okamoto, 2019.02, 2019.03
"""
import csv
import
matplotlib.pyplot as plt
import numpy
as np
import pymc3
as pm
import arviz
as az
fin_name =
input('Data file name(*.csv) = ')
with
open(fin_name, 'r') as f:
data = [v for v in
csv.reader(f)]
fout_name =
input('Output file name(*.txt) = ')
fout =
open(fout_name, 'w')
fout.write('Data
file name = {}\n'.format(fin_name))
print('data
= \n', data)
N =
len(data) - 1
M =
len(data[0]) - 1
Items = []
for j in range(M):
Items.append(data[0][j + 1])
X =
np.empty((N, M), dtype = int)
CaseID = []
for i in
range(1, N + 1):
CaseID.append(data[i][0])
for j in range(M):
X[i
- 1][j] = int(data[i][j + 1])
print('X =')
fout.write('\nX
=\n')
i = 0
for v in X:
print(CaseID[i], ' ',end =
'')
fout.write('{}
'.format(CaseID[i]))
i += 1
print(v)
fout.write('{}\n'.format(v))
N = len(X)
print('N =
', N)
fout.write('\nN
= {}\n'.format(N))
M =
len(X[0])
print('M =
', M)
fout.write('M
= {}\n'.format(M))
NStg =
int(input('Number of stages = '))
print('NStg
= ', NStg)
a_Diri =
np.ones((M, NStg + 1))
print('a_Diri
=\n', a_Diri)
a_coeff =
np.zeros((NStg, NStg + 1))
for i in
range(NStg):
for j
in range(i + 1):
a_coeff[i][j]
= 1
print('a_coeff
=\n', a_coeff)
a_p =
np.ones(NStg)
print('a_p
=\n', a_p)
with
pm.Model() as stage_model:
prob_steps
= pm.Dirichlet('prob_steps', a = a_Diri, shape = (M, NStg + 1))
coeff
= pm.math.constant(a_coeff, 'coeff')
prob_response_
= pm.math.matrix_dot(coeff, prob_steps.T)
prob_response
= pm.Deterministic('prob_response', prob_response_)
p_stages
= pm.Categorical('p_stages', p = a_p, shape = N)
p_at_stage
= pm.Deterministic('p_at_stage', prob_response[p_stages])
y =
pm.Bernoulli('y', p = prob_response[p_stages], observed = X)
trace
= pm.sample()
summary =
pm.summary(trace)
print('summary...\n',
summary)
fout.write('summary...\n{}'.format(summary))
ary_prob_response
= trace['prob_response']
print('\n\nP(correct|item,stage)\n')
fout.write('\n\nP(correct|item,stage)\n\n')
print('{0:
>10s}'.format(' '), end = '')
fout.write('{0:
>10s}'.format(' '))
for i in
range(NStg):
print("{0:>10s}".format('Stage-{}'.format(i)),
end = '')
fout.write("{0:>10s}".format('Stage-{}'.format(i)))
print('
')
fout.write('\n')
for j in
range(M):
print('{0:>10s}'.format(Items[j]),
end = '')
fout.write('{0:>10s}'.format(Items[j]))
for s
in range(NStg):
print('{0:
>10.5f}'.format(ary_prob_response[:,s,j].mean()), end = '')
fout.write('{0:
>10.5f}'.format(ary_prob_response[:,s,j].mean()))
print('
')
fout.write('\n')
ary_p_stages
= trace['p_stages']
ary_p_stages
= ary_p_stages.T
print("\nMean
stages and their SD's\n")
fout.write("\nMean
stages and their SD's\n\n")
print('{0:
>11s}{1:>10s}{2:>10s}'.format(' ', 'Mean', 'SD'))
fout.write('{0:
>11s}{1:>10s}{2:>10s}\n'.format(' ', 'Mean', 'SD'))
for i in
range(len(ary_p_stages)):
print('{0:>11s}{1:>10.3f}{2:>10.3f}'.format(CaseID[i],
ary_p_stages[i].mean(),
ary_p_stages[i].std()))
fout.write('{0:>11s}{1:>10.3f}{2:>10.3f}\n'.format(CaseID[i],
ary_p_stages[i].mean(),
ary_p_stages[i].std()))
v_waic =
az.waic(trace, scale = 'deviance')
print('\nWAIC
statistics...\n', v_waic, '\n')
fout.write('\nWAIC
statistics...\n{}\n\n'.format(v_waic))
print('\nWAIC
= {0:.3f} for the number of stages
= {1}'.format(v_waic['waic'], NStg))
fout.write('\nWAIC
= {0:.3f} for the number of stages
= {1}\n'.format(v_waic['waic'], NStg))
print('Output
file name =', fout_name)
fout.close()
上のスクリプト（リスト２）を実行すると、入力データファイル名と計算結果の出力用ファイル名の設定が求められる（図４）。
図４
図４では、入力データファイル名として図１のLibreOffice Calcで作成されたCSV形式のファイルData01.csvが設定され、計算結果の出力用ファイル名としてResults.txtが設定されている。
出力用ファイル名を設定するとデータの読み込みが始まり、回答者数Nと項目数Mの出力後、ステージ数（カテゴリー数）の設定が求められる（図５）。
図５
図５では、ステージ数として３が設定されている。ステージ数の設定後、PyMC3によるサンプリングが始まる。サンプリング終了後、計算結果がリスト１の場合と同様に以下のように出力される。上のリスト１の結果と完全には一致しないのは、ランダムサンプリングに基づく推定だからである。
P(correct|item,stage)
Stage-0 Stage-1 Stage-2
X1 0.10897 0.83819 0.94009
X2 0.10898 0.83755 0.94136
X3 0.10994 0.83308 0.93881
X4 0.11381 0.83792 0.94322
X5 0.10874 0.83533 0.94079
X6 0.05724 0.16366 0.88769
X7 0.06045 0.16525 0.88490
X8 0.06016 0.16684 0.88858
X9 0.06126 0.17013 0.88683
X10 0.05933 0.16181 0.88655
Mean stages
and their SD's
Mean SD
psn1 0.000 0.000
psn2 0.000 0.000
psn3 0.000 0.000
psn4 0.000 0.000
psn5 0.007 0.083
psn6 0.000 0.000
psn7 0.000 0.000
psn8 1.002 0.045
psn9 1.000 0.000
psn10 1.000 0.000
psn11 1.000 0.000
psn12 1.000 0.000
psn13 1.000 0.000
psn14 1.000 0.000
psn15 2.000 0.000
psn16 2.000 0.000
psn17 1.995 0.071
psn18 2.000 0.000
psn19 2.000 0.000
psn20 2.000 0.000
psn21 2.000 0.000
WAIC
statistics...
Computed from 1000 by 210 log-likelihood
matrix
Estimate SE
deviance_waic 55.72 1.58
p_waic
2.79 -
WAIC =
55.724 for the number of stages = 3
Output file
name = ResultsRD.txt
モデル適合度の指標WAICが出力されているので、モデルの比較、すなわち段階数（カテゴリ数）を変えて適合度の比較ができる。
出力の最後に、出力ファイル名（上の場合は、Results.txt）が出力されている。出力ファイルは、計算終了後、開いて見ることができる。
引用文献
生澤雅夫（1965）．学童の興味発達の潜在構造分析法による研究．教育心理学研究、13、6−17．
生澤雅夫（1972）．知能発達の基本構造．学位論文（京都大学）．
生澤雅夫（1974）．潜在クラス分析におけるGreen解と反応反転．行動計量学、2、1−10．
岡本安晴（2008）．能力カテゴリに基づくIRTのベイズ的分析．日本テスト学会第６回大会発表論文抄録集、98−101．
岡本安晴（2009）．離散型および連続型IRTモデルの比較．日本テスト学会第７回発表論文抄録集、142−145．
Y. Okamoto “Item Response Theory with Latent Classes”, Faculty of Integrated Arts and Social Sciences journal, 24, 103-107, 2013. Japan Women’s University. (Click on http://id.nii.ac.jp/1133/00001684/)