Bayesian Analysis of
Reliability of Ordered Categorical Items Scales
with Correction
for Shrinkage
Yasuharu Okamoto,
2023.06, 2023.08
Reliability coefficient of a test score is given by the precision of the test score measurement (McDonald, 1999, p. 63), which is given by the following equation
Reliability
coefficient 
is defined based on the following
model:
where 
is examinee 
’s
score on the 
-th
item (
),
is the mean of ,
is the factor loading of item ,
is examinee ’s
factor score whose distribution is the standard normal one 
,
and 
is a residual whose distribution is a
normal one. The summed score
of examinee is given by the following equation.
where
is a true score
and
is an error.
In
the case of categorical items, the continuous variable 
is discretized to categorical
variables 
(Okamoto, 2013). That is, for category boundaries 
, we
have
where
and K is the number of categories.
Hence, we have
where 
is the cumulative distribution function
for the standard normal random variable with the convention 
and 
.
The likelihood
function for data 
is given by
and the posterior
distribution by
with
a prior distribution 
.
For
the estimates of parameters and 
, reliability
coefficient 
is calculated by eq. (1).
It
is well known that shrinkage may occur in Bayesian analysis. Correction for
shrinkage is done in the scripts shown in the website.
To
investigate relation between the observed score
and
the true latent value
,
correlation
coefficient R and coefficient of determination 
are calculated (Okamoto, 2016, 2017).
To
show relation between observed values 
and estimates of the factor 
visually, a scatter plot is displayed.
Reliability coefficient 
is also calculated.
Actual
installation of a model depends on the number of categories K (Okamoto,
2013). Hence, the scripts are developed for K=2, K=3,
and K>3, respectively.
A
simple way of installing Python with Anaconda is explained on this Website.
Scripts in the website were developed in PyStan 3. For information about how to install and use PyStan 3, check the following websites
http://y-okamoto-psy1949.la.coocan.jp/Python/en2/wsl_ubuntu/
and
http://y-okamoto-psy1949.la.coocan.jp/Python/en2/TryingPyStan3/
Files used at the website are archived in the files files8_G3Cat.zip for K>3, files8_3Cat.zip for K=3, and files8_2Cat.zip for K=2. They can be freely downloaded and used.
In this case of K>3, the Stan script for sampling parameter values from posterior distribution is as follows (see Okamoto, 2013, p. 156):
data {
int N;
int M;
int K;
array[N, M] int<lower =
1, upper = K> X;
vector[K-2] a;
}
transformed data {
real C1;
real CK_1;
//
// These values specify
the origin and unit of the scale
// cf. Okamoto (2013).
Behaviormetrika, 40, 149-168.
//
C1 = -1;
CK_1 = 1;
}
parameters {
array[M] real mu;
array[M] real<lower =
0.0> Lmd;
array[M] real<lower =
0.0> psy;
simplex[K-2] Ck;
array[N] real F;
}
transformed parameters {
array[N, M] vector[K] p;
array[N, M] vector[K-1]
cumP;
array[K-1] real C;
C[1] = C1;
C[K-1] = CK_1;
for (k in 1:(K-3)) {
C[k+1] = C[k] + Ck[k] * 2.0;
}
for (i in 1:N)
for
(j in 1:M){
cumP[i][j][K-1] = normal_cdf((C[K-1] - (mu[j] + Lmd[j] * F[i])) / psy[j]
| 0.0, 1.0);
cumP[i][j][1] = normal_cdf((C[1] - (mu[j] + Lmd[j] * F[i])) / psy[j] |
0.0, 1.0);
for (k in 1:(K-3)) {
cumP[i][j][k+1] = normal_cdf((C[k+1] - (mu[j] + Lmd[j] * F[i])) / psy[j]
| 0.0, 1.0);
}
p[i][j][1] = cumP[i][j][1];
for (k in 1:(K-2)) {
p[i][j][k+1] = cumP[i][j][k+1] - cumP[i][j][k];
}
p[i][j][K] = 1.0 - cumP[i][j][K-1];
}
}
model {
for (j in 1:M) {
mu[j] ~ normal(0.0, 10.0);
psy[j] ~
exponential(0.1);
Lmd[j] ~ exponential(0.1);
}
Ck ~ dirichlet(a);
F ~ normal(0.0, 1.0);
for (i in 1:N) {
for
(j in 1:M) {
X[i][j] ~ categorical(p[i][j]);
}
}
}
generated quantities {
real sdF;
real mnF;
real sumF;
real ssF;
array[N] real cF;
array[M] real cLmd;
sumF = 0.0;
ssF = 0.0;
for (i in 1:N) {
sumF
+= F[i];
ssF
+= pow(F[i], 2);
}
mnF = sumF / N;
sdF = pow((ssF - N *
pow(mnF, 2)) / (N - 1), 0.5);
for (i in 1:N) {
cF[i] = (F[i] - mnF) / sdF;
}
for (j in 1:M) {
cLmd[j] = Lmd[j] * sdF;
}
}
The Python script, which uses this Stan script to analyze the data, is shown in Listing 1.
When you run the script on Windows, you must install Ubuntu on Windows, and run the script on Ubuntu on Windows. That is, you run the script in a Python 3.11 virtual environment of Anaconda on Ubuntu (WSL: Windows Subsystem for Linux) on Windows. For more information, check this website, and this website.
Run the script Reli_G3Cat.py in Listing 1, then the name of the input data file is asked. Set the input file name as follows.
(py311) … /files8_G3Cat$ python Reli_G3Cat.py
Data File Name (*.xlsx) = DataN5.xlsx
The input data file is expected to be an Excel data file (*.xlsx) (Figure 1.1).
Figure 1.1
In the first row, the names of the variables (items) are set. The first column is used for case identification, and any value (string) is OK. In the second row, set the number of categories K in the second column, i.e., the cell B2. In Figure 1.1, “5” is set. From the third row, responses on the categorical items are set. Values on the items should be integer values from 1 to K, where K is the number of categories. The data in Figure 1.1 is for the case of K=5.
After the input data file name is set,
sampling by Stan starts. When the sampling by Stan ends, posterior
distributions of 
s are displayed (Figure 1.2).
Figure 1.2
Click the X icon at the upper right corner of
the windows to close the window. The next window shows posterior distributions
of 
s (Figure 1.3).
Figure 1.3
Close the window in Figure 1.3, then the
window in Figure 1.4 appears and the posterior distributions of 
s are presented.
Figure 1.4
Close the window in Figure 1.4, then the
next window appears and the posterior distribution of is presented (Figure
1,5).
Figure 1.5
Close the window in Figure 1.5, then the
window, where the posterior distributions of
R and are drawn, appears
(Figure 1.6). 
is the correlation
coefficient of the observed values 
s and the true values s, and is the square of and represents the
coefficient of determination of s and s (Okamoto, 2016, 2017).
Figure 1.6
When the window in Figure 1.6 is closed, a
scatter gram of observed scores and estimates of is displayed (Figure
1.7).
Figure 1.7
Close the window of Figure 1.7, the script ends.
After the script ends, the following message is shown on the terminal.
Results.txt was saved.
(py311) … /files8_G3Cat$
The contents of the file Results.txt are as follows.
Data file = DataN5.xlsx
K = 5, N = 100, M = 10
1: [1 1 2 1 2 2 2 2 2 3]
2: [1 1 1 1 1 1 1 1 2 2]
3: [2 2 1 2 3 2 3 2 3 3]
.
.
.
98: [2 3 3 4 3 5 4 3 4 4]
99: [1 1 2 2 2 3 1 2 3 3]
100: [1 1 1 1 1 1 3 2 2 3]
Summary...
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu[0] -0.611 0.077 -0.758 -0.470 0.005 0.003 247.0 525.0 1.01
mu[1] -0.542 0.086 -0.702 -0.384 0.005 0.004 247.0 543.0 1.02
mu[2] -0.380 0.086 -0.548 -0.227 0.006 0.004 238.0 467.0 1.02
.
.
.
cF[98] -0.850 0.181 -1.175 -0.494 0.002 0.001 8479.0 3107.0 1.00
cF[99] -1.266 0.186 -1.627 -0.938 0.002 0.001 8960.0 2995.0 1.00
sdF 0.948 0.070 0.827 1.083 0.004 0.003 348.0 677.0 1.01
mnF 0.010 0.101 -0.174 0.197 0.008 0.005 175.0 345.0 1.02
mu[1] = -0.6118
mu[2] = -0.54057
mu[3] = -0.37889
mu[4] = -0.2674
mu[5] = -0.14473
mu[6] = 0.0005217
mu[7] = 0.090897
mu[8] = 0.28648
mu[9] = 0.48015
mu[10] = 0.44828
Lambda[1] = 0.59608
Lambda[2] = 0.67004
Lambda[3] = 0.67247
Lambda[4] = 0.76804
Lambda[5] = 0.72232
Lambda[6] = 0.75365
Lambda[7] = 0.65455
Lambda[8] = 0.78646
Lambda[9] = 0.74405
Lambda[10] = 0.67585
Psy[1] = 0.30844
Psy[2] = 0.34215
Psy[3] = 0.32863
Psy[4] = 0.35737
Psy[5] = 0.317
Psy[6] = 0.34649
Psy[7] = 0.30951
Psy[8] = 0.43274
Psy[9] = 0.33779
Psy[10] = 0.27708
omega(mean) = 0.97692
omega(Med.) = 0.97699
omega(Q1) = 0.97572 omega(Q3) = 0.97822
R^2(Med.) = 0.9605
R^2(Q1) = 0.95609 R^2(Q3) = 0.96455
R(Med.) = 0.98005
R(Q1) = 0.9778 R(Q3) = 0.98212
alpha = 0.96631
standardized alpha = 0.96836
When K=3, no parameters between 
and 
, which are set to be -1 and 1, respectively (see, Okamoto, 2013, p. 156), are present. The Stan script
is as follows:
data {
int N;
int M;
array[N, M] int<lower =
1, upper = 3> X;
}
transformed data {
real C1;
real C2;
//
// These values specify
the origin and unit of the scale
// cf. Okamoto (2013).
Behaviormetrika, 40, 149-168.
//
C1 = -1;
C2 = 1;
}
parameters {
array[M] real mu;
array[M] real<lower =
0.0> Lmd;
array[M] real<lower =
0.0> vpsy;
array[N] real F;
}
transformed parameters {
array[M] real psy;
array[N, M] vector[3] p;
array[N, M] real cumP;
for (j in 1:M) {
psy[j] = 0.01 + vpsy[j];
}
for (i in 1:N)
for
(j in 1:M){
cumP[i][j] = normal_cdf((C2 - (mu[j] + Lmd[j] * F[i])) / psy[j] | 0.0,
1.0);
p[i][j][1] = normal_cdf((C1 - (mu[j] + Lmd[j] * F[i])) / psy[j] | 0.0,
1.0);
p[i][j][2] = cumP[i][j] - p[i][j][1];
p[i][j][3]
= 1.0 - cumP[i][j];
}
}
model {
for (j in 1:M) {
mu[j] ~ normal(0.0, 10.0);
vpsy[j] ~ exponential(0.1);
Lmd[j] ~ exponential(0.1);
}
F ~ normal(0.0, 1.0);
for (i in 1:N) {
for (j
in 1:M) {
X[i][j] ~ categorical(p[i][j]);
}
}
}
generated quantities {
real sdF;
real mnF;
real sumF;
real ssF;
array[N] real cF;
array[M] real cLmd;
sumF = 0.0;
ssF = 0.0;
for (i in 1:N) {
sumF
+= F[i];
ssF
+= pow(F[i], 2);
}
mnF = sumF / N;
sdF = pow((ssF - N *
pow(mnF, 2)) / (N - 1), 0.5);
for (i in 1:N) {
cF[i] = (F[i] - mnF) / sdF;
}
for (j in 1:M) {
cLmd[j] = Lmd[j] * sdF;
}
}
The Python script, which uses this Stan script to analyze the data, is shown in Listing 2.
When you run the script on Windows, you must install Ubuntu on Windows, and run the script on Ubuntu on Windows. That is, you run the script in a Python 3.11 virtual environment of Anaconda on Ubuntu (WSL: Windows Subsystem for Linux) on Windows. For more information, check this website, and this website.
Run the script Reli_3Cat.py in Listing 2, the input file name is asked as follows.
(py311) … /files8_3Cat$ python Reli_3Cat.py
Data File Name (*.xlsx) = Data.xlsx
The input data file is expected to be an Excel data file (*.xlsx) (Figure 2.1).
・
・
・
Figure 2.1
In the first row, the names of the
variables (Items) are set. The first column is used for case identification,
and any value (string) is OK. From the second column, responses on the
categorical items are set. Values on the items should be integer values from 1
to K, where K is the number of categories, that is, K=3. Hence, category values
are 1, 2, or 3. After the input file name is set, sampling by Stan starts. When
the sampling by Stan ends, posterior distributions of s are displayed (Figure 2.2).
Figure 2.2
Click the X icon at the upper right corner
of the windows to close the window. The next window shows posterior
distributions of s are displayed (Figure 2.3).
Figure 2.3
Close the window in Figure 2.3, then the window
in Figure 2.4 appears and the posterior distributions of s are presented.
Figure 2.4
Close the window in Figure 2.4, then the
next window appears and the posterior distribution of is presented (Figure
2.5).
Figure 2.5
Close the window in Figure 2.5, then the
window, where the posterior distributions of
R and are drawn, appears
(Figure 2.6). is the correlation
coefficient of the observed values s and the true values s, and is the square of and represents the
coefficient of determination of s and s (Okamoto, 2016, 2017).
Figure 2.6
When the window in Figure 2.6 is closed, a
scatter gram of observed scores s and estimates of s is shown (Figure 2.7).
Figure 2.7
Close the window of Figure 2.7, then the program ends. After the program ends, the terminal display shows the following messages.
Results.txt was saved.
(py311) … /files8_3Cat$
The contents of the file Results.txt are as follows.
Data file = Data.xlsx
VarLabels =
['X1' 'X2' 'X3' 'X4' 'X5' 'X6' 'X7'
'X8' 'X9' 'X10']
N = 100, M = 10
Data...
1: [1 1 1 2 3 2 2 2 3 2]
2: [1 1 2 2 2 2 2 3 3 1]
3: [1 1 2 2 2 2 2 3 3 3]
.
.
.
98: [2 3 2 3 3 3 3 3 3 3]
99: [1 2 2 1 2 2 2 2 3 3]
100: [1 1 1 1 2 1 2 2 2 2]
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu[0] -1.783 0.317 -2.358 -1.195 0.016 0.012 341.0 892.0 1.01
mu[1] -1.588 0.358 -2.258 -0.930 0.020 0.014 316.0 724.0 1.01
mu[2] -0.623 0.240 -1.065 -0.177 0.016 0.011 220.0 534.0 1.02
.
.
.
cF[97] 1.104 0.290 0.578 1.649 0.004 0.003 5397.0 3141.0 1.00
cF[98] -0.089 0.182 -0.438 0.241 0.002 0.003 6161.0 2973.0 1.00
cF[99] -0.610 0.181 -0.954 -0.275 0.002 0.002 6409.0 2626.0 1.00
mu[1] = -1.7546
mu[2] = -1.5558
mu[3] = -0.61515
mu[4] = -0.50779
mu[5] = -0.085098
mu[6] = 0.19538
mu[7] = 0.77405
mu[8] = 1.1675
mu[9] = 2.5391
mu[10] = 2.1659
Lambda[1] = 1.8078
Lambda[2] = 2.0748
Lambda[3] = 1.8036
Lambda[4] = 2.0777
Lambda[5] = 2.2282
Lambda[6] = 2.4247
Lambda[7] = 2.502
Lambda[8] = 2.5311
Lambda[9] = 3.5708
Lambda[10] = 2.605
Psy[1] = 0.85983
Psy[2] = 1.2953
Psy[3] = 1.079
Psy[4] = 1.2426
Psy[5] = 0.95208
Psy[6] = 0.73502
Psy[7] = 0.88215
Psy[8] = 1.1486
Psy[9] = 1.631
Psy[10] = 1.1614
omega(mean) = 0.97659
omega(Med.) = 0.97688
omega(Q1) = 0.9743 omega(Q3) = 0.97915
R^2(Med.) = 0.89435
R^2(Q1) = 0.88037 R^2(Q3) = 0.9065
R(Med.) = 0.9457
R(Q1) = 0.93828 R(Q3) = 0.9521
alpha = 0.93625
standardized alpha = 0.93771
When K=2, only one parameter for the category
boundary is present and set to be 0. Parameters 
s and 
s are set as follows by technical reasons (see Okamoto, 2013, p. 158).
and 
.
The Stan script is as follows:
data {
int N;
int M;
array[N, M] int<lower =
1, upper = 2> X;
}
parameters {
array[M] real mu;
real<lower =
0.0> psy;
array[N] real F;
}
transformed parameters {
array[N, M] real p;
real sdF;
real mnF;
real sumF;
real ssF;
array[N] real cF;
sumF = 0.0;
ssF = 0.0;
for (i in 1:N) {
sumF
+= F[i];
ssF
+= pow(F[i], 2);
}
mnF = sumF / N;
sdF = pow((ssF - N *
pow(mnF, 2)) / (N - 1), 0.5);
for (i in 1:N) {
cF[i] = (F[i] - mnF) / sdF;
}
for (i in 1:N)
for
(j in 1:M){
p[i][j] = normal_cdf((mu[j] + cF[i]) / psy | 0.0, 1.0);
}
}
model {
psy ~ exponential(0.1);
mu ~ normal(0.0, 10.0);
F ~ normal(0.0, 1.0);
for (i in 1:N) {
for
(j in 1:M) {
(X[i][j] - 1) ~ bernoulli(p[i][j]);
}
}
}
The Python script, which uses this Stan script to analyze the data, is shown in Listing 3 in the latter half of this website.
When you run the script on Windows, you must install Ubuntu on Windows, and run the script on Ubuntu on Windows. That is, you run the script in a Python 3.11 virtual environment of Anaconda on Ubuntu (WSL: Windows Subsystem for Linux) on Windows. For more information, check this website, and this website.
Run the script Reli_2Cat.py in Listing 3, the input file name is asked as follows.
(py311) … /files8_2Cat$ python Reli_2Cat.py
Data File Name (*.xlsx) = Data.xlsx
The input data file is expected to be an Excel data file (*.xlsx) (Figure 3.2).
・
・
・
Figure 3.1
In the first row, the names of the
variables (Items) are set. The first column is used for case identification,
and any value (string) is OK. From the second column, responses on the
categorical items are set. Values on the items should be integer values from 1
to K, where K is the number of categories. In this case of K=2, category values
are 1 or 2. After the input file name is set, sampling by Stan starts. When the
sampling by Stan ends, posterior distributions of s are displayed (Figure 3.2).
Figure 3.2
Click the X icon at the upper right corner of
the windows to close. The window which displays the posterior distribution of appears (Figure 3.3).
Figure 3.3
Close the windows of Figure 3.3, The window
in Figure 3.4 appears and the posterior distributions of is presented .
Figure 3.4
Close the window in Figure 3.4, then the
window, where the posterior distributions of
R and are drawn, appears (Figure
3.5). is the correlation
coefficient of the observed values s and the true values s, and is the square of and represents the
coefficient of determination of s and s (Okamoto, 2016, 2017).
Figure 3.5
Close the window in Figure 3.5, then a
scatter plot of observed scores and estimates of is displayed (Figure
3.6).
Figure 3.6
Close the window of Figure 3.6, then the program ends.
After the program ends, the terminal display shows the following messages.
Results.txt was saved.
(py311) … /files8_2Cat$
The contents of the file Results.txt are as follows.
Data file = Data.xlsx
VarLables =
{VarLables}
N = 100, M = 10
X =
1 [1 1 2 1 1 1 2 2 2 2]
2 [1 2 1 2 1 1 1 2 2 2]
3 [1 1 1 1 1 1 2 1 2 1]
.
.
.
98 [2 2 1 1 1 2 2 1 2 2]
99 [1 1 1 1 1 1 1 1 2 1]
100 [1 1 1 1 1 1 1 1 1 2]
Summary...
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu[0] -1.192 0.137 -1.450 -0.929 0.002 0.001 6084.0 3189.0 1.00
mu[1] -1.059 0.134 -1.307 -0.811 0.002 0.001 6591.0 3170.0 1.00
mu[2] -0.736 0.121 -0.959 -0.513 0.001 0.001 6776.0 3268.0 1.00
.
.
.
cF[99] -0.928 0.375 -1.700 -0.287 0.004 0.004 8655.0 2869.0 1.00
mnF -0.002 0.097 -0.177 0.186 0.003 0.002 1132.0 1856.0 1.00
sdF 0.999 0.070 0.870 1.128 0.002 0.001 1788.0 2445.0 1.00
mu[1] = -1.1871
mu[2] = -1.0553
mu[3] = -0.73387
mu[4] = -0.26835
mu[5] = -0.23873
mu[6] = 0.028156
mu[7] = 0.15127
mu[8] = 0.32501
mu[9] = 0.48954
mu[10] = 0.84849
Psy = 0.63552
omega(mean) = 0.96078
omega(Med.) = 0.96118
omega(Q1) = 0.95694 omega(Q3) = 0.96507
R^2(Med.) = 0.85445
R^2(Q1) = 0.83922 R^2(Q3) = 0.8688
R(Med.) = 0.92436
R(Q1) = 0.91609 R(Q3) = 0.9321
alpha = 0.88414
standardized alpha = 0.88536
Listing 1 Python script in the case of K>3 (Reli_G3Cat.py). The files are contained in the archived file files8_G3Cat.zip , which can be freely downloaded and used.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
import stan
import seaborn as sb
import arviz as az
pd.options.display.max_rows = 10000
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
data = pd.read_excel(xlsx)
fout = open('Results.txt', 'w')
fout.write(f'Data file = {flnm}\n')
VarLabels = data.columns[1:].values;
print(VarLabels)
print('data.vlaues...\n',
data.values[:5])
K = int(data.values[0,1])
print('K =', K)
X = np.array(data.values[1:, 1:],
dtype='int')
N = len(X)
M = len(X[0])
print('N = ', N)
print('M = ', M)
print('X =\n', X)
fout.write('\n')
fout.write(f'K = {K}, N = {N}, M = {M}\n')
for i, v in enumerate(X):
fout.write(f'{i+1:>5d}: {v}\n')
Data = {'N': N, 'M': M, 'X': X, 'K':
K, 'a': [1.0]*(K-2)}
def f_init():
return dict(mu = [0.0]*M,
psy = [1.0]*M, F = [0.0]*N)
with open('Reli_G3Cat.stan', 'r') as
fstan:
sm =
stan.build(fstan.read(), data = Data)
fit = sm.sample(num_chains = 4, init =
[f_init()]*4)
print(fit)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data, var_names =
['mu', 'cLmd', 'psy', 'Ck', 'cF',
'sdF', 'mnF'] )
print('Summary =\n',smry)
fout.write('\n\n')
fout.write('Summary...\n')
fout.write(smry.__str__())
fout.write('\n')
d_fm = fit.to_frame()
Mus = []
for j in range(M):
Mus.append(d_fm[f'mu.{j+1}'])
Mus = np.array(Mus).T
for i in range(M):
sb.kdeplot(Mus.T[i])
plt.title('$\mu$', fontsize = 16)
plt.show()
print()
fout.write('\n')
for j in range(M):
print('mu[{0:}] =
{1:.5}'.format(j+1, np.percentile(Mus.T[j], 50)))
fout.write('mu[{0:}] =
{1:.5}\n'.format(j+1, np.percentile(Mus.T[j], 50)))
Lmds = []
for j in range(M):
Lmds.append(d_fm[f'cLmd.{j+1}'])
Lmds = np.array(Lmds).T
for i in range(M):
sb.kdeplot(Lmds.T[i])
plt.title('$\lambda$', fontsize = 24)
plt.show()
print()
fout.write('\n')
for j in range(M):
print('Lambda[{0:}] =
{1:.5}'.format(j+1, np.percentile(Lmds.T[j], 50)))
fout.write('Lambda[{0:}] =
{1:.5}\n'.format(j+1, np.percentile(Lmds.T[j], 50)))
Psys = []
for j in range(M):
Psys.append(d_fm[f'psy.{j+1}'])
Psys = np.array(Psys).T
for j in range(M):
sb.kdeplot(Psys.T[j])
plt.title('$\psi$', fontsize = 24)
plt.show()
print()
fout.write('\n')
for j in range(M):
print('Psy[{0:}] =
{1:.5}'.format(j+1, np.percentile(Psys.T[j], 50)))
fout.write('Psy[{0:}] =
{1:.5}\n'.format(j+1, np.percentile(Psys.T[j], 50)))
rho = [] # omega
for i in range(len(Lmds)):
ss_Lmd = Lmds[i].sum() ** 2
s_Psys2 = (Psys[i] **
2).sum()
rho.append(ss_Lmd / (ss_Lmd
+ s_Psys2))
rho_Q1, rho_Med, rho_Q3 =
np.percentile(rho, [25, 50, 75])
rho_Mean = np.array(rho).mean()
fout.write('\n')
print('\nomega(mean) =
{0:.5}'.format(rho_Mean))
fout.write('\nomega(mean) =
{0:.5}\n'.format(rho_Mean))
print('omega(Med.) =
{0:.5}'.format(rho_Med))
fout.write('omega(Med.) =
{0:.5}\n'.format(rho_Med))
print('omega(Q1) = {0:.5} omega(Q3) =
{1:.5}'.format(rho_Q1, rho_Q3))
fout.write('omega(Q1) = {0:.5} omega(Q3) =
{1:.5}\n'.format(rho_Q1, rho_Q3))
sb.kdeplot(rho)
plt.title('Mean = {0:.3}\nQ1 = {1:.3},
Med.= {2:.3}, Q3 = {3:.3}'.
format(rho_Mean, rho_Q1, rho_Med, rho_Q3), fontsize = 14)
plt.xlabel('$\omega$', fontsize = 16)
plt.show()
Fs = []
for i in range(N):
Fs.append(d_fm[f'cF.{i+1}'])
Fs = np.array(Fs).T
Xtot = X.sum(axis = 1)
Rs = []
R2s = []
for t in range(len(Fs)):
r = np.corrcoef(Xtot,
Fs[t])[0][1]
Rs.append(r)
R2s.append(r**2)
R_Q1, R_med, R_Q3 = np.percentile(Rs,
[25, 50, 75])
R2_Q1, R2_med, R2_Q3 =
np.percentile(R2s, [25, 50, 75])
sb.kdeplot(Rs, label = 'R')
sb.kdeplot(R2s, label = '$R^2$')
plt.title('R(Med.) = {0:.3} $R^2$(Med.) =
{1:.3}'.format(R_med, R2_med), fontsize = 16)
plt.legend()
plt.show()
FsMeds = np.median(Fs, axis=0)
plt.plot(FsMeds, Xtot, 'o')
plt.xlabel('F', fontsize = 16)
plt.ylabel('X', fontsize = 16)
plt.show()
fout.write('\n')
print('\nR^2(Med.) =
{0:.5}'.format(R2_med))
fout.write('\nR^2(Med.) =
{0:.5}\n'.format(R2_med))
print('R^2(Q1) = {0:.5} R^2(Q3) =
{1:.5}'.format(R2_Q1, R2_Q3))
fout.write('R^2(Q1) = {0:.5} R^2(Q3) =
{1:.5}\n'.format(R2_Q1, R2_Q3))
print('R(Med.) =
{0:.5}'.format(R_med))
fout.write('R(Med.) =
{0:.5}\n'.format(R_med))
print('R(Q1) = {0:.5} R(Q3) =
{1:.5}'.format(R_Q1, R_Q3))
fout.write('R(Q1) = {0:.5} R(Q3) =
{1:.5}'.format(R_Q1, R_Q3))
S = np.cov(X, rowvar = False)
VarX = np.var(X.sum(axis = 1))
alpha_coef = (M/(M-1)) * (1 -
(np.diag(S).sum() / VarX))
print('\nalpha =
{0:.5}'.format(alpha_coef))
fout.write('\n')
fout.write('\nalpha =
{0:.5}\n'.format(alpha_coef))
Cor = np.corrcoef(X, rowvar = False)
r_bar = (Cor.sum() -
np.diag(Cor).sum()) / (M * (M - 1))
std_alpha = M * r_bar / (1.0 + (M - 1)
* r_bar)
print('standardized alpha =
{0:.5}'.format(std_alpha))
fout.write('standardized alpha =
{0:.5}'.format(std_alpha))
fout.close()
print('Results.txt was saved.')
Listing 2 Python script in the case of K=3. The files are contained in the archived file files8_3Cat.zip , which can be freely downloaded and used.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
import stan
import seaborn as sb
import arviz as az
pd.options.display.max_rows = 10000
fout = open('Results.txt', 'w')
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
data = pd.read_excel(xlsx)
fout.write(f'Data file = {flnm}\n')
VarLabels = data.columns[1:].values;
print(VarLabels)
fout.write(f'VarLabels =
\n{VarLabels}\n')
X = data.values[:, 1:]
N = len(X)
M = len(X[0])
print('N = ', N)
print('M = ', M)
fout.write(f'\nN = {N}, M = {M}\n')
fout.write('\n\nData...\n')
for i, v in enumerate(X):
fout.write(f'{i+1:>5d}: {v}\n')
Data = {'N': N, 'M': M, 'X': X}
def f_init():
return dict(mu = [0.0]*M,
vpsy = [1.0]*M, F = [0.0]*N)
with open('Reli_3Cat.stan', 'r') as
fstan:
sm =
stan.build(fstan.read(), data = Data)
fit = sm.sample(num_chains = 4, init =
[f_init()]*4)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data, var_names =
['mu', 'Lmd', 'psy', 'sdF', 'mnF', 'cF'])
print('Summary =\n',smry)
fout.write('\n')
fout.write(f'{smry}\n')
d_fm = fit.to_frame()
print(fit)
Mus = []
for j in range(M):
Mus.append(d_fm[f'mu.{j+1}'])
Mus = np.array(Mus).T
for i in range(M):
sb.kdeplot(Mus.T[i])
plt.title('$\mu$', fontsize = 16)
plt.show()
print()
fout.write('\n')
for j in range(M):
print('mu[{0:}] =
{1:.5}'.format(j+1, np.percentile(Mus.T[j], 50)))
fout.write('mu[{0:}] =
{1:.5}\n'.format(j+1, np.percentile(Mus.T[j], 50)))
Lmds = []
for j in range(M):
Lmds.append(d_fm[f'cLmd.{j+1}'])
Lmds = np.array(Lmds).T
for i in range(M):
sb.kdeplot(Lmds.T[i],)
plt.title('$\lambda$', fontsize = 24)
plt.show()
print()
fout.write('\n')
for j in range(M):
print('Lambda[{0:}] =
{1:.5}'.format(j+1, np.percentile(Lmds.T[j], 50)))
fout.write('Lambda[{0:}] =
{1:.5}\n'.format(j+1, np.percentile(Lmds.T[j], 50)))
Psys = []
for j in range(M):
Psys.append(d_fm[f'psy.{j+1}'])
Psys = np.array(Psys).T
for j in range(M):
sb.kdeplot(Psys.T[j])
plt.title('$\psi$', fontsize = 24)
plt.show()
print()
fout.write('\n')
for j in range(M):
print('Psy[{0:}] =
{1:.5}'.format(j+1, np.percentile(Psys.T[j], 50)))
fout.write('Psy[{0:}] =
{1:.5}\n'.format(j+1, np.percentile(Psys.T[j], 50)))
rho = []
for i in range(len(Lmds)):
ss_Lmd = Lmds[i].sum() ** 2
s_Psys2 = (Psys[i] **
2).sum()
rho.append(ss_Lmd / (ss_Lmd
+ s_Psys2))
rho_Q1, rho_Med, rho_Q3 =
np.percentile(rho, [25, 50, 75])
rho_Mean = np.array(rho).mean()
fout.write('\n')
print('\nomega(mean) =
{0:.5}'.format(rho_Mean))
fout.write('\nomega(mean) =
{0:.5}\n'.format(rho_Mean))
print('omega(Med.) =
{0:.5}'.format(rho_Med))
fout.write('omega(Med.) =
{0:.5}\n'.format(rho_Med))
print('omega(Q1) = {0:.5} omega(Q3) =
{1:.5}'.format(rho_Q1, rho_Q3))
fout.write('omega(Q1) = {0:.5} omega(Q3) =
{1:.5}\n'.format(rho_Q1, rho_Q3))
sb.kdeplot(rho)
plt.title('Mean = {0:.3}\nQ1 = {1:.3},
Med.= {2:.3}, Q3 = {3:.3}'.
format(rho_Mean, rho_Q1, rho_Med, rho_Q3), fontsize = 14)
plt.xlabel('$\omega$', fontsize = 16)
plt.show()
Xtot = X.sum(axis = 1)
Fs = []
for i in range(N):
Fs.append(d_fm[f'cF.{i+1}'])
Fs = np.array(Fs).T
Rs = []
R2s = []
for t in range(len(Fs)):
r = np.corrcoef(Xtot,
Fs[t])[0][1]
Rs.append(r)
R2s.append(r**2)
R_Q1, R_med, R_Q3 = np.percentile(Rs,
[25, 50, 75])
R2_Q1, R2_med, R2_Q3 =
np.percentile(R2s, [25, 50, 75])
sb.kdeplot(Rs, label = 'R')
sb.kdeplot(R2s, label = '$R^2$')
plt.title('R(Med.) = {0:.3} $R^2$(Med.) =
{1:.3}'.format(R_med, R2_med), fontsize = 16)
plt.show()
FsMeds = np.median(Fs, axis=0)
plt.plot(FsMeds, Xtot, 'o')
plt.xlabel('F', fontsize = 16)
plt.ylabel('X', fontsize = 16)
plt.show()
fout.write('\n')
print('\nR^2(Med.) =
{0:.5}'.format(R2_med))
fout.write('\nR^2(Med.) =
{0:.5}\n'.format(R2_med))
print('R^2(Q1) = {0:.5} R^2(Q3) = {1:.5}'.format(R2_Q1,
R2_Q3))
fout.write('R^2(Q1) = {0:.5} R^2(Q3) =
{1:.5}\n'.format(R2_Q1, R2_Q3))
print('R(Med.) =
{0:.5}'.format(R_med))
fout.write('R(Med.) =
{0:.5}\n'.format(R_med))
print('R(Q1) = {0:.5} R(Q3) =
{1:.5}'.format(R_Q1, R_Q3))
fout.write('R(Q1) = {0:.5} R(Q3) =
{1:.5}'.format(R_Q1, R_Q3))
S = np.cov(X, rowvar = False)
VarX = np.var(X.sum(axis = 1))
alpha_coef = (M/(M-1)) * (1 -
(np.diag(S).sum() / VarX))
print('\nalpha =
{0:.5}'.format(alpha_coef))
fout.write('\n')
fout.write('\nalpha =
{0:.5}\n'.format(alpha_coef))
Cor = np.corrcoef(X, rowvar = False)
r_bar = (Cor.sum() -
np.diag(Cor).sum()) / (M * (M - 1))
std_alpha = M * r_bar / (1.0 + (M - 1)
* r_bar)
print('standardized alpha =
{0:.5}'.format(std_alpha))
fout.write('standardized alpha =
{0:.5}'.format(std_alpha))
fout.close()
print('Results.txt was saved.')
Listing 3 Python script in the case of K=2. The files are contained in the archived file files8_2Cat.zip , which can be freely downloaded and used.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
import stan
import seaborn as sb
import arviz as az
pd.options.display.max_rows = 10000
fout = open('Results.txt', 'w')
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
data = pd.read_excel(xlsx)
fout.write(f'Data file = {flnm}\n')
VarLabels = data.columns[1:].values;
print(VarLabels)
X = data.values[:, 1:]
print(X[:10])
N = len(X)
M = len(X[0])
print('N = ', N)
print('M = ', M)
fout.write('\n')
fout.write('VarLables
=\n{VarLables}\n')
fout.write(f'N = {N}, M = {M}\n')
fout.write('\nX =\n')
for i, v in enumerate(X):
fout.write(f'{i+1:>5d}
{v}\n')
Data = {'N': N, 'M': M, 'X': X}
def f_init():
return dict(mu = [0.0]*M,
psy = 1.0,
F = np.random.normal(size = N))
with open('Reli_2Cat.stan', 'r') as
fstan:
sm =
stan.build(fstan.read(), data = Data)
fit = sm.sample(num_chains = 4,init =
[f_init()]*4)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data, var_names =
['mu', 'psy', 'cF', 'mnF', 'sdF'])
print('Summary =\n',smry)
fout.write('\nSummary...\n')
fout.write(f'{smry}\n')
d_fm = fit.to_frame()
print('d_fm...keys...\n', d_fm.keys())
print(fit)
Mus = []
for j in range(M):
Mus.append(d_fm[f'mu.{j+1}'])
Mus = np.array(Mus).T
for i in range(M):
sb.kdeplot(Mus.T[i])
plt.title('$\mu$', fontsize = 16)
plt.show()
print()
fout.write('\n')
for j in range(M):
print('mu[{0:}] =
{1:.5}'.format(j+1, np.percentile(Mus.T[j], 50)))
fout.write('mu[{0:}] =
{1:.5}\n'.format(j+1, np.percentile(Mus.T[j], 50)))
Psys = np.array(d_fm['psy'])
sb.kdeplot(Psys)
psy_Q1, psy_Med, psy_Q3 =
np.percentile(Psys, [25, 50, 75])
plt.title(f'Mean =
{np.mean(Psys):.3f}\nQ1 = {psy_Q1:.3f},
Med.= {psy_Med:.3f}' +
f', Q3 = {psy_Q3:.3f}' ,
fontsize = 16)
plt.xlabel('$\psi$', fontsize = 20)
plt.tight_layout()
plt.show()
print()
fout.write('\n')
print('Psys =\n', Psys)
print('Psy =
{0:.5f}\n'.format(np.percentile(Psys, 50)))
fout.write('Psy =
{0:.5f}\n'.format(np.percentile(Psys, 50)))
rho = []
for i in range(len(Psys)):
rho.append(M**2 / (M**2 +
M*(Psys[i]**2)))
rho_Q1, rho_Med, rho_Q3 =
np.percentile(rho, [25, 50, 75])
rho_Mean = np.array(rho).mean()
print('\nomega(mean) =
{0:.5}'.format(rho_Mean))
fout.write('\nomega(mean) =
{0:.5}\n'.format(rho_Mean))
print('omega(Med.) =
{0:.5}'.format(rho_Med))
fout.write('omega(Med.) =
{0:.5}\n'.format(rho_Med))
print('omega(Q1) = {0:.5} omega(Q3) =
{1:.5}'.format(rho_Q1, rho_Q3))
fout.write('omega(Q1) = {0:.5} omega(Q3) =
{1:.5}\n'.format(rho_Q1, rho_Q3))
sb.kdeplot(rho)
plt.title('Mean = {0:.3}\nQ1 = {1:.3},
Med.= {2:.3}, Q3 = {3:.3}'.
format(rho_Mean, rho_Q1, rho_Med, rho_Q3), fontsize = 14)
plt.xlabel('$\omega$', fontsize = 16)
plt.show()
Xtot = X.sum(axis = 1)
Fs = []
for i in range(N):
Fs.append(d_fm[f'cF.{i+1}'])
Fs = np.array(Fs).T
Rs = []
R2s = []
for t in range(len(Fs)):
r = np.corrcoef(Xtot,
Fs[t])[0][1]
Rs.append(r)
R2s.append(r**2)
R_Q1, R_med, R_Q3 = np.percentile(Rs,
[25, 50, 75])
R2_Q1, R2_med, R2_Q3 =
np.percentile(R2s, [25, 50, 75])
sb.kdeplot(Rs, label = 'R')
sb.kdeplot(R2s, label = '$R^2$')
plt.title('R(Med.) = {0:.3} $R^2$(Med.) =
{1:.3}'.format(R_med, R2_med), fontsize = 16)
plt.legend(fontsize = 16)
plt.show()
FsMeds = np.median(Fs, axis=0)
plt.plot(FsMeds, Xtot, 'o')
plt.xlabel('F', fontsize = 16)
plt.ylabel('X', fontsize = 16)
plt.show()
fout.write('\n')
print('\nR^2(Med.) =
{0:.5}'.format(R2_med))
fout.write('\nR^2(Med.) =
{0:.5}\n'.format(R2_med))
print('R^2(Q1) = {0:.5} R^2(Q3) =
{1:.5}'.format(R2_Q1, R2_Q3))
fout.write('R^2(Q1) = {0:.5} R^2(Q3) = {1:.5}\n'.format(R2_Q1,
R2_Q3))
print('R(Med.) =
{0:.5}'.format(R_med))
fout.write('R(Med.) =
{0:.5}\n'.format(R_med))
print('R(Q1) = {0:.5} R(Q3) =
{1:.5}'.format(R_Q1, R_Q3))
fout.write('R(Q1) = {0:.5} R(Q3) =
{1:.5}'.format(R_Q1, R_Q3))
S = np.cov(X, rowvar = False)
VarX = np.var(X.sum(axis = 1))
alpha_coef = (M/(M-1)) * (1 -
(np.diag(S).sum() / VarX))
print('\nalpha =
{0:.5}'.format(alpha_coef))
fout.write('\n')
fout.write('\nalpha =
{0:.5}\n'.format(alpha_coef))
Cor = np.corrcoef(X, rowvar = False)
r_bar = (Cor.sum() -
np.diag(Cor).sum()) / (M * (M - 1))
std_alpha = M * r_bar / (1.0 + (M - 1)
* r_bar)
print('standardized alpha =
{0:.5}'.format(std_alpha))
fout.write('standardized alpha =
{0:.5}'.format(std_alpha))
fout.close()
print('Results.txt was saved.')
McDonald, R. P. (1999). Test theory: A unified treatment. London: Lawrence Erlbaum Associates, Publishers.
Okamoto, Y.
(2013). A Direct Bayesian Estimation of Reliability. Behaviormetrika, 2013, 40, 149-168.