順序カテゴリ項目尺度の信頼性係数ω（オメガ）
ベイズ分析による直接推定
順序カテゴリ項目から構成される尺度の信頼性係数を求める。各項目は、K個のカテゴリ値をとるものとし、１からKまでの整数値でカテゴリ値を表す。
項目
の回答者
の潜在変数
を考え、この値によってカテゴリ値が決まるとする。
McDonald (1999)に従い、１因子モデル
を設定する。
、
は項目に対する定数、
は回答者の因子得点であり、
であるとする。
は、いわゆる誤差項であり、と独立な確率変数で、
とする。このとき、テスト得点の和
の信頼性係数
は次式（McDonald (1999)の
）で与えられる（岡本、2014）。
カテゴリ境界
を
と設定し、の値ととの関係より、回答選択肢（カテゴリ値）
が次式
により決まるものとする。このとき、、、、
、の値に対して
であり、
となる。ここで、
は累積標準正規分布関数で、
とおく。また、の次元の単位と原点を次式の制約で決める。
ここで、
という制約を置いている。すなわち、カテゴリ数が３以上の場合について考える。
以上の仮定の下で、データに対する尤度関数は次式で与えられる。
したがって、事後確率は
で与えられる。
は事前確率分布で、十分に平らな分布を用いる（Gelman et al., 2021）。
上記の確率モデルのStanスクリプトをリスト１に示す。
リスト１のStanスクリプトを用いて、上記の確率モデルのベイズ分析を行い、式（２）のを求めるPythonスクリプトをリスト２に示す。
リスト１とリスト２のスクリプトファイル、およびサンプルデータファイルは、圧縮ファイルomega_files.zipとしてまとめた。ダウンロード展開して、自由に利用されたい。
PythonによるStanの利用は、岡本（2019）で説明している。
Pythonの準備については、このウェブサイトで説明している。
以下のリスト２の実行例では、Anaconda/Windows 10上のPython 3.7仮想環境においてPyStan 2.17.1.0を利用した。
リスト２のスクリプトは次のコマンドで実行できる。
(py37) PS
D:\omega_files> python CalcOmegaMain.py
上のコマンドで実行すると、入力データファイルの設定が求められる。
(py37) PS
D:\omega_files> python CalcOmegaMain.py
Input file =
入力データファイルは、図１に示す形式でExcelファイルとして用意する。
図１
各項目の回答値（カテゴリ値）かは、２列目から設定する。１列目は、各ケースの名前で、数値でも文字列でもよい。
１行目は、1列目に適当な文字列を設定し、２列目からは各項目のカテゴリ数Kを設定する。
２行目は、各列の変数名である。
３行目は、項目の選択・非選択の指示であり、１が設定されている項目が分析の対象になる。図１では、すべて１が設定されているので、すべての項目が取り込まれて分析されてる。
４行目から、各ケースのデータを１行ずつ設定していく。
以下のように、入力データファイル名を設定して、Enterキーを押すと、計算が始まる。
下の例では、図１のExcelファイル名が設定されている。
(py37) PS D:\omega_files> python
CalcOmegaMain.py
Input file = DataCat5.xlsx
MCMCサンプリングが終了すると、以下のように統計量が表示される。
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
・
・
・
その後、ωの事後分布が表示される（図２）。
図２
図２のWindowを閉じると、因子と項目和
との相関係数ｒの事後分布のグラフが表示される（図３）。
図３
図３のWindowを閉じると、プログラムの実行終了となる。
端末画面には、以下のように出力されている。
FigOmega.png was saved.
FigCors,png was saved.
Output.txt was saved.
(py37) PS D:\omega_files>
ファイルFigOmega.pngに図２のグラフ、ファイルFigCors,pngに図３のグラフが保存されている。
ファイルOutput.txtには、計算結果の出力が以下のように書き出されている。
・
・
・
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
・
・
・
各項目のLambdaの事後分布の第一四分位数Q1、中央値Median、第三四分位数Q3が出力されている。
項目Item-3のLambdaの値が、他のものより小さいので、この項目を除いて分析するためのファイルを図４のように用意してみる。
図４
３行目４列目のItem-3に対する選択支持の値が０に設定されていることに注意。
このファイルを入力データファイルとして、以下のようにプログラムを実行する。
(py37) PS D:\omega_files> python
CalcOmegaMain.py
Input file = DataCat5del3.xlsx
実行結果は、図５のようになる。
図５
項目を除いた結果、ωとｒの値は少し下がったようである。これは、Lambdaの小さい項を除いて期待される＋の効果よりも、項目数が減ることのーの効果の方が大きかったようである。
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.
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版
Okamoto,
Y. (2013). “A Direct Bayesian Estimation of
Reliability” Behaviormetrika,
40, 149-168.
リスト１ Stanスクリプト（calc_omega.stan）
data {
int n_items;
int n_persons;
int<lower = 3> K;
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]);
}
}
}
リスト２ Pythonスクリプト（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.')