Yasuharu Okamoto, 2017.7

Analysis of Polytomous Items by Item Response Theory

A Stan script with Python

 

The well-known model for polytomous items is the graded response model (Samejima, 1969). For binary items, visit this website.

Let the number of categories of items be K, where categories are denoted by integers 1 to K. Denote a response (category) on item  by , then

Set probabilities of  as follows:

where

Parameters s and s appear in equations in the form , and for any value , we have

That is, the model is not identifiable. To make the model identifiable, set

 

For the above model, the following Stan script is prepared.

 

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

}

=================================================================

 

A Python script using the above Stan script is this. The script files and an example data file were archived into a file PrgIRTpol.zip, which can be down loaded by clicking the name PrgIRTpol.zip.

 

=========================================================

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

=======================================================

 

Run this Python script, then file names of input data and output data are asked (Figure 1).

Figure 1

 

A file name of the output data is any text file name. The input data file should be prepared in the following format shown in Figure 2.

Figure 2

 

On the line next to the line with slash / at the head, the number of data (persons) is written (200 in Figure 2). On the next line to this line, the number of items is written (10 in Figure 2), and on the following line, the number of categories is written (6 in Figure 2). After these 3 lines, the responses of persons are written, one persons data per line. Each line starts with data ID, then responses are written with one or more than one spaces as seperators.

After the output file name is set, calculation starts. The numbers of persons, items, and categories are displayed, then calculation by Stan starts (Figure 3).

Figure 3

 

Figure 4 shows that Stan took about 5 minutes to sample in a chain for this condition. Total time of sampling is the sum of these times for chains when the chains are run sequentially.

Figure 4

 

When the program ends, statistics by Stan is displayed (Figure 5).

Figure 5

 

Parts of the content of the output file are shown in Figures 6 and 7.

Figure 6

 

On the first line, the input data file name is printed out. Then, the read in data are printed out, so that you can check whether the intended data are correctly read in or not. After the read in data, the point estimates of parameters, mean and medians, are printed out (Figure 7).

Figure 7

 

As shown by Equation (1), any shifted values  and  for arbitrary constant  can also be accepted. If we want to put the center category boundary  at the origin 0, set  and transform s and s as follows:

 

 

 

Reference

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

 

 

Home