項目反応理論による２値項目の分析
Stanスクリプト
各項目が２値（「正答」と「誤答」、「Yes」と「No」など）をとる場合、第
項目の回答を
で表すと
とおける。各項目が３値以上の順序カテゴリ値をとる場合については別のウェブサイトで扱っている。
2値の場合、
とおく。詳しくは、＜岡本安晴「心理学データ分析と測定」、第９章；岡本安晴「計量心理学」、第６章＞などを参照されたい。
上のモデルに基づいてStanスクリプトを以下のように用意した。
＝＝＝＝＝＝＝＝＝＝＝＝ irt_bin.stan ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
data {
int Npsn;
//
Number of persons
int Nitm;
//
Number of items
int Ntot;
// Number of data, Ntot = Nprn * Nitm
int<lower = 1, upper =
Npsn> IDpsn[Ntot]; // Person ID
int<lower = 1, upper =
Nitm> IDitm[Ntot]; // Item ID
int<lower = 1, upper =
2> Res[Ntot];
//
Response-> 1: No,
2: Yes
}
parameters {
real<lower = 0.0> a[Nitm];
real
b[Nitm];
real
F[Npsn];
}
transformed parameters {
vector[2] p[Ntot];
for (i in 1:Ntot){
p[i][2] = inv_logit(1.7 * a[IDitm[i]] * (F[IDpsn[i]] - b[IDitm[i]]));
p[i][1] = 1.0 - p[i][2];
}
}
model {
for (i in 1:Npsn)
F[i]
~ normal(0.0, 1.0);
for (i in 1:Ntot)
Res[i] ~ categorical(p[i]);
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のStanスクリプトを用いるPyhtonスクリプトを以下のように用意した。プログラムのファイルは、サンプルデータファイルも含めて圧縮ファイルPrgIRTbin.zipにまとめた。ファイル名PrgIRTbin.zipをクリックしてダウンロードすることができる。なお、このPythonスクリプトでは能力の推定値は出力されないが、特性（能力）
（上のStanスクリプトではパラメータＦで表されている）の推定値も出力するスクリプトはこのウェブサイトの後半に示した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ main.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import pystan
def calc_mean_med( d
):
d.sort()
med = d[len(d) // 2]
sum = 0
for v in d:
sum
+= v
mean = sum / len(d)
return mean, med
#
# Prepare
the input data file
#
fn_in =
input("Input data file = ")
f_in = open(fn_in,
"r")
#
# Set the
contents of the input data file in the object data_f
#
data_f =
f_in.readlines()
f_in.close()
#
# Prepare
the output file
#
fn_out =
input("Output data file = ")
f_out = open(fn_out,
"w")
f_out.write("Input
data file = " + fn_in + "\n\n")
pos = 0
while True:
if len(data_f[pos]) > 0:
if
data_f[pos][0] == "/":
break
pos += 1
pos += 1
N = int(data_f[pos])
pos += 1
M = int(data_f[pos])
print('N = ', N,
' M = ', M)
X = []
for i in range(N):
X.append([])
pos += 1
temp_strs =
data_f[pos].split()
for j in range(M):
X[i].append(int(temp_strs[j + 1]))
IDpsn = []
IDitm = []
Res = []
for i in range(N):
for j in range(M):
IDpsn.append(i + 1)
IDitm.append(j + 1)
Res.append(X[i][j] + 1)
f_out.write('\nData...\n')
for i in range(N):
f_out.write('{0:>5}: '.format(i+1))
for j in range(M):
f_out.write(' {0}'.format(X[i][j]))
f_out.write('\n')
Data = {'Npsn': N,
'Nitm': M, 'Ntot': N * M, 'IDpsn': IDpsn, 'IDitm': IDitm, 'Res': Res}
fit =
pystan.stan(file = 'irt_bin.stan', data = Data,
seed = 999, pars = ["a", "b"], n_jobs = 1)
print(fit)
a = fit['a']
b = fit['b']
print('len(a) = ',
len(a))
print('len(b) = ',
len(b))
a_smpls = []
b_smpls = []
for i in range(M):
a_smpls.append([])
b_smpls.append([])
for v in a:
for i in range(M):
a_smpls[i].append(v[i])
for v in b:
for i in range(M):
b_smpls[i].append(v[i])
print('\nlen(a) = ',
len(a), ' len(a_smpls) = ',
len(a_smpls),
' len(a_smpls[0] = ', len(a_smpls[0]))
print('\nlen(b) = ',
len(b), ' len(b_smpls) = ',
len(b_smpls),
' len(b_smpls[0] = ',
len(b_smpls[0]))
f_out.write('\n')
for i in range(M):
mean_a, med_a =
calc_mean_med(a_smpls[i])
f_out.write('a[{0}] =
{1:<.5}(mean) {2:<.5}(med)\n'.
format(i+1, mean_a, med_a))
f_out.write('\n')
for i in range(M):
mean_b, med_b =
calc_mean_med(b_smpls[i])
f_out.write('b[{0}] =
{1:<.5}(mean)
{2:<.5}(med)\n'.
format(i+1, mean_b, med_b))
f_out.close()
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のPythonスクリプトを実行すると、入力データファイル名と出力データファイル名の設定が求められる（図１）。
図１
出力データ用ファイル名は適当なテキストファイル名を設定すればよいが、入力データファイルは図２に示す形式で用意する。
・
・
・
図２
先頭がスラッシュ / で始まる行の次の行にデータ数（回答者数）を書く。その次の行には項目数を書く。その次の行から、１行に1人の回答データを書いていく。まず、先頭に回答者ＩＤを書く。通し番号を書いておくとよいが、文字列として読み込まれるので通し番号以外でもよい。回答者ＩＤに続いて回答カテゴリ値（１または０）を空白文字で区切って並べていく。
入力データファイル名の入力の次に、出力用ファイル名の設定が求められる。出力用ファイル名には適当なテキストファイル名を設定する。出力用ファイル名を入力すると、計算が始まる。回答者数Ｎ，項目数Ｍの値が出力された後、Stanによる計算が始まる（図３）。
図３
回答者数300人、項目数10個の図２のデータの場合、図４では１つのchain当たりの計算時間が約１分になっている。Chainが１つずつ順番に実行される場合には、全体の計算時間はそれらのchainの計算時間の和になる。
図４
計算が終わると図５のようにStanによる計算の統計量が表示される。
図５
プログラムの終了後、出力ファイルを開くと図６、図７のようになっている。
図６
まず、読み込んだデータが通し番号付きで書き出されている。これにより、データが正しく読み込まれたかどうかの確認ができる。読み込みデータの書き出しに続いて、各パラメータの推定値が平均値と中央値で示されている（図７）。
図７
必要に応じてプログラムを書き換えれば、種々の分析が可能である。
特性（能力）
の推定値も出力するスクリプト
特性値
（StanスクリプトではパラメータＦで表されている）の推定値も図Ｔ１のように出力されるようにしたPythonスクリプトも作成した。
・
・
・
図Ｔ１
以下のPythonスクリプトでは、特性（能力）パラメータＦをリストthetaに読み込んで処理をしている。スクリプトとサンプルデータファイルは圧縮ファイルPrgIRTbin_theta.zipにまとめたので、ファイル名PrgIRTbin_theta.zipをクリックしてダウンロードできる。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ main_theta.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import pystan
def calc_mean_med( d ):
d.sort()
med = d[len(d) // 2]
sum = 0
for v in d:
sum
+= v
mean = sum / len(d)
return mean, med
#
# Prepare
the input data file
#
fn_in = input("Input data file =
")
f_in = open(fn_in, "r")
#
# Set the
contents of the input data file in the object data_f
#
data_f = f_in.readlines()
f_in.close()
#
# Prepare
the output file
#
fn_out = input("Output data file
= ")
f_out = open(fn_out, "w")
f_out.write("Input data file =
" + fn_in + "\n\n")
pos = 0
while True:
if len(data_f[pos]) > 0:
if
data_f[pos][0] == "/":
break
pos += 1
pos += 1
N = int(data_f[pos])
pos += 1
M = int(data_f[pos])
print('N = ', N, ' M = ', M)
X = []
for i in range(N):
X.append([])
pos += 1
temp_strs =
data_f[pos].split()
for j in range(M):
X[i].append(int(temp_strs[j + 1]))
IDpsn = []
IDitm = []
Res = []
for i in range(N):
for j in range(M):
IDpsn.append(i + 1)
IDitm.append(j +
1)
Res.append(X[i][j] + 1)
f_out.write('\nData...\n')
for i in range(N):
f_out.write('{0:>5}:
'.format(i+1))
for j in range(M):
f_out.write(' {0}'.format(X[i][j]))
f_out.write('\n')
Data = {'Npsn': N, 'Nitm': M, 'Ntot':
N * M, 'IDpsn': IDpsn, 'IDitm': IDitm, 'Res': Res}
fit = pystan.stan(file =
'irt_bin.stan', data = Data,
seed = 999, pars = ["a", "b", "F"], n_jobs
= 1)
print(fit)
a = fit['a']
b = fit['b']
theta = fit['F']
print('len(a) = ', len(a))
print('len(b) = ', len(b))
print('len(theta)', len(theta))
print('len(theta[0]) = ',
len(theta[0]))
a_smpls = []
b_smpls = []
for i in range(M):
a_smpls.append([])
b_smpls.append([])
for v in a:
for i in range(M):
a_smpls[i].append(v[i])
for v in b:
for i in range(M):
b_smpls[i].append(v[i])
f_out.write('\n')
for i in range(M):
mean_a, med_a =
calc_mean_med(a_smpls[i])
f_out.write('a[{0}] =
{1:<.5}(mean)
{2:<.5}(med)\n'.
format(i+1, mean_a, med_a))
f_out.write('\n')
for i in range(M):
mean_b, med_b =
calc_mean_med(b_smpls[i])
f_out.write('b[{0}] =
{1:<.5}(mean)
{2:<.5}(med)\n'.
format(i+1, mean_b, med_b))
f_out.write('\ntheta...')
theta_smpls = []
for i in range(N):
theta_smpls.append([])
for v in theta:
for i in range(N):
theta_smpls[i].append(v[i])
for i in range(N):
mean_theta, med_theta =
calc_mean_med(theta_smpls[i])
f_out.write('\nPerson-{0:<7}
{1:>10.5}(mean)
{2:>10.5}(med)'.format(i, mean_theta, med_theta))
f_out.close()
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
プログラムの実行終了後、出力ファイルを開くと図Ｔ２のようになっている。図Ｔ１に示されている能力の推定値（パラメータＦの値として示されている）がthetaの平均値と中央値で出力されている。
図Ｔ２