Up

項目反応理論による多値項目の分析

Stanスクリプト

各項目が３値以上の順序カテゴリ値をとる場合について考える。各項目の値が2値の場合については別のウェブサイトで扱っている。多値の場合のモデルとしてgraded response model (Samejima, 1969)と呼ばれているものがよく知られているが、以下のようなものである。

カテゴリの数をＫとし、カテゴリ値を１からＫまでの整数で表す。第項目の選択カテゴリをで表すと

である。例えば、Ｋ＝６のときは、「１：大反対」、「２：反対」、「３：どちらかといえば反対」、「４：どちらかといえば賛成」、「５：賛成」、「６：大賛成」という具合である。

カテゴリの選択確率を次のようにおく。

ここで、

である。項目特性を表すパラメータがとである。カテゴリ値の境界に対応するパラメータがである。モデル式では、の形で表れているので、パラメタ値が不定になる。パラメタ値のidentificationが可能になるように項目１のパラメタを

とおく。このとき、カテゴリ値の境界のパラメータは項目１を基準にして設定され、難易度パラメータは項目１を基準とする項目の相対的な難易度を表す。他の値を基準にするときは、任意の値に対して

が成り立つので、基準にする値に合わせて各に同じ数を足し、各から同じ数を引けばよい。

以上の設定の下でStanスクリプトを以下のように用意した。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ irt_pol.sta ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

data {

int <lower = 3> K; // Number of categories

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 = K> Res[Ntot]; // Response-> an integer from 1 to K

}

parameters {

real<lower = 0.0> a[Nitm];

real b[Nitm - 1];

ordered[K - 1] c;

real F[Npsn];

}

transformed parameters {

vector[K] p[Ntot];

vector[K] shft[Nitm];

for (k in 1:(K-1)){

shft[1][k] = c[k];

for (j in 2:(Nitm)){

shft[j][k] = shft[1][k] + b[j - 1];

}

}

for (i in 1:Ntot){

p[i][K] = inv_logit(1.7 * a[IDitm[i]] * (F[IDpsn[i]] - shft[IDitm[i]][K - 1]));

p[i][1] = 1.0 - inv_logit(1.7 * a[IDitm[i]] * (F[IDpsn[i]] - shft[IDitm[i]][1]));

for (k in 2:(K - 1)){

p[i][k] = inv_logit(1.7 * a[IDitm[i]] * (F[IDpsn[i]] - shft[IDitm[i]][k-1]))

-inv_logit(1.7 * a[IDitm[i]] * (F[IDpsn[i]] - shft[IDitm[i]][k]));

}

}

}

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スクリプトを以下のように用意した。プログラムのファイルは、サンプルデータファイルも含めて圧縮ファイルPrgIRTpol.zipにまとめた。ファイル名PrgIRTpol.zipをクリックしてダウンロードすることができる。なお、このPythonスクリプトでは能力の推定値は出力されないが、特性（能力）（上のStanスクリプトではパラメータＦで表されている）の推定値も出力するスクリプトはこのウェブサイトの後半に示した。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

import pystan

from pystan import StanModel

import pickle

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

pos += 1

K = int(data_f[pos])

print('N = ', N, ' M = ', M, ' K = ', K)

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

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 = {'K': K, 'Npsn': N, 'Nitm': M, 'Ntot': N * M, 'IDpsn': IDpsn, 'IDitm': IDitm, 'Res': Res}

fit = pystan.stan(file = 'irt_pol.stan', data = Data, seed = 999,

pars = ['a', 'b', 'c'], n_jobs = 1)

print(fit)

a = fit['a']

b = fit['b']

c = fit['c']

a_smpls = []

b_smpls = []

c_smpls = []

for i in range(M):

a_smpls.append([])

for i in range(M-1):

b_smpls.append([])

for k in range(K-1):

c_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-1):

b_smpls[i].append(v[i])

for v in c:

for i in range(K-1):

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

print('\nlen(c) = ', len(c), ' len(c_smpls) = ', len(c_smpls),

' len(c_smpls[0] = ', len(c_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')

f_out.write('b[1] = 0.0 (indentification condition)\n')

for i in range(M-1):

mean_b, med_b = calc_mean_med(b_smpls[i])

f_out.write('b[{0}] = {1:<.5}(mean) {2:<.5}(med)\n'.

format(i+2, mean_b, med_b))

f_out.write('\n')

for i in range(K-1):

mean_c, med_c = calc_mean_med(c_smpls[i])

f_out.write('c[{0}] = {1:<.5}(mean) {2:<.5}(med)\n'.

format(i+1, mean_c, med_c))

f_out.close()

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

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

図１

出力データ用ファイル名は適当なテキストファイル名設定すればよいが、入力データファイルは図２に示す形式で用意する。

・

・

・

図２

先頭がスラッシュ / で始まる行の次の行にデータ数（回答者数）を書く。その次の行には項目数を書き、次の行にカテゴリ数を書く。カテゴリ数を書いた次の行から、１行に1人の回答データを書いていく。まず、先頭に回答者ＩＤを書く。通し番号を書いておくとよいが、文字列として読み込まれるので通し番号以外でもよい。回答者ＩＤに続いて回答カテゴリ値を空白文字で区切って並べていく。

入力データファイル名の入力の次に、出力用ファイル名の設定が求められる。出力用ファイル名には適当なテキストファイル名を設定する。出力用ファイル名を入力すると、計算が始まる。回答者数Ｎ，項目数Ｍ，カテゴリ数Ｋの値が出力された後、Stanによる計算が始まる（図３）。

図３

回答者数200人、項目数10個、カテゴリ数６個の図２のデータの場合、図４では１つのchain当たりの計算時間が５分になっている。Chainが１つずつ順番に実行される場合には、全体の計算時間はそれらのchainの計算時間の和になる。

図４

計算が終わると図５のようにStanによる計算の統計量が表示される。

図５

プログラムの終了後、出力ファイルを開くと図６、図７のようになっている。

図６

まず、読み込んだデータが通し番号付きで書き出されている。これにより、データが正しく読み込まれたかどうかの確認ができる。読み込みデータの書き出しに続いて、各パラメータの推定値が平均値と中央値で示されている（図７）。

図７

必要に応じてプログラムを書き換えれば、種々の分析が可能である。

図７の出力では、パラメータ値の同定のための条件の基準で算出された値が示されているが、およびの値を式（１）により任意の値によりおよびにシフトすることができる。いま、カテゴリ境界の中央の値を０とする値にシフトするためにとおいて計算すると以下のようになる。

特性（能力）の推定値も出力するスクリプト

特性値（StanスクリプトではパラメータＦで表されている）の推定値も図Ｔ１のように出力されるようにしたPythonスクリプトも作成した。

・

・

・

図Ｔ１

以下のPythonスクリプトでは、特性（能力）パラメータＦをリストthetaに読み込んで処理をしている。スクリプトとサンプルデータファイルは圧縮ファイルPrgIRTpol_theta.zipにまとめたので、ファイル名PrgIRTpol_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])

pos += 1

K = int(data_f[pos])

print('N = ', N, ' M = ', M, ' K = ', K)

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

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 = {'K': K, 'Npsn': N, 'Nitm': M, 'Ntot': N * M, 'IDpsn': IDpsn, 'IDitm': IDitm, 'Res': Res}

fit = pystan.stan(file = 'irt_pol.stan', data = Data, seed = 999,

pars = ['a', 'b', 'c', 'F'], n_jobs = 1)

print(fit)

a = fit['a']

b = fit['b']

c = fit['c']

theta = fit['F']

a_smpls = []

b_smpls = []

c_smpls = []

for i in range(M):

a_smpls.append([])

for i in range(M-1):

b_smpls.append([])

for k in range(K-1):

c_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-1):

b_smpls[i].append(v[i])

for v in c:

for i in range(K-1):

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

f_out.write('b[1] = 0.0 (indentification condition)\n')

for i in range(M-1):

mean_b, med_b = calc_mean_med(b_smpls[i])

f_out.write('b[{0}] = {1:<.5}(mean) {2:<.5}(med)\n'.

format(i+2, mean_b, med_b))

f_out.write('\n')

for i in range(K-1):

mean_c, med_c = calc_mean_med(c_smpls[i])

f_out.write('c[{0}] = {1:<.5}(mean) {2:<.5}(med)\n'.

format(i+1, mean_c, med_c))

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の平均値と中央値で出力されている。

図Ｔ２

引用文献

Samejima, F. (1969). Estimation of latent ability using a response pattern of graded scores. Psychometric Monograph, No. 17, Richmond: Psychometric Society.

Up