Up

潜在クラス項目反応理論

発達段階理論による尺度構成モデル

現代テスト理論の項目反応理論（Item Response Model: IRT）では、正答確率を特性値の関数として

（１）

と表す。ここで、特性値は連続変数である。しかし、Piagetの発達段階理論に見られるように、心理学的概念あるいは特性は心理特性の段階を表すカテゴリ変数の場合もある。このとき、正答確率は心理特性の段階の関数として

（２）

と表される。生澤はモデル（２）を用いた研究を報告しており、学童の興味発達の潜在構造分析（生澤、1965）、知能発達の基本構造（生澤、1972）、潜在クラス分析におけるGreen解の研究（生澤、1974）を上げることができる。

モデル（２）による確率モデルのベイズ分析を岡本（2008、2009、2013、2014）は報告しているが、Metropolis-Hasting within GibbsアルゴリズムをC++で実装したものが用いられている。

本ウェブサイトでは、ベイズ分析を行うスクリプトをStanで開発してみたところ、以下に示すように簡単なスクリプトになることが分かった。

スクリプトを示す前に、モデル（２）について少し説明を加えておく。

まず、モデル（２）で表される正答確率は、段階が進むにつれて増加する条件、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

である。

以下に示すStanスクリプトでは、p_stepはsimplex[NStg+1]型のM次元配列として宣言されている。NStgはカテゴリ数K、Mは項目数である。

段階理論あるいは潜在クラスモデルでは、各回答者の属するクラスはカテゴリ変数であるが、Stanではパラメータとしてカテゴリ変数を設定することはできない。これに対応するため、パラメータの確率の周辺化を行い、潜在クラスを直接パラメータとして扱うことを避ける。次の周辺化を行い、回答者の正答確率を設定する。

以下のスクリプトでは、確率をsimplex[NStg]型のN次元配列thetaで表している。Nは回答者の数である。

である。K＝NStg＝カテゴリ数であることに注意。

なお、カテゴリ変数もパラメータとして扱うことのできるソフトとしてPyMC3があるが、PyMC3によるスクリプトは別のウェブサイトで説明する。

リスト１に上で説明したモデルのサンプルスクリプトを示すが、スクリプトファイルは圧縮ファイルIRTLClassStanFiles.zipとして用意した。また、データをファイルから入力する一般的な使用を想定したサンプルスクリプトは、このウェブサイトの後半にリスト２として用意した。PythonでのStan（PyStan）の利用の説明は＜岡本安晴「いまさら聞けないPythonでデータ分析――多変量解析、ベイズ分析（PyStan、PyMC）――」丸善出版＞で行っている。

リスト１　発達段階理論用尺度構成モデルのデモ用サンプルスクリプト

"""

Yasuharu Okamoto, 2018.11, 2019.02

"""

import pystan

import pickle

import matplotlib.pyplot as plt

import seaborn as sb

import numpy as np

X = ([[1]*5 + [1]*5] * 7) + ([[2]*5 + [1]*5] * 7) + \

([[2]*5 + [2]*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

print('NStg = ', NStg)

a = np.ones(NStg + 1)

stan_data = {'N': N, 'M': M, 'NStg': NStg, 'X': X, 'a': a}

stan_code = """

//

// Binary Responses for NStg Latent Classes

//

data {

int N; // Number of subjects

int M; // Number of items

int NStg; // Number of stages

int X[N, M]; // Responses: 1 -> failure, 2 -> success

vector[NStg + 1] a; // Parameters for Dirichlet distributions

}

parameters {

simplex[NStg+1] p_step[M];

simplex[NStg] theta[N];

}

transformed parameters {

real p[M, NStg];

real p_person[N, M];

vector[2] p_person_cat[N, M];

for (i in 1:M) {

for (stg in 1:NStg) {

p[i][stg] = 0.0;

for (s in 1:stg){

p[i][stg] += p_step[i][s];

}

}

}

for (j in 1:N) {

for (i in 1:M){

p_person[j][i] = 0.0;

for (stg in 1:NStg){

p_person[j][i] = p_person[j][i] + p[i][stg] * theta[j][stg];

}

}

}

for (j in 1:N) {

for (i in 1:M) {

p_person_cat[j][i][2] = p_person[j][i];

p_person_cat[j][i][1] = 1.0 - p_person_cat[j][i][2];

}

}

}

model {

for (i in 1:M){

p_step[i] ~ dirichlet(a);

}

for (j in 1:N) {

for (i in 1:M){

X[j][i] ~ categorical(p_person_cat[j][i]);

}

}

}

"""

sm = pystan.StanModel(model_code = stan_code)

fit = sm.sampling(data = stan_data, n_jobs = 1)

print(fit)

p_samples = fit['p']

p_item_stage = np.empty((M, NStg))

for i in range(M):

for s in range(NStg):

p_item_stage[i][s] = p_samples[:,i,s].mean()

print('\n\nP(correct|item,stage)\n')

print('{0:>7s}{1:>9s}{2:>9s}{3:>9s}'.format(' ', 'Stage-0', 'Stage-1', 'Stage-2'))

for i in range(M):

print('Item-{0} '.format(i), end = '')

for s in range(NStg):

print('{0:>9.5f}'.format(p_item_stage[i][s]), end = '')

print(' ')

theta = fit['theta']

p_stage = np.empty((N, NStg))

for j in range(N):

for s in range(NStg):

p_stage[j][s] = theta[:,j,s].mean()

print('\n\nP(stage|person)\n')

print('{0:>10s}{1:>9s}{2:>9s}{3:>9s}'.format(' ', 'stage-0', 'stage-1', 'stage-2'))

for j in range(N):

print('Person-{0: <3d}'.format(j), end = '')

for s in range(NStg):

print('{0:>9.5f}'.format(p_stage[j][s]), end = '')

print(' ')

リスト１のスクリプトを実行すると、各回答者の反応データが次のように出力される。

X =

[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]

[2, 2, 2, 2, 2, 1, 1, 1, 1, 1]

[2, 2, 2, 2, 2, 1, 1, 1, 1, 1]

[2, 2, 2, 2, 2, 1, 1, 1, 1, 1]

[2, 2, 2, 2, 2, 1, 1, 1, 1, 1]

[2, 2, 2, 2, 2, 1, 1, 1, 1, 1]

[2, 2, 2, 2, 2, 1, 1, 1, 1, 1]

[2, 2, 2, 2, 2, 1, 1, 1, 1, 1]

[2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

[2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

[2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

[2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

[2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

[2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

[2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

１０項目、２１人分のデータであり、各人のデータが各１行に表示される。正答は２で、誤答は１で表されている。

Stanによるサンプリングが終了すると、パラメータのサンプリング値から以下の数値が出力される。

まず、の平均値が、項目ごとに各ステージにおける値として以下のように表示される。

P(correct|item,stage)

Stage-0 Stage-1 Stage-2

Item-0 0.16381 0.76938 0.92920

Item-1 0.16745 0.77199 0.92807

Item-2 0.16629 0.77151 0.92881

Item-3 0.16538 0.77150 0.92966

Item-4 0.16494 0.77240 0.92997

Item-5 0.07067 0.22874 0.83420

Item-6 0.07121 0.22893 0.83373

Item-7 0.07105 0.23268 0.83208

Item-8 0.07175 0.22785 0.83237

Item-9 0.07178 0.23026 0.83233

ステージ０のときは、どの項目も正答確率は0.1前後と低く、ステージ１では、Item-0からItem-4までは正答確率は0.8に近いが、Item-5からItem-9までは約0.2で低い。ステージ２になるとどの項目の正答確率も0.8以上と高い値になる。

次に各回答者がどのステージにいるかを表す確率の平均値がステージごとに以下のように表示される。

P(stage|person)

stage-0 stage-1 stage-2

Person-0 0.75444 0.15627 0.08929

Person-1 0.75426 0.15666 0.08908

Person-2 0.75434 0.15862 0.08704

Person-3 0.75016 0.16200 0.08783

Person-4 0.75008 0.15896 0.09096

Person-5 0.75240 0.15717 0.09042

Person-6 0.75063 0.16133 0.08804

Person-7 0.21065 0.57399 0.21536

Person-8 0.21214 0.57046 0.21740

Person-9 0.21686 0.56770 0.21543

Person-10 0.21738 0.56599 0.21663

Person-11 0.21470 0.57053 0.21477

Person-12 0.21601 0.56789 0.21611

Person-13 0.21403 0.57072 0.21525

Person-14 0.08895 0.16085 0.75020

Person-15 0.08913 0.16293 0.74794

Person-16 0.08851 0.16073 0.75076

Person-17 0.08788 0.16256 0.74956

Person-18 0.08862 0.16061 0.75078

Person-19 0.08778 0.15900 0.75322

Person-20 0.08760 0.16251 0.74989

Person-0からPerson-6まではステージ０である確率が約0.75と高く、Person-7からPerson-13まではステージ１の確率が約0.57と一番高い。Person-14からPerson-20までは、ステージ２の確率が約0.75と一番高い。

入力データをファイルとして用意する場合

リスト１のデモスクリプトでは、データXはプログラム内で生成されている。しかし実データを分析する場合は、データをファイルとして用意してプログラムで読み込めるようにしておくと使い易い。CSV形式で用意されたデータを読み込むスクリプト例を作成した（リスト２）。CSV形式のファイルはExcelなどで作成できる。図１は、上のリスト１におけるデータXをExcelに入力したものである。１行目に変数名を設定して、２行目以降に１行に1人のデータを書いていく。第１欄に回答者のIDとしての文字列を設定し（数値は実行時にエラーとなる）、第２欄目以降に各項目の解答の正誤を記入する。正答は２、誤答は１である。

図１

ファイルの保存のときに、「ファイルの種類」欄で拡張子がCSVであるものを選べば、CSV形式のファイルとして保存される（図２）。

図２

リスト２に示すスクリプト例では、StanスクリプトもファイルIRTLatentClass.stanとして用意され、Pythonスクリプトの実行時に読み込まれている。Pythonスクリプト、Stanスクリプト、サンプルデータファイルは圧縮ファイルIRTLClassStanFilesRD.zipとして用意した。

リスト2　データをファイルから読み込むスクリプト例

"""

Yasuharu Okamoto, 2018.11, 2019.03

"""

import csv

import pystan

import pickle

import matplotlib.pyplot as plt

import seaborn as sb

import numpy as np

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))

#raise Exception()

NStg = int(input('Number of stages = '))

print('NStg = ', NStg)

a = np.ones(NStg + 1)

stan_data = {'N': N, 'M': M, 'NStg': NStg, 'X': X, 'a': a}

sm = pystan.StanModel(file = 'IRTLatentClass.stan')

fit = sm.sampling(data = stan_data, n_jobs = 1)

p_samples = fit['p']

p_item_stage = np.empty((M, NStg))

for i in range(M):

for s in range(NStg):

p_item_stage[i][s] = p_samples[:,i,s].mean()

print('\n\nP(correct|item,stage)\n')

fout.write('\n\nP(correct|item,stage)\n\n')

print('{0:>7s}'.format('Item'), end = '')

fout.write('{0:>7s}'.format('Item'))

for j in range(NStg):

print('{0:>9s}'.format('Stage-{}'.format(j)), end = '')

fout.write('{0:>9s}'.format('Stage-{}'.format(j)))

print()

fout.write('\n')

for i in range(M):

print('{0:>7s}'.format(Items[i]), end = '')

fout.write('{0:>7s}'.format(Items[i]))

for s in range(NStg):

print('{0:>9.5f}'.format(p_item_stage[i][s]), end = '')

fout.write('{0:>9.5f}'.format(p_item_stage[i][s]))

print(' ')

fout.write('\n')

theta = fit['theta']

p_stage = np.empty((N, NStg))

for j in range(N):

for s in range(NStg):

p_stage[j][s] = theta[:,j,s].mean()

print('\n\nP(stage|person)\n')

fout.write('\n\nP(stage|person)\n\n')

print('{0:>10s}'.format('Person'), end = '')

fout.write('{0:>10s}'.format('Person'))

for s in range(NStg):

print('{0:>9s}'.format('Stage-{}'.format(s)), end = '')

fout.write('{0:>9s}'.format('Stage-{}'.format(s)))

print()

fout.write('\n')

for j in range(N):

print('{0:>10s}'.format(CaseID[j]), end = '')

fout.write('{0:>10s}'.format(CaseID[j]))

for s in range(NStg):

print('{0:>9.5f}'.format(p_stage[j][s]), end = '')

fout.write('{0:>9.5f}'.format(p_stage[j][s]))

print(' ')

fout.write('\n')

fout.close()

print('\nOutput file name =', fout_name)

上のスクリプト（リスト２）を実行すると、入力データファイル名と計算結果の出力用ファイル名の設定が求められる（図３）。

図３

図３では、入力データファイル名として図１のExcelで作成されたCSV形式のファイルData.csvが設定され、計算結果の出力用ファイル名としてResults.txtが設定されている。

出力用ファイル名を設定するとデータの読み込みが始まり、回答者数Nと項目数Mの出力後、ステージ数（カテゴリー数）の設定が求められる（図４）。

図４

図４では、ステージ数として３が設定されている。ステージ数の設定後、Stanスクリプトのコンパイルビルドと実行が始まる。Stanスクリプトの実行後、計算結果がリスト１の場合と同様に以下のように出力される。上のリスト１の結果と完全には一致しないのは、ランダムサンプリングに基づく推定だからである。

P(correct|item,stage)

Item Stage-0 Stage-1 Stage-2

X1 0.16341 0.77434 0.93005

X2 0.16621 0.77096 0.92992

X3 0.16507 0.77125 0.92960

X4 0.16388 0.77241 0.92974

X5 0.16708 0.77287 0.93038

X6 0.07277 0.23114 0.83089

X7 0.07045 0.22672 0.83532

X8 0.07036 0.22967 0.83383

X9 0.07011 0.23023 0.83443

X10 0.07131 0.22905 0.83218

P(stage|person)

Person Stage-0 Stage-1 Stage-2

psn1 0.75121 0.16025 0.08854

psn2 0.75147 0.16164 0.08688

psn3 0.74933 0.16278 0.08788

psn4 0.75372 0.15932 0.08696

psn5 0.75555 0.15774 0.08670

psn6 0.75256 0.15832 0.08911

psn7 0.75171 0.15846 0.08983

psn8 0.21270 0.56814 0.21916

psn9 0.21607 0.57165 0.21228

psn10 0.21437 0.57045 0.21518

psn11 0.21450 0.57009 0.21541

psn12 0.21640 0.56648 0.21712

psn13 0.21701 0.56655 0.21644

psn14 0.20983 0.57570 0.21447

psn15 0.08893 0.16141 0.74965

psn16 0.08765 0.16090 0.75146

psn17 0.08799 0.15971 0.75230

psn18 0.08824 0.15822 0.75354

psn19 0.08691 0.16451 0.74858

psn20 0.08964 0.15892 0.75144

psn21 0.08789 0.16576 0.74634

Output file name = Results.txt

>>>

上の計算結果は、出力ファイル（上の場合は、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/)

岡本安晴（2014）．心理学データ分析と測定．勁草書房．Pp. 140-143「離散型項目反応モデル」

Up