項目反応理論のパラメータ推定のシミュレーション
事前分布の選択
項目反応理論のパラメータ推定において、バイアスを避けるために多くの研究者は無情報事前分布（フラットな事前分布）を用いると説明されている（Johnson & Sinharay, 2016）。このことを確かめるために行ったシミュレーションでは、以下に示すように、一様分布でバイアスが認められた。また、標準正規分布を事前分布とした場合もバイアスが認められたが、標準偏差を適切に調整した正規分布を事前分布とした場合はバイアスは認められなかった。実施したシミュレーションを、データの生成から順に説明する。使用したスクリプトは、付録に示した。ファイルは、irtsimufiles.zipにまとめた。
データの生成
データ生成のモデルとして１パラメータのIRTモデルを用いた。すなわち、特性値
の解答者
の項目
に対する回答
を次式の確率で与えた。
モデル(1)は、
のときRaschモデルと呼ばれている（de Ayala,
2009）。データ生成においても、とした。
項目数n_itemsは、10項目とした。項目数が少ない場合、例えばn_items=1の場合は、回答パターンは
との2つしかなく、特性値の推定値も２つになる。項目数を増やすと、事前分布の影響が弱くなる。事前分布が推定値の分布に及ぼす影響を検討するためには項目数は10個ぐらいが適当と考え、項目数n_items=10とした。
の値は、最小値
、最大値
として等間隔に設定した。特性値は標準正規分布に従うものとして生成した。
である。人数n_personsは、n_persons=1000人とした。
上の条件でデータを生成するスクリプトとしてリスト１に示すものを作成した。
リスト１のスクリプトを実行すると、図１に示すExcelデータファイルが作成される。
図１ シミュレーションで作成したデータ
１列目に回答者のラベル、２列目に式（２）により生成された特性値（ここでは、特性値の真値と呼ぶ）、３列目以降に回答（１または０）が設定されている。
作成されたデータは、乱数を用いているので、スクリプトの実行ごとに変わる。
作成されたデータ（図１）に対して、１パラメータモデルを適用して、パラメータ値を求める。
分析モデルとベイズ分析
１パラメータ項目反応モデル（１）により、図１のデータを分析する。式（１）は
と
との２種類のパラメータがあるが、は尺度の単位に関わるものであり項目にわたって共通である。項目間の違いはによって表されている。式（１）は、本質的に１パラメータであると考えられる。
モデル（１）の原点と単位を設定するために、特性値の事前分布を標準正規分布とする。
以上の設定の下でベイズ分析（Gelman et al., 2014; 岡本, 2019）を行うStanスクリプトをリスト２のように用意した。
リスト２のStanスクリプトを用いて、図１のデータのベイズ分析を行うPythonスクリプトをリスト３のように用意した。リスト３は、CmdStanPyのインストールされた環境で実行するが、CmdStanPyの説明のウェブサイトも用意している。
リスト３のスクリプトを実行すると、図１のデータに対してMCMCサンプリングによるベイズ分析が行われる。MCMCサンプリング後、トレース図が表示される（図２）。
図２
サンプリングに問題はなさそうである。
得られたMCMCサンプルの分析を進める。
MCMCサンプルの分析
リスト３のスクリプトによる図１のデータのMCMCサンプルの分析を、リスト４のスクリプトで行った。
特性値の推定値を事後分布の中央値とした（Gelman et al., 2021）。
特性値の推定値に「individual theta est.」、特性値の真値（図１の第２列（theta））の値に「org.」、標準正規分布に「std normal distribution」のラベルの付けられたグラフを図３の左側のグラフに示した。右側のグラフは、特性値の真値と特性値の推定値の散布図である。
図３
パラメータの事後分布を図４の左のグラフに、パラメータの事後分布を図４の右側のグラフに示す。
図４
パラメータの推定値として事後分布の中央値を採用した（Gelman et al., 2021）。
得られたパラメータの推定値を項目のパラメータ値として、個人の回答パターンに対して個人の特性値を求めた場合の特性値について検討する。
推定された個人特性の検討
上のベイズ分析で得られたパラメータの推定値を項目の特性を表すパラメータ値としたときの個人の回答パターンに対する個人特性値の推定値について調べる。
使用するモデルは１パラメータモデル（１）式
である。項目パラメータの値に原点と単位の情報が含まれているので、個人特性値の事前分布に一様分布
を設定することができる。一様分布のサポートは、十分広く取っておけばよい。
個人特性値は標準正規分布に従っているという情報を事前分布とすることも考えられる。
事前情報として標準正規分布を設定するモデルをリスト５に、一様分布を設定するモデルをリスト６のようにスクリプトを作成した。
個人回答パターンのデータを、図１のデータから１行ずつとり、３列目からの値を回答パターンのデータとして個人ごとの特性値のベイズ分析をリスト５およびリスト６のスクリプトにより行うスクリプトをリスト７のように作成した。
リスト７によるシミュレーションの結果を表示するために、リスト８のスクリプトを作成した。
リスト８のスクリプトを実行すると、図５のグラフが表示される。
上のグラフの横軸
は、図１のデータにおける第２列（theta）の値で、真値として扱っている。事前分布として標準正規分布を設定した推定値は「normal
prior」で、一様分布を設定したものは「uniform prior」でラベル付けしている。
図５
図５の左下のグラフには、標準正規分布を事前分布とする特性値の分布（青）、一様分布を事前分布とする特性値の分布（橙）とラベル「org」の付けられている真値の分布（緑）が描かれている。一様分布を事前分布とする推定値の分布は真値の分布より広がっており、標準正規分布を事前分布とする推定値は真値の分布より収縮している。一様分布の場合、推定値が広がっているのは、回答パターンが全項目が１の場合、あるいは全項目が０の場合は尤度関数が発散するためと思われる。事前分布として正規分布を設定した場合は、尤度関数の発散は正規分布の裾により抑えられる。
標準正規分布を事前分布としたときの推定値の標準偏差SD(norm priot)と真値の標準偏差SD(org)の比が
であるので、事前分布として標準偏差が
の正規分布を用いて個人別の特性値推定を行ってみる。
標準偏差がの正規分布を用いた場合
事前分布の標準偏差をとした正規分布を用いてみる。
Stanスクリプトはリスト９に示すものである。リスト９のスクリプトを用いて、個人単位の特性値推定のスクリプトを作成したものをリスト10に示す。
リスト10のスクリプトによるシミュレーション結果を、リスト11のスクリプトにより表示すると図６のようになる。
図６
個人単位の特性値の推定値の分布は、真値の分布によく重なっているように見える。
de Ayala, R. J. (2009) The theory and practice of item response theory. The Guilford Press.
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, 3rd. Ed. CRC Press. (free download)
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and Other Stories. Cambridge University Press.
Johnson, M. S., & Sinharay, S. (2016) Bayesian estimation. In van der Linden, W. J.(ed.) Handbook of item response theory, Vol. two: Statistical tools. CRC Press.
岡本安晴（2014） 心理学データ分析と測定．勁草書房（電子書籍版あり）
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版
ファイルは、irtsimufiles.zipにまとめた。
import numpy
as np
import scipy.stats
as ss
import pandas as pd
n_items = 10
a_lst = np.array([1]*n_items)
b_lst = np.linspace(-2,
2, n_items)
print(a_lst)
print(b_lst)
ab_lst = {}
for i, ab in
enumerate(zip(a_lst, b_lst)):
ab_lst['a'+f'{i}'+'b'+f'{i}'] = (ab)
for v in ab_lst:
print(v, ab_lst[v])
resp = {'person': [], 'theta':[]}
for v in ab_lst:
resp[v] = []
rng = np.random.default_rng()
n_persons = 1000
i_person = 0
for i in
range(n_persons):
t = ss.norm.rvs()
#tv
= []
resp['person'].append(f'p{i+1}')
resp['theta'].append(t)
for
v in ab_lst:
a = ab_lst[v][0]
b = ab_lst[v][1]
if ss.uniform.rvs()
< ss.norm.cdf([a * (t - b)]):
resp[v].append(1)
else:
resp[v].append(0)
df_pd = pd.DataFrame(resp)
print(df_pd)
df_pd.to_excel('sample_data.xlsx', index=False)
リスト２ 生成したデータのベイズ分析で用いるStanスクリプト（irt.stan）
functions {
real ck_v(real
v) {
real
tv;
tv =
v;
if
(tv < 1.0e-9) {
tv = 1.0e-9;
}
if
(tv > 1.0-1.0e-9) {
tv = 1.0 - 1.0e-9;
}
return tv;
}
}
data {
int n_persons;
int n_items;
array[n_persons,
n_items] int X;
}
parameters {
array[n_persons]
real t;
real<lower=0,
upper=10> a;
array[n_items]
real<lower=-10, upper=10> b;
}
model {
for (i
in 1:n_persons){
t[i] ~ normal(0, 1);
}
a ~ uniform(0, 10);
for (j in 1:n_items) {
b[j]
~ uniform(-10, 10);
}
for (i
in 1:n_persons) {
for
(j in 1:n_items) {
X[i][j] ~ bernoulli(ck_v(normal_cdf(a * (t[i] - b[j]) | 0, 1)));
}
}
}
リスト３ 生成したデータのベイズ分析のPythonスクリプト（anal_data.py）
from cmdstanpy import CmdStanModel
import numpy
as np
import scipy.stats
as ss
import pandas as pd
import matplotlib.pyplot
as plt
import arviz
as az
import seaborn as sb
import pickle
dframe = pd.read_excel('sample_data.xlsx')
print(dframe)
var_nm = dframe.keys()
print(var_nm)
var_nm = var_nm[1:]
print(var_nm)
raw_data = dframe.values
irt_data = raw_data[:,2:]
print(np.mean(irt_data, axis=0))
sm = CmdStanModel(stan_file = 'irt.stan')
Data = {'n_persons':len(irt_data), 'n_items':len(irt_data[0]), 'X':irt_data}
fit = sm.sample(data=Data)
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
# Arviz(az) InferenceData
# トレースプロット
az.plot_trace(inf_data)
plt.tight_layout()
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame
with open('mcmcsamples.pkl'
, 'wb') as f:
pickle.dump(fit, f)
リスト４ ベイズ分析の結果のチェック（check_mcmcsample.py）
import pickle
import pandas as pd
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
import seaborn as sb
dframe = pd.read_excel('sample_data.xlsx')
print(dframe)
n_persons = len(dframe['person'])
theta_org = dframe['theta'].values
print('theta_org
=\n', theta_org[:5], theta_org[-5:])
var_nm = dframe.keys()
var_nm = var_nm[1:]
n_items = len(var_nm) - 1
raw_data = dframe.values
irt_data = raw_data[:,2:]
with open('mcmcsamples.pkl'
, 'rb') as f:
fit
= pickle.load(f)
theta = []
for i in
range(n_persons):
theta.append(fit[f't[{i+1}]'])
t_means = np.median(theta,
axis=1)
print(len(t_means))
tcoord = np.linspace(-3,
3, 1000)
ycoord = ss.Normal(mu=0,
sigma=1).pdf(tcoord)
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(tcoord, ycoord, label='std normal distribution')
sb.kdeplot(t_means,
label='individual theta est.')
sb.kdeplot(theta_org,
label='orig.')
plt.xlabel(r'$\theta$', fontsize=16)
plt.legend()
plt.subplot(1,2,2)
plt.plot(theta_org, t_means, 'o', markerfacecolor='none')
plt.plot([-3, 3],[-3, 3], c='r')
plt.xlabel(r'$\theta$', fontsize=14)
plt.ylabel(r'$\widehat\theta$',
fontsize=16)
plt.show()
a_smpl = fit['a']
a_est = np.median(a_smpl)
print('a_est
=', a_est)
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
sb.kdeplot(a_smpl,
label=f'Med.={a_est:.3f}')
plt.legend(fontsize=12)
plt.title('posterior distribution of a', fontsize=16)
plt.subplot(1,2,2)
b_smpl = []
for j in range(n_items):
b_smpl.append(fit[f'b[{j+1}]'])
b_est = np.median(b_smpl, axis=1)
print('b_est
=', b_est)
for j in range(n_items):
sb.kdeplot(b_smpl[j], label=f'b[{j+1}]={b_est[j]:.3f}')
plt.title('Posterior distributions of b[j] and
Medians', fontsize=14)
plt.legend(loc='upper left', bbox_to_anchor=(1,
1))
plt.tight_layout()
plt.show()
with open('ab_est.pkl',
'wb') as f:
pickle.dump((a_est, b_est), f)
リスト５ 個人データの標準正規分析を事前分布とするベイズ分析のためのStanスクリプト（ind_normal.stan）
functions {
real ck_v(real
v) {
real
tv;
tv =
v;
if
(tv < 1.0e-9) {
tv = 1.0e-9;
}
if
(tv > 1.0-1.0e-9) {
tv = 1.0 - 1.0e-9;
}
return tv;
}
}
data {
int n_items;
array[n_items]
int X;
real a;
array[n_items]
real b;
}
parameters {
real t;
}
model {
t ~ normal(0, 1);
for (j in 1:n_items) {
X[j]
~ bernoulli(ck_v(normal_cdf(a * (t - b[j]) | 0, 1)));
}
}
リスト６ 個人データの一様分布を事前分布とするベイズ分析のためのStanスクリプト（ind_uniform.stan）
functions {
real ck_v(real
v) {
real
tv;
tv =
v;
if
(tv < 1.0e-9) {
tv = 1.0e-9;
}
if
(tv > 1.0-1.0e-9) {
tv = 1.0 - 1.0e-9;
}
return tv;
}
}
data {
int n_items;
array[n_items]
int X;
real a;
array[n_items]
real b;
}
parameters {
real<lower=-10,
upper=10> t;
}
model {
t ~ uniform(-10, 10);
for (j in 1:n_items) {
X[j]
~ bernoulli(ck_v(normal_cdf(a * (t - b[j]) | 0, 1)));
}
}
リスト７ 個人単位にベイズ分析を行うPythonスクリプト（anal_individual.py）
from cmdstanpy import CmdStanModel
import numpy
as np
import scipy.stats
as ss
import pandas as pd
import matplotlib.pyplot
as plt
import arviz
as az
import seaborn as sb
import pickle
import time
dframe = pd.read_excel('sample_data.xlsx')
print(dframe)
var_nm = dframe.keys()
print(var_nm)
var_nm = var_nm[1:]
print(var_nm)
raw_data = dframe.values
print(raw_data)
irt_data = raw_data[:,2:]
print(irt_data)
print(np.mean(irt_data, axis=0))
with open('ab_est.pkl',
'rb') as f:
(a_est, b_est) = pickle.load(f)
print('a_est
=', a_est)
print('len(b_est) =', len(b_est))
theta_normal = []
theta_uniform = []
n_persons = len(irt_data)
for i in
range(n_persons):
if i % 10 == 0:
print(f'\n{i}/{n_persons}\n')
ck_t = time.time()
while time.time() < ck_t
+ 3:
pass
X_ind = irt_data[i]
sm = CmdStanModel(stan_file = 'ind_normal.stan')
Data
= {'n_items':len(b_est),
'X':X_ind, 'a':a_est, 'b':b_est}
fit
= sm.sample(data=Data, show_progress=False,
max_treedepth=20, adapt_delta=0.99)
fit
= fit.draws_pd()
theta_normal.append(np.median(fit['t']))
sm1
= CmdStanModel(stan_file =
'ind_uniform.stan')
fit1
= sm1.sample(data=Data, show_progress=False,
max_treedepth=20, adapt_delta=0.99)
fit1
= fit1.draws_pd()
theta_uniform.append(np.median(fit1['t']))
sb.kdeplot(theta_normal,
label='normal', alpha=0.5)
sb.kdeplot(theta_uniform,
label='uniform', alpha=0.5)
plt.legend()
plt.show()
with open('theta_esta.pkl',
'wb') as f:
pickle.dump((theta_normal, theta_uniform), f)
リスト８ 個人単位のベイズ分析の事前分布による比較（comp_distri_est.py）
import numpy
as np
import scipy.stats
as ss
import pandas as pd
import matplotlib.pyplot
as plt
import seaborn as sb
import pickle
dframe = pd.read_excel('sample_data.xlsx')
print(dframe)
theta_org = dframe['theta'].values
with open('mcmcsamples.pkl'
, 'rb') as f:
fit_mcmc = pickle.load(f)
n_persons = len(dframe['person'])
raw_data = dframe.values
irt_data = raw_data[:,2:]
n_persons = len(irt_data)
with open('theta_esta.pkl',
'rb') as f:
(theta_normal,
theta_uniform) = pickle.load(f)
theta_org = np.array(theta_org)
theta_normal = np.array(theta_normal)
theta_uniform = np.array(theta_uniform)
sse_normal = np.mean((theta_normal - theta_org)**2)
sse_uniform = np.mean((theta_uniform - theta_org)**2)
plt.figure(figsize=(10,6))
plt.subplot(2,4, (2,3))
plt.plot(theta_org,theta_normal,
'o',
markerfacecolor='none', label='normal prior')
plt.plot(theta_org, theta_uniform, 'o',
markerfacecolor='none', label='uniform prior')
plt.plot([-6, 6], [-6, 6], color='m',
label='y=x')
plt.xticks(np.arange(-6,
6.1, 2))
plt.xlabel(r'$\theta$', fontsize=
16)
plt.ylabel(r'$\widehat{\theta}$',
fontsize= 18)
plt.title(f'sse_normal={sse_normal:.3f}, sse_uniform={sse_uniform:.3f}')
plt.legend()
plt.subplot(2,4,(5,6))
plt.title('Distributions', fontsize=16)
sb.kdeplot(theta_normal,
label='normal prior')
sb.kdeplot(theta_uniform,
label='uniform prior')
sb.kdeplot(theta_org,
label='org.')
plt.xlabel(r'$\theta, \widehat{\theta}$',
fontsize=18)
plt.legend()
plt.tight_layout()
theta_org_s = theta_org[np.argsort(theta_org)]
theta_normal_s = theta_normal[np.argsort(theta_normal)]
theta_uniform_s = theta_uniform[np.argsort(theta_uniform)]
n_data = len(theta_org)
cnt_ary = [1]*n_data
cum_org = np.cumsum(cnt_ary) / n_data
plt.subplot(2,4,(7,8))
plt.plot(theta_normal_s,
cum_org, label='normal prior')
plt.plot(theta_uniform_s,
cum_org, label='uniform prior')
plt.plot(theta_org_s,
cum_org, label='org.')
plt.title('Cumulative distributions', fontsize=16)
plt.xlabel(r'$\theta, \widehat{\theta}$',
fontsize=18)
plt.legend(fontsize=14)
plt.tight_layout()
plt.show()
sd_org = np.std(theta_org)
sd_normal = np.std(theta_normal)
print('normal/pop =', sd_normal/sd_org)
リスト９ SD=1/0.88による個人単位のベイズ分析のためのStanスクリプト（ind_normal_88.stan）
functions {
real ck_v(real
v) {
real
tv;
tv =
v;
if
(tv < 1.0e-9) {
tv = 1.0e-9;
}
if
(tv > 1.0-1.0e-9) {
tv = 1.0 - 1.0e-9;
}
return tv;
}
}
data {
int n_items;
array[n_items]
int X;
real a;
array[n_items]
real b;
}
parameters {
real t;
}
model {
t ~ normal(0, 1/0.88);
for (j in 1:n_items) {
X[j]
~ bernoulli(ck_v(normal_cdf(a * (t - b[j]) | 0, 1)));
}
}
リスト10 SD=1/0.88による個人単位のベイズ分析（anal_individual_n88.py）
from cmdstanpy import CmdStanModel
import numpy
as np
import scipy.stats
as ss
import pandas as pd
import matplotlib.pyplot
as plt
import arviz
as az
import seaborn as sb
import pickle
import time
dframe = pd.read_excel('sample_data.xlsx')
print(dframe)
var_nm = dframe.keys()
print(var_nm)
var_nm = var_nm[1:]
print(var_nm)
raw_data = dframe.values
print(raw_data)
irt_data = raw_data[:,2:]
print(irt_data)
print(np.mean(irt_data, axis=0))
with open('ab_est.pkl',
'rb') as f:
(a_est, b_est) = pickle.load(f)
print('a_est
=', a_est)
print('len(b_est) =', len(b_est))
theta_normal = []
theta_uniform = []
n_persons = len(irt_data)
for i in
range(n_persons):
if i % 10 == 0:
print(f'\n{i}/{n_persons}\n')
ck_t = time.time()
while
time.time() < ck_t + 3:
pass
X_ind = irt_data[i]
sm = CmdStanModel(stan_file = 'ind_normal_88.stan')
Data
= {'n_items':len(b_est),
'X':X_ind, 'a':a_est, 'b':b_est}
fit
= sm.sample(data=Data, show_progress=False,
max_treedepth=20, adapt_delta=0.99)
fit
= fit.draws_pd()
theta_normal.append(np.median(fit['t']))
"""
sm1
= CmdStanModel(stan_file =
'ind_uniform.stan')
fit1
= sm1.sample(data=Data, show_progress=False,
max_treedepth=20, adapt_delta=0.99)
fit1
= fit1.draws_pd()
theta_uniform.append(np.median(fit1['t']))
"""
sb.kdeplot(theta_normal,
label='normal_88', alpha=0.5)
#sb.kdeplot(theta_uniform,
label='uniform', alpha=0.5)
plt.legend()
plt.show()
with open('theta_est_n88.pkl', 'wb') as f:
pickle.dump(theta_normal, f)
リスト11 SD=1/0.88による個人単位のベイズ分析との比較（comp_distri_est_n88.py）
import numpy
as np
import scipy.stats
as ss
import pandas as pd
import matplotlib.pyplot
as plt
import seaborn as sb
import pickle
dframe = pd.read_excel('sample_data.xlsx')
print(dframe)
theta_org = dframe['theta']
with open('theta_est_n88.pkl', 'rb') as f:
theta_normal_n88 = pickle.load(f)
plt.figure(figsize=(12,5))
plt.subplot(1,2,2)
sb.kdeplot(theta_normal_n88, label='prior sd=1/0.88')
sb.kdeplot(theta_org,
label='org.')
plt.xlabel(r'$\theta, \widehat{\theta}$',
fontsize=18)
plt.title('Distributions', fontsize=16)
plt.legend(fontsize=14)
plt.tight_layout()
theta_org = np.array(theta_org)
theta_normal_n88 = np.array(theta_normal_n88)
theta_org_s = theta_org[np.argsort(theta_org)]
theta_normal_n88_s = theta_normal_n88[np.argsort(theta_normal_n88)]
plt.subplot(1,2,1)
n_data = len(theta_org)
cnt_ary = [1]*n_data
cum_prop = np.cumsum(cnt_ary) / n_data
plt.plot(theta_normal_n88_s, cum_prop, label='prior sd=1/0.88')
plt.plot(theta_org_s,
cum_prop, label='org.')
plt.xlabel(r'$\theta, \widehat{\theta}$',
fontsize=18)
plt.title('Cumulative distributions', fontsize=16)
plt.legend(fontsize=14)
plt.tight_layout()
plt.show()