Reliability
Coefficients of Categorical Item Scales
Omega, alpha, and
coefficient of determination
Yasuharu Okamoto,
2020.05
Reliability coefficient of a test score is given by the precision of the test score measurement (McDonald, 1999, p. 63), which is calculated by the following equaiton
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 a standard normal one 
,
and 
is a residual whose distribution is
normal. 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 parameters and 
,
reliability coefficient 
is calculated by eq. (1).
To
investigate relation between the observed score
and
the true latent value
,
correlation
coefficient R and coefficient of determination 
are calculated (Okamoto,
2016, 2017).
Reliability coefficient 
is also calculated.
Actual
installation of a model is depend on the number of categories K. 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 2. Scripts in PyStan 3 are shown in this website.
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;
int<lower = 1, upper =
K> X[N, M];
vector[K-2] a;
}
transformed data {
real C1;
real CK_1;
C1 = -1;
CK_1 = 1;
}
parameters {
real mu[M];
real<lower = 0.0>
Lmd[M];
real<lower = 0.0>
psy[M];
simplex[K-2] Ck;
real F[N];
}
transformed parameters {
vector[K] p[N, M];
vector[K-1] cumP[N, M];
real C[K-1];
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, 100.0);
psy[j] ~ exponential(0.01);
Lmd[j] ~ exponential(0.01);
}
Ck ~ dirichlet(a);
for (i in 1:N) {
F[i]
~ normal(0.0, 1.0);
}
for (i in 1:N) {
for
(j in 1:M) {
X[i][j] ~
categorical(p[i][j]);
}
}
}
The Python script, which uses this Stan script to analyze the data, is shown in Listing 1 in the latter half of this website.
Run the script, then the name of the input data file is asked (Figure 1.1).
Figure 1.1
The input data file is expected to be an Excel data file (*.xlsx) (Figure 1.2).
.
.
.
Figure 1.2
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 recorded. Values on the items should be integer values from 1 to K, where K is the number of categories. The data in Figure 1.2 is for the case of K=5.
After the input data file name is set, the number of response categories K is asked (Figure 1.3).
Figure 1.3
After the value of K is set, sampling by
Stan starts. When the sampling by Stan ends, posterior distributions of 
s are displayed (Figure 1.4).
Figure 1.4
Click the X icon at the upper right corner
of the windows to close. The next window shows posterior distributions of 
s (Figure 1.5).
Figure 1.5
Close the window in Figure 1.5, then the
window in Figure 1.6 appears and the posterior distributions of 
s are presented (Figure 1.6).
Figure 1.6
Close the window in Figure 1.6, then the
next window appears and the posterior distribution of is presented (Figure
1,7).
Figure 1.7
Close the window in Figure 1.7, then the
window, where the posterior distributions of
R and are drawn, appears
(Figure 1.8). 
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.8
When the window in Figure 1.8 is closed, the program ends. After the program ends, the terminal display shows the following messages.
mu[1] = -0.59881
mu[2] = -0.52904
mu[3] = -0.36531
mu[4] = -0.25324
mu[5] = -0.13544
mu[6] = 0.013651
mu[7] = 0.10306
mu[8] = 0.30145
mu[9] = 0.49264
mu[10] = 0.46195
Lambda[1] = 0.63104
Lambda[2] = 0.70976
Lambda[3] = 0.71212
Lambda[4] = 0.81548
Lambda[5] = 0.7646
Lambda[6] = 0.79926
Lambda[7] = 0.69422
Lambda[8] = 0.83381
Lambda[9] = 0.79007
Lambda[10] = 0.7185
Psy[1] = 0.30787
Psy[2] = 0.3421
Psy[3] = 0.32802
Psy[4] = 0.35826
Psy[5] = 0.31646
Psy[6] = 0.34562
Psy[7] = 0.31062
Psy[8] = 0.43288
Psy[9] = 0.33922
Psy[10] = 0.27833
omega(mean) = 0.97924
omega(Med.) = 0.97947
omega(Q1) = 0.97702 omega(Q3) = 0.98171
R^2(Med.) = 0.96031
R^2(Q1) = 0.95617 R^2(Q3) = 0.96423
R(Med.) = 0.97995
R(Q1) = 0.97784 R(Q3) = 0.98195
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;
int<lower = 1, upper =
3> X[N, M];
}
transformed data {
real C1;
real C2;
C1 = -1;
C2 = 1;
}
parameters {
real mu[M];
real<lower = 0.0>
Lmd[M];
real<lower = 0.0>
psy[M];
real F[N];
}
transformed parameters {
vector[3] p[N, M];
real cumP[N, M];
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, 100.0);
psy[j] ~ exponential(0.01);
Lmd[j] ~ exponential(0.01);
}
for (i in 1:N) {
F[i]
~ normal(0.0, 1.0);
}
for (i in 1:N) {
for
(j in 1:M) {
X[i][j] ~ categorical(p[i][j]);
}
}
}
The Python script, which uses this Stan script to analyze the data, is shown in Listing 2 in the latter half of this website.
Run the script, then the name of the input data file is asked (Figure 2.1).
Figure 2.1
The input data file is expected to be an Excel data file (*.xlsx) (Figure 2.2).
・
・
・
Figure 2.2
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 recorded. Values on the items should be integer values
from 1 to K, where K is the number of categories. In this case of K=3, 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.3).
Figure 2.3
Click the X icon at the upper right corner
of the windows to close. The next window shows posterior distributions of s are displayed (Figure 2.4).
Figure 2.4
Close the window in Figure 2.4, then the
window in Figure 2.5 appears and the posterior distributions of s are presented.
Figure 2.5
Close the window in Figure 2.5, then the
next window appears and the posterior distribution of is presented (Figure
2.6).
Figure 2.6
Close the window in Figure 2.6, then the
window, where the posterior distributions of
R and are drawn, appears
(Figure 2.7). 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.7
When the window in Figure 2.7 is closed, the program ends. After the program ends, the terminal display shows the following messages.
mu[1] = -1.7334
mu[2] = -1.535
mu[3] = -0.59251
mu[4] = -0.4819
mu[5] = -0.058653
mu[6] = 0.23374
mu[7] = 0.81563
mu[8] = 1.219
mu[9] = 2.6697
mu[10] = 2.2209
Lambda[1] = 1.9298
Lambda[2] = 2.2258
Lambda[3] = 1.9386
Lambda[4] = 2.2223
Lambda[5] = 2.393
Lambda[6] = 2.5833
Lambda[7] = 2.676
Lambda[8] = 2.7033
Lambda[9] = 3.9394
Lambda[10] = 2.8044
Psy[1] = 0.86738
Psy[2] = 1.3084
Psy[3] = 1.086
Psy[4] = 1.2532
Psy[5] = 0.95873
Psy[6] = 0.7304
Psy[7] = 0.89208
Psy[8] = 1.1605
Psy[9] = 1.6746
Psy[10] = 1.1748
omega(mean) = 0.97897
omega(Med.) = 0.97937
omega(Q1) = 0.97624 omega(Q3) = 0.9822
R^2(Med.) = 0.89283
R^2(Q1) = 0.87786 R^2(Q3) = 0.90551
R(Med.) = 0.9449
R(Q1) = 0.93694 R(Q3) = 0.95158
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;
int<lower = 1, upper =
2> X[N, M];
}
parameters {
real mu[M];
real<lower = 0.0> psy;
real F[N];
}
transformed parameters {
real p[N, M];
for (i in 1:N)
for
(j in 1:M){
p[i][j] = normal_cdf((mu[j] + F[i]) / psy, 0.0, 1.0);
}
}
model {
psy ~ exponential(0.01);
for (j in 1:M) {
mu[j] ~ normal(0.0, 100.0);
}
for (i in 1:N) {
F[i]
~ 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.
Run the script, then the name of the input data file is asked (Figure 3.1).
Figure 3.1
The input data file is expected to be an Excel data file (*.xlsx) (Figure 3.2).
・
・
・
Figure 3.2
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 recorded. 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.3).
Figure 3.3
Click the X icon at the upper right corner
of the windows to close. The window in Figure 3.4 appears and the posterior distributions
of is presented (Figure
3.4).
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
When the window in Figure 3.5 is closed, the program ends. After the program ends, the terminal display shows the following messages.
mu[1] = -1.2527
mu[2] = -1.1176
mu[3] = -0.77675
mu[4] = -0.29035
mu[5] = -0.26226
mu[6] = 0.025895
mu[7] = 0.15128
mu[8] = 0.33489
mu[9] = 0.51242
mu[10] = 0.88316
omega(mean) = 0.95671
omega(Med.) = 0.95745
omega(Q1) = 0.95134 omega(Q3) = 0.96292
R^2(Med.) = 0.8555
R^2(Q1) = 0.83992 R^2(Q3) = 0.86943
R(Med.) = 0.92493
R(Q1) = 0.91647 R(Q3) = 0.93243
alpha = 0.88414
standardized alpha = 0.88536
Listing 1 Python script in the case of K>3. 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 pystan
import seaborn as sb
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
data = pd.read_excel(xlsx)
VarLabels = data.columns[1:].values;
print(VarLabels)
X = data.values[:, 1:]
N = len(X)
M = len(X[0])
print('N = ', N)
print('M = ', M)
K = int(input('Number of categories =
'))
print('K = ', K)
Data = {'N': N, 'M': M, 'X': X, 'K':
K, 'a': [1.0]*(K-2)}
sm = pystan.StanModel(file =
'Reli_G3Cat.stan')
fit = sm.sampling(data = Data, pars =
['mu', 'Lmd', 'psy', 'Ck', 'F'], n_jobs = 1)
print(fit)
Mus = fit['mu']
for i in range(M):
sb.kdeplot(Mus.T[i])
plt.title('$\mu$', fontsize = 16)
plt.show()
print()
for j in range(M):
print('mu[{0:}] =
{1:.5}'.format(j+1, np.percentile(Mus.T[j], 50)))
Lmds = fit['Lmd']
for i in range(M):
sb.kdeplot(Lmds.T[i])
plt.title('$\lambda$', fontsize = 24)
plt.show()
print()
for j in range(M):
print('Lambda[{0:}] =
{1:.5}'.format(j+1, np.percentile(Lmds.T[j], 50)))
Psys = fit['psy']
for j in range(M):
sb.kdeplot(Psys.T[j])
#plt.hist(Psys)
plt.title('$\psi$', fontsize = 24)
plt.show()
print()
for j in range(M):
print('Psy[{0:}] =
{1:.5}'.format(j+1, np.percentile(Psys.T[j], 50)))
Fs = fit['F']
rho = []
for i in range(len(Lmds)):
#rho.append(M**2 / (M**2 +
M*(Psys[i]**2)))
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()
print('\nomega(mean) =
{0:.5}'.format(rho_Mean))
print('omega(Med.) =
{0:.5}'.format(rho_Med))
print('omega(Q1) = {0:.5} omega(Q3) =
{1:.5}'.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)
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()
print('\nR^2(Med.) =
{0:.5}'.format(R2_med))
print('R^2(Q1) = {0:.5} R^2(Q3) =
{1:.5}'.format(R2_Q1, R2_Q3))
print('R(Med.) =
{0:.5}'.format(R_med))
print('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))
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))
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 pystan
import seaborn as sb
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
data = pd.read_excel(xlsx)
VarLabels = data.columns[1:].values;
print(VarLabels)
X = data.values[:, 1:]
N = len(X)
M = len(X[0])
print('N = ', N)
print('M = ', M)
Data = {'N': N, 'M': M, 'X': X}
sm = pystan.StanModel(file =
'Reli_3Cat.stan')
fit = sm.sampling(data = Data, pars =
['mu', 'Lmd', 'psy', 'F'], n_jobs = 1)
print(fit)
Mus = fit['mu']
for i in range(M):
sb.kdeplot(Mus.T[i])
plt.title('$\mu$', fontsize = 16)
plt.show()
print()
for j in range(M):
print('mu[{0:}] =
{1:.5}'.format(j+1, np.percentile(Mus.T[j], 50)))
Lmds = fit['Lmd']
for i in range(M):
sb.kdeplot(Lmds.T[i],)
plt.title('$\lambda$', fontsize = 24)
plt.show()
print()
for j in range(M):
print('Lambda[{0:}] =
{1:.5}'.format(j+1, np.percentile(Lmds.T[j], 50)))
Psys = fit['psy']
for j in range(M):
sb.kdeplot(Psys.T[j])
#plt.hist(Psys)
plt.title('$\psi$', fontsize = 24)
plt.show()
print()
for j in range(M):
print('Psy[{0:}] =
{1:.5}'.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()
print('\nomega(mean) =
{0:.5}'.format(rho_Mean))
print('omega(Med.) =
{0:.5}'.format(rho_Med))
print('omega(Q1) = {0:.5} omega(Q3) =
{1:.5}'.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 = fit['F']
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()
print('\nR^2(Med.) =
{0:.5}'.format(R2_med))
print('R^2(Q1) = {0:.5} R^2(Q3) =
{1:.5}'.format(R2_Q1, R2_Q3))
print('R(Med.) =
{0:.5}'.format(R_med))
print('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))
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))
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 pystan
import seaborn as sb
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
data = pd.read_excel(xlsx)
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)
Data = {'N': N, 'M': M, 'X': X}
def f_init():
return dict(mu = [0.0]*M,
psy = 1.0, F = [0.0]*N)
sm = pystan.StanModel(file = 'Reli_2Cat.stan')
fit = sm.sampling(data = Data, pars =
['mu', 'psy', 'F'], init = f_init, n_jobs = 1)
print(fit)
Mus = fit['mu']
for i in range(M):
sb.kdeplot(Mus.T[i])
plt.title('$\mu$', fontsize = 16)
plt.show()
print()
for j in range(M):
print('mu[{0:}] =
{1:.5}'.format(j+1, np.percentile(Mus.T[j], 50)))
Psys = fit['psy']
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))
print('omega(Med.) =
{0:.5}'.format(rho_Med))
print('omega(Q1) = {0:.5} omega(Q3) =
{1:.5}'.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 = fit['F']
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()
print('\nR^2(Med.) =
{0:.5}'.format(R2_med))
print('R^2(Q1) = {0:.5} R^2(Q3) =
{1:.5}'.format(R2_Q1, R2_Q3))
print('R(Med.) =
{0:.5}'.format(R_med))
print('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))
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))
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.