Item Response Theory of Stage Theory in PyMC3

Item Response Theory of Stage Theory

A simple sample script in PyMC3

Yasuharu Okamoto

Item response theory (IRT) usually assumes that the ability is continuous variable, but in psychology, it is not rare to assume that ability is represented by stages. A famous stage theory is one by Piaget. Okamoto (2013) presented a model of stage theory by extension of item response theory. In this website, we present a simple script in PyMC3 of stage theory with a slight modification of Okamoto (2013)’s model. For general introduction to latent class modeling, see Levy and Mislevy (2016).

Scripts in this site were developed on Ubuntu. In Winodws, PyMC3 scripts runs on Jupyter Notebook. Using PyMC3 in Jupyter Notebook in Windows is easy and explained in this website.

In Item Response Model (IRT), probability of correct response is determined by the trait of a person:

So, the probability is represented as

But, for a latent stage model, this probability is represented as

Latent stage model assumes that probability of correct response increases as stage proceeds. That is,

(3) If Stage k -> Stage (k+1), then P(Correct|Stage k) <= P(Correct|Stage k+1)

Put the number of stages K. Denote the increment of probabilities of correct responses between stages as p_step, that is,

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)

Then we have

P(Correct|Stage k)=p_step[1]+…+p_step[k]

p_step[k] ≥ 0

p_step[1]+…+p_step[K+1]=1

In the case of K=3, we have the following matrix equation

Equation (4) is used in the following code where K is not restricted to be 3:

prob_response_ = pm.math.matrix_dot(coeff, prob_steps.T)

prob_step represents p_step, which is declared as an object of Dirichlet type as follows

prob_steps = pm.Dirichlet('prob_steps', a = a_Diri, shape = (M, NStg + 1))

where NStg and M denote the number of stages K and the number of items.

Probabilities of stages to which persons belong to are represented as an object of Categorical type and denoted by p_stages.

p_stages = pm.Categorical('p_stages', p = a_p, shape = N)

where N is the number of persons.

The complete script of the program in Listing 1 is a demonstration for a simple data, where K=3, M=10, and N=21.

More general version is presented later in this web site.

The compressed file IRTLClassPyMC3.zip of the script in Listing 1 can be freely downloaded.

Listing 1 A Sample program.

"""

Yasuharu Okamoto, 2019.02, 2020.03

"""

import matplotlib.pyplot as plt

import numpy as np

import pymc3 as pm

import arviz as az

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

([[1]*5 + [1]*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 # Number of stages

print('NStg = ', NStg)

a_Diri = np.ones((M, NStg + 1))

print('a_Diri =\n', a_Diri)

a_coeff = np.zeros((NStg, NStg + 1))

for i in range(NStg):

for j in range(i + 1):

a_coeff[i][j] = 1

print('a_coeff =\n', a_coeff)

a_p = np.ones(NStg)

print('a_p =\n', a_p)

with pm.Model() as stage_model:

prob_steps = pm.Dirichlet('prob_steps', a = a_Diri, shape = (M, NStg + 1))

coeff = pm.math.constant(a_coeff, 'coeff')

prob_response_ = pm.math.matrix_dot(coeff, prob_steps.T)

prob_response = pm.Deterministic('prob_response', prob_response_)

p_stages = pm.Categorical('p_stages', p = a_p, shape = N)

p_at_stage = pm.Deterministic('p_at_stage', prob_response[p_stages])

y = pm.Bernoulli('y', p = prob_response[p_stages], observed = X)

trace = pm.sample()

summary = pm.summary(trace)

print('summary...\n', summary)

ary_prob_response = trace['prob_response']

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

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

for i in range(NStg):

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

print(' ')

for j in range(M):

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

for s in range(NStg):

print('{0: >10.5f}'.format(ary_prob_response[:,s,j].mean()), end = '')

print(' ')

ary_p_stages = trace['p_stages']

ary_p_stages = ary_p_stages.T

print("\nMean stages and their SD's\n")

print('{0: >11s}{1:>10s}{2:>10s}'.format(' ', 'Mean', 'SD'))

for i in range(len(ary_p_stages)):

print('person-{0: <4d}{1:>10.3f}{2:>10.3f}'.format(i, ary_p_stages[i].mean(),

ary_p_stages[i].std()))

v_waic = az.waic(trace, scale = 'deviance')

print(v_waic)

print('\nWAIC = {0:.3f} for the number of stages = {1}'.format(v_waic['waic'], NStg))

Run the script in Listing 1, the following output will be displayed.

X =

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

[1, 1, 1, 1, 1, 0, 0, 0, 0, 0]

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

The data are shown in the form of 10 columns (items) and 21 rows (persons). Correct and wrong responses are denoted by 1 and 0, respectively.

After sampling by PyMC3, mean values of will be displayed as follows.

P(correct|item,stage)

Stage-0 Stage-1 Stage-2

Item-0 0.11299 0.83643 0.94046

Item-1 0.10882 0.83649 0.94223

Item-2 0.11266 0.83717 0.94057

Item-3 0.12034 0.83957 0.94359

Item-4 0.10975 0.83419 0.93951

Item-5 0.06003 0.16121 0.89432

Item-6 0.05926 0.16450 0.88910

Item-7 0.05962 0.16370 0.88480

Item-8 0.05887 0.16182 0.88767

Item-9 0.05980 0.16246 0.88645

In Stage-0, all probabilities of correct responses are low values. In Stage-1, probabilities of correct responses for Item-0 to Item-4 are high, but for Item-5 to Item-9 probabilities of correct responses are low. In Stage-2, all probabilities of correct responses are high.

The following output shows mean values of stages (stages are represented by integers), which persons are in.

Mean stages and their SD's

Mean SD

person-0 0.000 0.000

person-1 0.000 0.000

person-2 0.000 0.000

person-3 0.000 0.000

person-4 0.000 0.000

person-5 0.000 0.000

person-6 0.003 0.055

person-7 1.000 0.000

person-8 1.000 0.000

person-9 1.000 0.000

person-10 1.001 0.032

person-11 1.000 0.000

person-12 1.000 0.000

person-13 1.000 0.000

person-14 2.000 0.000

person-15 2.000 0.000

person-16 2.000 0.000

person-17 2.000 0.000

person-18 2.000 0.000

person-19 2.000 0.000

person-20 2.000 0.000

The results show that Person-0 to Person-6 are in Stage-0, Person-7 to Preson-13 in Stage-1, and Person-14 to Person-20 in Stage-2.

Summary statistics of sampling by PyMC3 are as follows.

Multiprocess sampling (2 chains in 2 jobs)

CompoundStep

>NUTS: [prob_steps]

>CategoricalGibbsMetropolis: [p_stages]

Sampling 2 chains, 0 divergences: 100%|█| 2000/2000 [00:06<00:00, 316.30draws/s]

summary...

mean sd hpd_3% ... ess_bulk ess_tail r_hat

p_stages[0] 0.000 0.000 0.000 ... 1000.0 1000.0 NaN

p_stages[1] 0.000 0.000 0.000 ... 1000.0 1000.0 NaN

p_stages[2] 0.000 0.000 0.000 ... 1000.0 1000.0 NaN

p_stages[3] 0.000 0.000 0.000 ... 1000.0 1000.0 NaN

p_stages[4] 0.000 0.000 0.000 ... 1000.0 1000.0 NaN

... ... ... ... ... ... ... ...

p_at_stage[20,5] 0.887 0.102 0.707 ... 625.0 430.0 1.0

p_at_stage[20,6] 0.888 0.096 0.719 ... 1307.0 575.0 1.0

p_at_stage[20,7] 0.892 0.095 0.720 ... 1020.0 494.0 1.0

p_at_stage[20,8] 0.888 0.099 0.709 ... 1061.0 417.0 1.0

p_at_stage[20,9] 0.892 0.092 0.723 ... 1029.0 693.0 1.0

[301 rows x 11 columns]

Script, which reads in data from a file

The script in Listing 1 contains data in it. The script in Listing 2 reads in data from a file.

Data file should be prepared in the format as shown in Figure 1, and be saved as a CSV file.

In the first row, labels are set, and in the first column case IDs set. Correct responses are coded as 1, and wrong responses as 0.

Figure 1

When saving the data file, a dialog box like shown in Figure 2 will be presented.

Figure 2

Click on the inverted triangle enclosed by a red circle in Figure 2, then a dialog box shown in Figure 3 will be presented.

Select the menu item “Text CSV(.csv)”, then the dialog box will become as shown in Figure 2.

Figure 3

The script which reads in data from a file is shown in the following Listing 2. The script file and a sample data file are archived into the compressed file IRTLClassPyMC3RD.zip, which can be freely downloaded. Files for Jupyter Notebook are also prepared, and archived into the compressed file IRTLClassPyMC3RD_jn.zip, which can also be freely downloaded.

Listing 2 Script for reading in data from a csv file.

"""

Yasuharu Okamoto, 2019.02, 2019.03

"""

import csv

import matplotlib.pyplot as plt

import numpy as np

import pymc3 as pm

import arviz as az

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

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

print('NStg = ', NStg)

a_Diri = np.ones((M, NStg + 1))

print('a_Diri =\n', a_Diri)

a_coeff = np.zeros((NStg, NStg + 1))

for i in range(NStg):

for j in range(i + 1):

a_coeff[i][j] = 1

print('a_coeff =\n', a_coeff)

a_p = np.ones(NStg)

print('a_p =\n', a_p)

with pm.Model() as stage_model:

prob_steps = pm.Dirichlet('prob_steps', a = a_Diri, shape = (M, NStg + 1))

coeff = pm.math.constant(a_coeff, 'coeff')

prob_response_ = pm.math.matrix_dot(coeff, prob_steps.T)

prob_response = pm.Deterministic('prob_response', prob_response_)

p_stages = pm.Categorical('p_stages', p = a_p, shape = N)

p_at_stage = pm.Deterministic('p_at_stage', prob_response[p_stages])

y = pm.Bernoulli('y', p = prob_response[p_stages], observed = X)

trace = pm.sample()

summary = pm.summary(trace)

print('summary...\n', summary)

fout.write('summary...\n{}'.format(summary))

ary_prob_response = trace['prob_response']

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

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

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

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

for i in range(NStg):

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

fout.write("{0:>10s}".format('Stage-{}'.format(i)))

print(' ')

fout.write('\n')

for j in range(M):

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

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

for s in range(NStg):

print('{0: >10.5f}'.format(ary_prob_response[:,s,j].mean()), end = '')

fout.write('{0: >10.5f}'.format(ary_prob_response[:,s,j].mean()))

print(' ')

fout.write('\n')

ary_p_stages = trace['p_stages']

ary_p_stages = ary_p_stages.T

print("\nMean stages and their SD's\n")

fout.write("\nMean stages and their SD's\n\n")

print('{0: >11s}{1:>10s}{2:>10s}'.format(' ', 'Mean', 'SD'))

fout.write('{0: >11s}{1:>10s}{2:>10s}\n'.format(' ', 'Mean', 'SD'))

for i in range(len(ary_p_stages)):

print('{0:>11s}{1:>10.3f}{2:>10.3f}'.format(CaseID[i],

ary_p_stages[i].mean(), ary_p_stages[i].std()))

fout.write('{0:>11s}{1:>10.3f}{2:>10.3f}\n'.format(CaseID[i],

ary_p_stages[i].mean(), ary_p_stages[i].std()))

v_waic = az.waic(trace, scale = 'deviance')

print('\nWAIC statistics...\n', v_waic, '\n')

fout.write('\nWAIC statistics...\n{}\n\n'.format(v_waic))

print('\nWAIC = {0:.3f} for the number of stages = {1}'.format(v_waic['waic'], NStg))

fout.write('\nWAIC = {0:.3f} for the number of stages = {1}\n'.format(v_waic['waic'], NStg))

print('Output file name =', fout_name)

fout.close()

Run the script in Listing 2, the data file name and output file name are asked to be set (Figure 4).

Figure 4

After setting the names, the numbers of persons and items are displayed.

Before starting sampling by PyMC3, the number of stages is asked to be set (Figure 5). Because the number of stages is set at execution, we can try various values for number of stages. In Figure5, 3 is set.

Figure 5

After sampling, results as shown in the following will be displayed.

P(correct|item,stage)

Stage-0 Stage-1 Stage-2

X1 0.10897 0.83819 0.94009

X2 0.10898 0.83755 0.94136

X3 0.10994 0.83308 0.93881

X4 0.11381 0.83792 0.94322

X5 0.10874 0.83533 0.94079

X6 0.05724 0.16366 0.88769

X7 0.06045 0.16525 0.88490

X8 0.06016 0.16684 0.88858

X9 0.06126 0.17013 0.88683

X10 0.05933 0.16181 0.88655

Mean stages and their SD's

Mean SD

psn1 0.000 0.000

psn2 0.000 0.000

psn3 0.000 0.000

psn4 0.000 0.000

psn5 0.007 0.083

psn6 0.000 0.000

psn7 0.000 0.000

psn8 1.002 0.045

psn9 1.000 0.000

psn10 1.000 0.000

psn11 1.000 0.000

psn12 1.000 0.000

psn13 1.000 0.000

psn14 1.000 0.000

psn15 2.000 0.000

psn16 2.000 0.000

psn17 1.995 0.071

psn18 2.000 0.000

psn19 2.000 0.000

psn20 2.000 0.000

psn21 2.000 0.000

WAIC statistics...

Computed from 1000 by 210 log-likelihood matrix

Estimate SE

deviance_waic 55.72 1.58

p_waic 2.79 -

WAIC = 55.724 for the number of stages = 3

Output file name = ResultsRD.txt

Fitness index WAIC is displayed. Using WAIC, we can compare model fitness for various values of stages.

The output file can be opened after execution of the script.

Reference

Levy, R. & Mislevy, R. J. (2016). Bayesian psychometric modeling. CRC Press

Okamoto, Y “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/)