Stan Script for Latent Class Item Response Theory
Yasuharu Okamoto, 2019.04
Item Response Theory (IRT) sets probability of response as being dependent on some trait, that is,
Responses is usually assumed to be binary, e.g., 1 (positive response +) or 0 (negative response -), and trait to be quantitative and continuous.
However, in some theory of psychology, trait is qualitative and discrete. For example, Piaget proposed Stage Theory of Development. In this stage theory, stages (latent classes) are ordered, and probabilities of eq. (1) are increasing functions of stages. For this ordered type of latent class IRT, sample scripts are presented at a website in Stan or a website in PyPC3.
In this website, latent classes are not assumed to be ordered, as explained below. Latent class models, in which classes are not ordered, are also adopted in psychology. For example, Levenson (2014) presents discrete view of emotions, and Meehl (2004) explains a taxon for schizophrenia.
Probability of response 1, which is denoted by +, on an item i by a person of a latent class c is represented as follows (cf. Levy & Mislevy, 2016)
The script in this website uses MCMC sampling, which needs some restriction to avoid label switching. The script of this website sets the following restriction:
Each item belongs to some class, that is, suppose that item i belong to class c. For this case, probability of response + is supposed to be highest for class c.
The Stan script in this website, data are prepared so that items are grouped according to the class to that they belongs. Grouping of items may be done based on related researches, or by principal components analysis (PCA) of items. An example of analysis of items by PCA is presented in the last section of this website. In any case, the first analysis by the script in this website would be considered as preparatory one, on which items would be re-selected and grouped. After this pre-analysis, the final analysis of the re-constructed items will be done by the same program in this website.
To develop the script in Stan, the following marginalization is employed.
The script in Stan is shown in Listing 1 below.
After parameters 
s are estimated, for a person with response pattern {
}, the posterior distribution of class is given by
where designates the response
on item i, and P(class) is a prior distribution for the population to which the
person belongs.
Note that
and
Listing 1. Stan script for the latent class IRT.
//
// Binary Responses
for NCls Latent Classes
//
Yasuharu Okamoto, 2019.04
//
data {
int N;
// Number of subjects
int NCls;
// Number of classes
int Msub[NCls]; // Number of items in each class
int M;
// Total number of items =
sum(Msub)
int X[N, M];
// Responses: 0 -> failure, 1 -> success
vector[NCls] a; // Parameter of the prior distribution of
the classes
}
transformed data{
int CumM[NCls];
CumM[1] = Msub[1];
for (s in 2:NCls){
CumM[s] = CumM[s-1] + Msub[s];
}
}
parameters {
simplex[NCls] theta[N];
real<lower = 0.0, upper =
1.0> pp[M, NCls];
}
transformed parameters {
real p[M, NCls];
real p_person[N, M];
vector[2] p_person_cat[N,
M];
for (i in 1:CumM[1]){
for
(s in 1:NCls){
if
(s == 1){
p[i][s]
= pp[i][1];
}
else{
p[i][s]
= pp[i][s] * pp[i][1];
}
}
}
for (s_g in 2:NCls){
for (i in (CumM[s_g-1]+1):CumM[s_g]){
for
(s in 1:NCls){
if
(s == s_g){
p[i][s]
= pp[i][s_g];
}
else{
p[i][s]
= pp[i][s] * pp[i][s_g];
}
}
}
}
for (j in 1:N) {
for
(i in 1:M){
p_person[j][i] = 0.0;
for (stg in 1:NCls){
p_person[j][i]
= p_person[j][i] + p[i][stg] * theta[j][stg];
}
}
}
}
model {
for (j in 1:N){
theta[j] ~ dirichlet(a);
}
for (i in 1:M){
for
(s in 1:NCls){
pp[i][s] ~ uniform(0.0, 1.0);
}
}
for (j in 1:N) {
for
(i in 1:M){
X[j][i] ~ bernoulli(p_person[j][i]);
}
}
}
The prior distribution of the classes, which is represented by theta in the Stan script in Listing 1, is given by Dirichlet distribution as follows,
for (j in 1:N){
theta[j] ~ dirichlet(a);
}
The Stan script in Listing 1 is used in the following Python script in Listing 2, where the file name of the Stan script in Listing 1 is LClass.stan.
Listing 2. Python script, which uses the Stan script in Listing 1.
"""
Yasuharu
Okamoto, 2019.04
"""
import csv
import pystan
import pickle
import matplotlib.pyplot as plt
import seaborn as sb
import numpy as np
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))
N = len(data) - 2
M = len(data[0]) - 1
ItemNames = data[0][1:]
ck_cnt = data[1][1:]
cnts = np.bincount(ck_cnt)
NCls = len(cnts) - 1
Msub = []
for s in range(NCls):
Msub.append(cnts[s+1])
X = np.empty((N, M), dtype = int)
CaseID = []
for i in range(2, N + 2):
CaseID.append(data[i][0])
for j in range(M):
X[i
- 2][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))
if M != sum(Msub):
print('\nError in number of
M, or M[]')
raise Exception()
print('\nNumber of items for each
class')
for j in range(NCls):
print('M[{0}] =
{1}'.format(j+1, Msub[j]))
fout.write('M[{0}] =
{1}\n'.format(j+1, Msub[j]))
print('NCls = ', NCls)
fout.write('NCls = {}'.format(NCls))
a = []
# Parameter of the prior
distribtuion of the classes
for j in range(NCls):
a.append(1.0);
stan_data = {'N': N, 'M': M, 'Msub':
Msub, 'NCls': NCls, 'a': a, 'X': X}
sm = pystan.StanModel(file =
'LClass.stan')
def f_init():
pp = []
pos_pp = -1
for s in range(NCls):
for
i in range(Msub[s]):
pp.append([])
pos_pp += 1
for si in range(NCls):
pp[pos_pp].append(0.8)
return dict({'pp': pp})
fit = sm.sampling(data = stan_data,
pars = ['p', 'theta'],
n_jobs = 1, init = f_init)
print(fit)
p_samples = fit['p']
p_item_class = np.empty((M, NCls))
for i in range(M):
for s in range(NCls):
p_item_class[i][s]
= p_samples[:,i,s].mean()
print('\n\nP(+|item,class)\n')
fout.write('\n\nP(+|item,class)\n\n')
print('{0: <9s}'.format('Item'),
end = '')
for s in range(NCls):
print('{0:>9s}'.format('Class-{}'.format(s+1)), end = '')
print()
fout.write('{0:
<9s}'.format('Item'))
for s in range(NCls):
fout.write('{0:>9s}'.format('Class-{}'.format(s+1)))
fout.write('\n')
for i in range(M):
print('{0:
<9s}'.format(ItemNames[i]), end = '')
fout.write('{0:
<9s}'.format(ItemNames[i]))
for s in range(NCls):
print('{0:>9.5f}'.format(p_item_class[i][s]), end = '')
fout.write('{0:>9.5f}'.format(p_item_class[i][s]))
print(' ')
fout.write('\n')
theta = fit['theta']
p_class = np.empty((N, NCls))
for j in range(N):
for s in range(NCls):
p_class[j][s] = theta[:,j,s].mean()
print('\n\nP(class|person)\n')
fout.write('\n\nP(class|person)\n\n')
print('{0: <10s}'.format('Person'),
end = '')
fout.write('{0:
<10s}'.format('Person'))
for s in range(NCls):
print('{0:>9s}'.format('class-{}'.format(s+1)),
end = '')
fout.write('{0:>9s}'.format('class-{}'.format(s+1)))
print()
fout.write('\n')
for j in range(N):
print('{0:
<10s}'.format(CaseID[j]), end = '')
fout.write('{0:
<10s}'.format(CaseID[j]))
for s in range(NCls):
print('{0:>9.5f}'.format(p_class[j][s]), end = '')
fout.write('{0:>9.5f}'.format(p_class[j][s]))
print(' ')
fout.write('\n')
fout.close()
print('\nOutput file {} was
saved.\n'.format(fout_name))
Files of scripts in Listing 1, and Listing 2 are archived as a zipped file GeneralSample.zip with the sample data file.
The prior distribution of the latent classes is the Dirichlet distribution with parameter values being all 1. The parameter values are given by
a = []
for j in range(NCls):
a.append(1.0);
This means that the prior distribution of the classes is non-informative. When some information about the distribution of the latent classes is available, you can set this information in the parameter. Brief explanation of Dirichlet distribution is presented at this website.
How to use this program is explained in the following, using PyStan2.17.1 on Python3.7.3 in Windows 10.
Enter Python command as follows (Figure 1), where LClassFLStan.py is the file name of the script in Listing 2.
Python
LClassFLStan.py
Figure 1
When the program starts, input data file name and output file name are asked. Input data file is prepared as a CSV file, which can be easily created by Excel, saving the file with a file name having file extension csv. An example of csv file is shown in Figure 2.
Figure 2
In the first row, variable names are set. The first column is for case identification. After the first column, names for item variables are set. Items are aggregated according to groups, to which they belong.. In Figure 2, items X1 to X5 are belong to class 1, items X6 to X10 to class 2, and items X11 to X15 to class 3. In the second row, numbers to which items are belong are set. Classes are numbered by consecutive inters. When the number of classes is K, then integers from 1 to K are used. Correspondence between classes and integers is arbitrary. However, integers in the second row of the datasheet must be in the ascending order from the left to the right. This means that order of groups of items are arbitrary, but integers in the second row must not be in descending order. An example of data file, which is the same as that in Figure 2 except the order of items, is presented in the later part of the website.
Individual data are set in rows after the second row in the format (CaseID, X1, …,X15). Responses are coded as 1 for positive responses and 0 for negative responses.
After setting the output file name as a text file, the program shows the input data in the form “X =”.
The number of cases and items are shown as N, and M (Figure 3).
The numbers of items in each class are displayed in the form of “M[class] =’, and the number of classes as NCls.
Then compiling of the Stan script starts.
Figure 3
When MCMC sampling ends, the following output is displayed.
P(+|item,class)
Item
Class-1 Class-2 Class-3
X1
0.86957 0.09526 0.09552
X2
0.87031 0.09407 0.09536
X3
0.86986 0.09638 0.09485
X4
0.87314 0.09552 0.09473
X5
0.87550 0.09411 0.09283
X6
0.09606 0.87285 0.09502
X7
0.09501 0.87515 0.09730
X8
0.09474 0.87236 0.09501
X9
0.09565 0.87481 0.09583
X10
0.09602 0.87205 0.09434
X11
0.09469 0.09580 0.87139
X12
0.09561 0.09603 0.87143
X13
0.09645 0.09479 0.87103
X14
0.09624 0.09527 0.86987
X15 0.09593 0.09514 0.87063
P(class|person)
Person class-1 class-2 class-3
1
0.81940 0.08976 0.09083
2
0.81776 0.09306 0.08918
3
0.81712 0.09074 0.09215
4
0.81704 0.09184 0.09112
5
0.81550 0.09169 0.09281
6
0.81737 0.09102 0.09161
7
0.81816 0.09058 0.09126
8
0.81905 0.08978 0.09117
9
0.81688 0.08935 0.09377
10
0.81889 0.09158 0.08954
11
0.09039 0.81836 0.09125
12
0.09222 0.81688 0.09090
13
0.09269 0.81408 0.09323
14
0.09212 0.81708 0.09080
15
0.09054 0.81837 0.09109
16
0.09219 0.81559 0.09222
17
0.09337 0.81575 0.09089
18
0.08968 0.81912 0.09120
19
0.09114 0.81759 0.09127
20
0.08858 0.81829 0.09313
21
0.09314 0.09143 0.81543
22
0.09130 0.09038 0.81832
23
0.09185 0.09114 0.81701
24
0.09179 0.09289 0.81532
25
0.09277 0.09094 0.81629
26
0.09077 0.09092 0.81831
27
0.09176 0.09030 0.81794
28
0.09338 0.09037 0.81625
29
0.09045 0.09182 0.81773
30
0.09076 0.09174 0.81749
Output file Results.txt was saved.
PS D:\yasuharu\www\LaCoocan\Python\en1\LClassIRTNoOrder\GeneralSample>
After heading P(+|item,class) , means
of the posterior distributions of 
are shown. Means of
these probabilities are highest in class 1 for items X to X5, in class 2 for
items X6 to X10, and in class 3 for items X11 to X15.
After heading P(class|person) , means
of the posterior distributions of 
, which are represented
in the script as theta[j][stg] , are shown. Means of persons 1 to 10 are highest for class 1,
means of persons 11 to 20 are highest for class 2, and means of persons 21 to
30 are highest for class 3.
These results are satisfactory.
Another
example of data set, in which groups of items of Figure 2 are permuted
Order of groups of items is arbitrary. The data file in Figure 2a is the same as that in Figure 2 except the order of groups of items.
Figure 2a
However, the integers in the second row must not be in descending order. So, the values in the second row are rewritten so that the integers are not in descending order. Run the Python script in Listing 2 for the data of Figure 2a, we have the following results.
Data file name(*.csv) = Data_a.csv
Output file name(*.txt) = Results_a.txt
X =
1
[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
2
[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
3
[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
4
[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
5
[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
6
[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
7
[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
8
[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
9
[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
10 [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
11 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]
12 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]
13 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]
14 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]
15 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]
16 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]
17 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]
18 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]
19 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]
20 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]
21 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]
22 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]
23 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]
24 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]
25 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]
26 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]
27 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]
28 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]
29 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]
30 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]
N = 30
M = 15
Number of items for each class
M[1] = 5
M[2] = 5
M[3] = 5
NCls = 3
INFO:pystan:COMPILING THE C++ CODE FOR
MODEL anon_model_dee74dfece458dbb4907ad69b4d4e9ad NOW.
.
.
.
P(+|item,class)
Item
Class-1 Class-2 Class-3
X6
0.87052 0.09513 0.09479
X7
0.87257 0.09632 0.09870
X8
0.87259 0.09339 0.09465
X9
0.87222 0.09547 0.09669
X10
0.87069 0.09469 0.09334
X11
0.09402 0.87234 0.09499
X12
0.09440 0.87058 0.09557
X13
0.09449 0.86981 0.09703
X14
0.09478 0.87117 0.09638
X15
0.09464 0.87158 0.09524
X1
0.09710 0.09714 0.87376
X2
0.09517 0.09518 0.87286
X3 0.09602 0.09622 0.87336
X4
0.09758 0.09536 0.87222
X5
0.09403 0.09504 0.87351
P(class|person)
Person class-1 class-2 class-3
1
0.09085 0.08969 0.81946
2
0.09192 0.09103 0.81705
3
0.09097 0.08927 0.81977
4
0.09118 0.09057 0.81825
5
0.09092 0.09299 0.81610
6
0.09331 0.09137 0.81532
7
0.09161 0.09152 0.81687
8
0.09183 0.09261 0.81556
9
0.09226 0.09194 0.81580
10
0.09161 0.09339 0.81500
11
0.81555 0.09203 0.09242
12
0.81794 0.09036 0.09170
13
0.81996 0.09074 0.08930
14
0.81653 0.09313 0.09035
15
0.81843 0.09078 0.09079
16
0.81405 0.09497 0.09098
17
0.81699 0.09204 0.09097
18
0.81755 0.09066 0.09179
19
0.81883 0.09129 0.08988
20
0.81800 0.09074 0.09126
21
0.09051 0.81821 0.09128
22
0.09187 0.81732 0.09081
23 0.09147 0.81862 0.08992
24
0.09281 0.81674 0.09045
25
0.09335 0.81631 0.09034
26
0.09263 0.81767 0.08970
27
0.09262 0.81616 0.09121
28
0.08983 0.81902 0.09115
29
0.09250 0.81778 0.08972
30
0.09328 0.81495 0.09177
Output file Results_a.txt was saved.
The results is essentially the same as for the data of Figure 2, except the class labels
Principal
Component Analysis of Items
Results of principal component analysis (PCA) of the items in the data shown in Figure 2, are as follows.
Scattergram of the pattern matrix is shown in Figure 4. We see three groups of items X1 to X5, items X6 to X10, and items X11 to X15, in each group items are completely overlapped.
Figure 4
These
overlapping of items are shown in the following output:
Pattern
matrix for the principal components =
comp.1 comp.2
X1 -0.00000 1.00000
X2 0.00000 1.00000
X3 -0.00000 1.00000
X4 -0.00000 1.00000
X5 -0.00000 1.00000
X6 -0.86603 -0.50000
X7 -0.86603 -0.50000
X8 -0.86603 -0.50000
X9 -0.86603 -0.50000
X10 -0.86603 -0.50000
X11 0.86603 -0.50000
X12 0.86603 -0.50000
X13 0.86603 -0.50000
X14 0.86603 -0.50000
X15 0.86603 -0.50000
Data in
Figure 2 is artificial one. In the case of real data, we would see rather
obscure grouping.
References
Levenson,
R. W. (2014). What is the underlying structure of emotions: An argument for
discrete emotions. In S. N.-Hoeksema, B. L. Fredrickson, G. R. Loftus, & C.
Lutz (2014). Atkinson & Hilgard”s Introduction to Psychology, 16th edition. Cengage Learning. pp. 406-407.
Levy, R.
& Mislevy, R. J. (2016). Bayesian
psychometric modeling. CRC Press
Meehl, P.
E. (2004). What’s in a taxon? Journal of
Abnormal Psychology, 2004, 113, 39-43.