Reliability Coefficient omega 
of an Ordinal Categorical Items Scale
Direct Bayesian Estimation
A reliability coefficient, McDonald’s omega (McDonald, 1999), is based on the following model (1).
denotes the value of
person 
on item 
.
denotes a factor score
of person . It is assumed to have the standard normal distribution.
The residual 
has the following
distribution
The score of person is given by the sum of
categorical values as follows
Based on this model, the reliability coefficient McDonald’s omeg (McDonald, 1999) is given by eq. (2)
In the case of an ordinal categorical items
scale, a categorical value of an item is given based on the continuous value of
.
Number of categorical values is denoted by 
, and a categorical value is represented by an integer from 1 to .
Let’s set the following model (Okamoto, 2013).
Category boundaries is set as follows.
Using these boundaries, a categorical
response 
is given based on as follows.
Because
we have
denotes the cumulative
function of the standard normal distribution, and it is assumed that
To estimate parameter values, put the following constraint
Furthermore, here, it is assumed that the number of categories is larger than 2
For the reason, see Okamoto (2013).
Now, we have the following likelihood function
Hence, we have the following posterior distribution
The prior distribution 
is set to be almost flat
(Gelman et al., 2021).
A Stan script for the above model is shown in Listing 1.
A Python script, which uses the script in Listing 2, is shown in Listing 2.
Script files and sample data files are archived into the ZIP file omega_files.zip, which can be freely downloaded and used.
The Python script in Listing 2 was run on Python 3.7 virtual environment on Anaconda on Windows 10.
For information on how to install Python and PyStan (the Python interface to Stan), check this website for Windows, or this website for Linux. This website is for Windows and for Linux.
The script in Listing 2 can be executed by the following code
(py37) PS
D:\omega_files> python CalcOmegaMain.py
CalcOmegaMain.py is the name of the script in Listing 2.
When the script starts, the name of the input data file was asked as follows
(py37) PS
D:\omega_files> python CalcOmegaMain.py
Input file =
An input data file should be prepared as an Excel file like one shown in Figure 1.
Figure 1
In the first row, the number of categories, which are the same numbers, is set.
In the second row, names of items are set.
In the third row, integer values 1 or 0 are set. Integer 1 means that this item should be included in the analysis, 0 means that this item is excluded from the analysis.
Data of the category responses are set from the fourth row on.
In the first column, case identification strings or numbers are set. In the analysis, these values are not used.
From the second column, category responses of a person are set in each row.
When the name of the input data file is set like the following, the program starts.
(py37) PS D:\omega_files> python
CalcOmegaMain.py
Input file = DataCat5.xlsx
After MCMC sampling ends, statistics of the sampling are shown as follows.
Elapsed Time: 75.129 seconds (Warm-up)
31.615 seconds (Sampling)
106.744 seconds (Total)
Inference for Stan model:
anon_model_4de137d01c5c96aab1aeacdbb8a86c8c.
4 chains, each with iter=2000;
warmup=1000; thin=1;
post-warmup draws per chain=1000,
total post-warmup draws=4000.
mean se_mean
sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu[0] 0.53 7.5e-4 0.04 0.45 0.5 0.53 0.56 0.6 2571 1.0
mu[1] 0.48 7.5e-4 0.04 0.41 0.46 0.48 0.5 0.55 2236 1.0
mu[2] 0.51 6.1e-4 0.04 0.43 0.48 0.51 0.53 0.58 4000 1.0
mu[3] 0.5 7.9e-4 0.04 0.42 0.47 0.5 0.52 0.57 2299 1.0
mu[4] 0.5 7.7e-4 0.04 0.42 0.47 0.5 0.53 0.58 2566 1.0
lmbd[0] 0.33 8.5e-4 0.04 0.25 0.3 0.33 0.36 0.41 2608 1.0
.
.
.
Then, posterior distribution of the omega 
is shown (Figure 2).
Figure 2
Close the window of Figure 2, posterior
distribution of correlation coefficient between and 
is shown (Figure 3).
Figure 3
Close the window of Figure 3, execution of the script ends.
In the terminal, messages like the following are displayed.
FigOmega.png was saved.
FigCors,png was saved.
Output.txt was saved.
(py37) PS D:\omega_files>
Files FigOmega.png, FigCors,png contain Figure 2 and 3,
respectively.
The file
Output.txt contains the following outputs.
.
.
.
Lambda...
Item-1: 0.299 0.328 0.357
Item-2: 0.290 0.317 0.347
Item-3: 0.276 0.306 0.337
Item-4: 0.296 0.324 0.353
Item-5: 0.295 0.325 0.356
.
.
.
For each item, the first quantile value Q1, the median, the third quantile value Q3 are shown in each row.
The values of Lambda of Item-3 are the smallest.
In data file of Figure 4, integer 0 in set for Item-3 to exclude from the analysis.
Figure 4
Then, run the script as follows.
(py37) PS D:\omega_files> python CalcOmegaMain.py
Input file = DataCat5del3.xlsx
File name DataCat5del3.xlsx is the name of the file shown in Figure 4.
Then, we get the results as shown in Figure 5.
Figure 5
Values, which are a little smaller than the values in Figures 2 and 3, are shown.
Maybe, effect of excluding the item with smaller lambda is cancelled by that of reducing the number of items.
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and Other Stories. Cambridge University Press.
McDonald, R. P. (1999). Test Theory: A unified treatment. Lawrence Erlbaum Associations, Publishers.
Okamoto,
Y. (2013). “A Direct Bayesian Estimation of
Reliability” Behaviormetrika,
40, 149-168.
Listing 1. Stan script for the model (calc_omega.stan)
//
//
Yasuharu Okamoto, 2023.02.01
//
data {
int n_items;
int n_persons;
int<lower = 3> K;
// Number of
categories is larger than 2
int X[n_persons,n_items];
}
transformed data {
vector[K-2] a_K;
for (k in 1:(K-2)) {
a_K[k] = 1.0;
}
}
parameters {
real F[n_persons];
real mu[n_items];
real<lower = 0.0>
lmbd[n_items];
real<lower = 0.0>
sgm[n_items];
simplex[K-2] theta;
}
transformed parameters {
real C[K-1];
simplex[K]
p[n_persons,n_items];
C[1] = 0.0;
C[K-1] = 1.0;
if (K > 3) {
for
(j in 2:(K-2)){
C[j] = C[j-1] + theta[j-1];
}
}
for (i in 1:n_persons) {
for
(j in 1:n_items) {
p[i][j][1] = normal_cdf(C[1],
mu[j] + lmbd[j] * F[i], sgm[j]);
p[i][j][K] = 1.0 - normal_cdf(C[K-1],
mu[j] + lmbd[j] * F[i], sgm[j]);
for (k in 2:(K-1)) {
p[i][j][k] = normal_cdf(C[k],
mu[j] + lmbd[j] * F[i], sgm[j]) -
normal_cdf(C[k-1],
mu[j] + lmbd[j] * F[i], sgm[j]);
}
}
}
}
model {
theta ~ dirichlet(a_K);
for (i in 1:n_persons) {
F[i]
~ normal(0.0, 1.0);
}
for (j in 1:n_items) {
lmbd[j] ~ exponential(0.1);
mu[j] ~ normal(0.0, 10.0);
sgm[j] ~ exponential(0.1);
}
for (i in 1:n_persons) {
for
(j in 1:n_items) {
X[i][j] ~ categorical(p[i][j]);
}
}
}
Listing 2. Python script file (CalcOmegaMain.py)
import pandas as pd
import numpy as np
import pystan
import platform
import pickle
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
in_flnm = input('Input file(*.xlsx) =
')
# Excel file
out_flnm = 'Output.txt'
out_f = open(out_flnm, 'w')
out_f.write('Input data file = ' +
in_flnm + '\n')
#
# Read in the Excel file
#
fl_xlsx = pd.ExcelFile(in_flnm)
data = pd.read_excel(fl_xlsx, header =
None)
data = data.values
K = data[0][1]
print('K =', K)
print(data[:10])
out_f.write('\ndata...\n')
for v in data:
out_f.write(v.__str__() +
'\n')
item_names = data[1]
print(item_names)
ck_items = data[2]
sel_items = ck_items == 1
print(sel_items)
item_names = item_names[sel_items]
print('Selected items...')
print(item_names)
out_f.write('\n\nThe selected item
names...' + item_names.__str__() + '\n')
data = np.array(data)
X = data[3:].T[1:].T
X = X[:, sel_items[1:]]
print(X[:3])
#
# MCMC
sampling by Stan
#
sm = pystan.StanModel(file =
'calc_omega.stan')
X = X.astype(dtype = int)
fit = sm.sampling(data =
{'n_items':len(X[0]), 'n_persons':len(X), 'K':K, 'X': X},
pars = ['mu', 'lmbd', 'sgm', 'F'],
n_jobs = 1 if platform.system() == 'Windows' else -1)
#
#
Results of MCMC sampling
#
print(fit)
out_f.write('\nfit...\n')
out_f.write(fit.__str__())
#
# Calculation on
parameter mu
#
print('mu =\n', fit['mu'[:10]])
out_f.write('\nMu (Q1, Med, Q3)...\n')
for j in range(len(item_names)):
out_f.write(' ' + item_names[j] + ': ')
mu_q1, mu_med, mu_q3 = np.percentile(fit['mu'].T[j],
[25, 50, 75])
out_f.write('{0:.3f} {1:.3f} {2:.3f}\n'.format(mu_q1, mu_med,
mu_q3))
#
#
Calculation on parameter lambda
#
print('lmbd =\n', fit['lmbd'][:10])
print('Lambda =\n', fit['lmbd'[:10]])
out_f.write('\nLambda (Q1, Med,
Q3)...\n')
for j in range(len(item_names)):
out_f.write(' ' + item_names[j] + ': ')
lmbd_q1, lmbd_med, lmbd_q3 =
np.percentile(fit['lmbd'].T[j], [25, 50, 75])
out_f.write('{0:.3f} {1:.3f} {2:.3f}\n'.format(lmbd_q1, lmbd_med,
lmbd_q3))
#
#
Calculation on parameter sigma
#
print('Sigma =\n', fit['sgm'[:10]])
out_f.write('\nSigma (Q1, Med,
Q3)...\n')
for j in range(len(item_names)):
out_f.write(' ' + item_names[j] + ': ')
sgm_q1, sgm_med, sgm_q3 =
np.percentile(fit['sgm'].T[j], [25, 50, 75])
out_f.write('{0:.3f} {1:.3f} {2:.3f}\n'.format(sgm_q1, sgm_med,
sgm_q3))
#
#
Calculation on omega
#
sum_lmbd = fit['lmbd'].sum(axis=1)
sum_sqr_sgm =
(fit['sgm']**2).sum(axis=1)
omegas = sum_lmbd**2 / (sum_lmbd**2 +
sum_sqr_sgm)
def calcMAPEst(samples, a = 0.05,
n_points = 1000):
"""
Estimation of the MAP estimate by kde
samples: Input
data for which MAP estimate is calculated.
a:
Each side of a/2 of samples is ignored, Default value is 0.05
n_points: number of
points used in kde. Default value
is 1000
"""
import numpy as np
import scipy.stats as ss
Lp, Up =
np.percentile(samples, [100 * a/2, 100 * (1 - a/2)])
coord = np.linspace(Lp, Up,
n_points)
est_pdf =
ss.gaussian_kde(samples).pdf(coord)
map_idx = np.argmax(est_pdf)
MAP_Est =
coord[map_idx]
return MAP_Est,
est_pdf[map_idx] # MAP estimate, the value at MAP estimate
MAP = calcMAPEst(omegas)[0]
print('MAP =', MAP)
Q1, Med, Q3 = np.percentile(omegas,
[25, 50, 75])
print(Q1, Med, Q3)
out_f.write('\nomega (MAP / Q1, Med,
Q3)...\n')
out_f.write(' {0:.3f}\n'.format(MAP))
out_f.write(' {0:.3f} {1:.3f} {2:.3f}\n'.format(Q1, Med, Q3))
sb.kdeplot(omegas, color = 'b')
plt.title('MAP = {0:.2f}\nQ1 =
{1:.2f} Med = {2:.2f} Q3 = {3:.2f}'.format(MAP, Q1, Med, Q3))
plt.xlabel('$\omega$', fontsize = 14)
plt.savefig('FigOmega.png')
plt.show()
#
#
Calculation on Pearson correlation coefficient
#
between factor F and the sum score of Xjs
#
SumX = X.sum(axis=1)
Cors = []
for Fvalues in fit['F']:
r = ss.pearsonr(SumX,
Fvalues)[0]
Cors.append(r)
r_map = calcMAPEst(Cors)[0]
r_q1, r_med, r_q3 =
np.percentile(Cors, [25, 50, 75])
out_f.write('\nCorrelation between F
and X (MAP / Q1, Med, Q3)...\n')
out_f.write(' {0:.3f}\n'.format(r_map))
out_f.write(' {0:.3f} {1:.3f} {2:.3f}\n'.format(r_q1, r_med,
r_q3))
sb.kdeplot(Cors)
plt.title('MAP = {0:.2f}\nQ1 =
{1:.2f} Med = {2:.2f} Q3 = {3:.2f}'.
format(r_map, r_q1, r_med, r_q3), fontsize = 14)
plt.xlabel('r', fontsize = 16)
plt.savefig('FigCors.png')
plt.show()
print('FigOmega.png was saved.')
print('FigCors,png was saved.')
out_f.close()
print(out_flnm + ' was saved.')