カテゴリ項目尺度の信頼性係数
信頼性係数
は次式（１）で与えられる(McDonald, 1999)
式（１）は、式（２）のモデルに基づくものである。
ここで、
は回答者
の項目
における回答カテゴリであり(
)、
はの平均値、
は項目の因子負荷量、
は回答者の因子得点で標準正規分布に従うものとし、
は残差（誤差）である。
いま、連続変量
が、カテゴリ境界
によって、以下のように離散化されてカテゴリ回答
が与えられるものとする(Okamoto, 2013)。
ここで、
であり、Kはカテゴリ回答のカテゴリ数である。
このとき、カテゴリ回答
に対して次式を設定する。
ここで、
は累積標準正規分布関数であり、
、 
.
とする。
このとき、尤度関数はデータに対して次式で与えられる。
したがって、事後確率は次式で与えられる。
ここで、
.は事前分布である。
パラメータと
から、式（１）によりが算出される。
本サイトのプログラムでは、測定値としての各項目の回答カテゴリ値の和
と回答者の特性の真値としての因子得点
,
の相関係数
c、およびその２乗である決定係数
も算出されている(Okamoto, 2016, 2017).
また、アルファ係数
も算出されているので、参考になるであろう。.
上のモデルでは、Kが２の場合、Kが３の場合、Kが３より多い場合に応じてプログラムの実装が異なるので、それぞれの場合に分けて扱う。
なお、PythonにおけるStan (PyStan)の使い方については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第10章などを参照されたい。
K＞３のときのStanスクリプト（ファイル名Reli_G3Cat.stan）を以下のように用意した。
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]);
}
}
}
両端のカテゴリ境界は
に固定されている（Okamoto, 2013）。
上のStanスクリプトを利用して信頼性係数
を求めるスクリプトはこのウェブサイトの後ろにリスト１として示した。
リスト１のプログラムを実行すると、入力データファイル名が求められる（図1.1）。
図1.1
入力データファイルはExcelファイルとして用意する（図1.2）。
・
・
・
図1.2
1行目に項目変数名を書く。第１列はデータの識別用変数であるが、これはデータの識別用で、通し番号でも文字列でもなんでもよい。第２列目から項目反応カテゴリ値のデータをかく。カテゴリ値は、カテゴリ数がKの場合、１からKまでの整数値で表す。また、逆転項目があれば、値を逆にして他の項目とカテゴリの方向を揃える。カテゴリ値は、測定概念の値と同順で、測定概念の特性が強いほどカテゴリ値も大きいとする。逆転項目のカテゴリ値を逆にするプログラムは、別のウェブサイトに用意した。
入力データファイル名が設定されてデータが読み込まれると、次にカテゴリ値の数Kの入力が求められる（図1.3）。
図1.3
図1.3では、データに合わせて５が入力されている。カテゴリ数が入力されるとStanによるサンプリングが始まる。サンプリングが終了すると、まず
の事後分布が表示される（図1.4）。
図1.4
図1.4のウィンドウの右上のX印アイコンをクリックして閉じると、次に
の事後分布が描かれる（図1.5）。
図1.5
図1.5のウィンドウを閉じると、
の事後分布が表示される（図1.6）。
図1.6
図1.6のウィンドウを閉じると、信頼性係数
（式（2.2.5）の値
）の事後分布と平均値、四分位数が表示される（図1.7）。
図1.7
図1.7のウィンドウを閉じると、測定対象の特性の真値と項目カテゴリ値の和（尺度得点）との相関係数および決定係数の事後分布が表示される（図1.8）。これは、測定値と真値との関係を表すものであるが、測定値はカテゴリ化による情報の損失があるので、真値に対する測定値の説明率を表す決定係数は、潜在変数の真値に対する説明率を表すより小さい値となる（Okamoto, 2017）。
図1.8
図1.8のウィンドウを閉じるとプログラムの終了となる。終了後、端末には以下の統計量が出力されている。
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
K＝３のときのモデルに対するStanスクリプト（ファイル名Reli_3Cat.stan）を以下のように用意した。
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]);
}
}
}
基本的には、上のK>3の場合と同じであるが、カテゴリ境界は
、 
.
の２個で、フリーのカテゴリ境界のパラメータはない。上のスクリプトを用いるPythonスクリプトはこのウェブサイトの後半にリスト２として示した。
リスト２のプログラムを実行すると、入力データファイル名が求められる（図2.1）。
図2.1
入力データファイルはExcelファイルとして用意する（図2.2）。
・
・
・
図2.2
1行目に項目変数名を描く。第１列はデータの識別用変数であるが、これはデータの識別用で、通し番号でも文字列でもなんでもよい。第２列目から項目反応カテゴリ値のデータをかく。カテゴリ値は、カテゴリ数がKの場合、１からKまでの整数値で表す。K＝３であるので、カテゴリ値は１、２、３のいずれかである。逆転項目があれば、値を逆にして他の項目とカテゴリの方向を揃える。カテゴリ値は、測定概念の値と同順で、測定概念の特性が強いほどカテゴリ値も大きいとする。逆転項目のカテゴリ値を逆にするプログラムは、別のウェブサイトに用意した。
入力データファイル名が設定されてデータが読み込まれると、Stanによるサンプリングが始まる。サンプリングが終了すると、まずの事後分布が表示される（図2.3）。
図2.3
図2.3のウィンドウの右上のX印アイコンをクリックして閉じると、次にの事後分布が描かれる（図2.4）。
図2.4
図2.4のウィンドウを閉じると、の事後分布が表示される（図2.5）。
図2.5
図2.5のウィンドウを閉じると、信頼性係数（式（2.2.5）の値）の事後分布と平均値、四分位数が表示される（図2.6）。
図2.6
図2.6のウィンドウを閉じると、測定対象の特性の真値と項目カテゴリ値の和（尺度得点）との相関係数および決定係数の事後分布が表示される（図2.7）。これは、測定値と真値との関係を表すものであるが、測定値はカテゴリ化による情報の損失があるので、真値に対する測定値の説明率を表す決定係数は、潜在変数の真値に対する説明率を表すより小さい値となる（Okamoto, 2017）。
図2.7
図2.7のウィンドウを閉じるとプログラムの終了となる。終了後、端末には以下の統計量が出力されている。
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
K＝２のときのモデルに対するStanスクリプト（ファイル名Reli_2Cat.stan）を以下のように用意した。
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]);
}
}
}
基本的には、上のK>3の場合と同じであるが、カテゴリ境界が以下のように固定されている（Okamoto, 2013）。
また、とにも以下の制約が置かれている。
上のスクリプトを用いるPythonスクリプトは、このウェブサイトの後半にリスト３として示した。
リスト３のプログラムを実行すると、入力データファイル名が求められる（図3.1）。
図3.1
入力データファイルはExcelファイルとして用意する（図3.2）。
・
・
・
図3.2
1行目に項目変数名を描く。第１列はデータの識別用変数であるが、これはデータの識別用で、通し番号でも文字列でもなんでもよい。第２列目から項目反応カテゴリ値のデータをかく。カテゴリ値は、カテゴリ数がKの場合、１からKまでの整数値で表す。K＝２であるので、カテゴリ値は１か２のいずれかである。逆転項目があれば、値を逆にして他の項目とカテゴリの方向を揃える。カテゴリ値は、測定概念の値と同順で、測定概念の特性が強いほどカテゴリ値も大きいとする。逆転項目のカテゴリ値を逆にするプログラムは、別のウェブサイトに用意した。
入力データファイル名が設定されてデータが読み込まれると、Stanによるサンプリングが始まる。サンプリングが終了すると、まずの事後分布が表示される（図3.3）。
図3.3
図3.3のウィンドウの右上のX印アイコンをクリックして閉じると、信頼性係数（式（2.2.5）の値）の事後分布と平均値、四分位数が表示される（図3.4）。
図3.4
図3.4のウィンドウを閉じると、測定対象の特性の真値と項目カテゴリ値の和（尺度得点）との相関係数および決定係数の事後分布が表示される（図3.5）。これは、測定値と真値との関係を表すものであるが、測定値はカテゴリ化による情報の損失があるので、真値に対する測定値の説明率を表す決定係数は、潜在変数の真値に対する説明率を表すより小さい値となる（Okamoto, 2017）。
図3.5
図3.5のウィンドウを閉じるとプログラムの終了となる。終了後、端末には以下の統計量が出力されている。
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
リスト１ カテゴリ数K＞３の場合のスクリプト。スクリプトファイルなどは圧縮ファイルfiles8_G3Cat.zipにまとめた。
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))
リスト２ カテゴリ数K＝3の場合のスクリプト。スクリプトファイルなどは圧縮ファイルfiles8_3Cat.zipにまとめた。
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))
リスト３ カテゴリ数K＝２の場合のスクリプト。スクリプトファイルなどは圧縮ファイルfiles8_2Cat.zipにまとめた。
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.