Up

心理尺度の信頼性係数ω

尺度の分析を、主因子法のアルゴリズムを適用する場合、直接最小2乗法を適用して特異値分解を利用する場合、ベイズ分析を用いる場合を扱う。

主因子法は途中の繰り返し計算が収束しないことがあるが、特異値分解は常に解が得られるようである。

ベイズ分析では正規分布が仮定されているが、特異値分解による方法の場合は平均が０、分散が１という位置とサイズについての仮定だけで、分布の形状についての仮定はない。

心理尺度の信頼性係数としてオメガω（McDonald, 1999）は項目jの得点の因子モデル

に基づいて算出される（岡本、2014、第８章）。ここで、は測定対象の値（真値）、は項目jに固有の特性など誤差項としてまとめられるものである。は平均０、分散１に標準化されたものを考えるが、は素点である。

項目が項目あるとし、それらの合計点がテスト得点として算出されるとする。

である。このとき、オメガは次式で与えられる。

の平均が０、分散が１であることから式（３）は

と書ける。式番号（6.20a）はMcDonald (1999)におけるものである。

さらに、各が他の変数と独立であるとすると、式（6.20a）は

と書ける。

各項目の真値の係数（因子負荷量）が同じであるとき、アルファ係数の公式が導かれる。

主因子法のアルゴリズムを用いる場合

モデル（１）におけるの推定を、主因子法のアルゴリズムで行う。すなわち、因子分析における主因子法の計算方法を共分散行列に適用して因子負荷量を推定しする。詳しくは、＜岡本安晴「心理学データ分析と測定」勁草書房＞の「第８章　心理学的特性を測る：心理尺度」を参照されたい。ベイズ分析のスクリプト例は、このウェブサイトの後半に示した。

主因子法のアルゴリズムを用いるコードは、以下のように簡単である。

psi2 = numpy.zeros(n_vars)

A = numpy.zeros(n_vars)

while True:

A0 = A.copy()

C = S.copy()

for i in range(n_vars):

C[i][i] -= psi2[i]

w, vctr = numpy.linalg.eigh(C)

A = vctr[:, -1] * (w[-1] ** 0.5)

for i in range(n_vars):

psi2[i] = S[i][i] - (A[i] ** 2)

diff = ((A0 - A) ** 2).sum()

if diff < 1.0e-7:

break

プログラムはサンプルデータファイルとともに圧縮ファイルfilesEstOmega.zipにまとめた。ダウンロード・展開すれば自由に利用できる。

なお、Jupyter Notebook用のプログラムも用意した。圧縮ファイルfilesEstOmega_jn.zipにまとめている。Jupyter Notebookでの実行画面例をここに示した。

プログラムはPythonで書かれているが、Pythonの簡単なインストール法と使用法を説明したウェブサイトも用意している。

圧縮ファイルをダウンロード展開して作成されるフォルダにあるプログラムEstOmegaMain.pyを実行すると、入力データファイル名の入力が求められる（図１）。

図１

入力データファイルはExcelファイルとして用意する。図２は、10項目Item_1からItem_10に対する50人分のデータを表している。

図2.a　入力データファイル例（Data10.xlsx）

図２.b　入力データファイル例（Data10.xlsx）

１行目には変数名が書かれている。第１列目はデータの識別値であり、数値でも文字列でも何でもよい。識別値であるが、同じ値が複数回用いられてもよい。

２行目には、その変数が信頼性係数の計算において含めるかどうかを示すものである。整数値１が設定されている変数が尺度を構成する変数として含められ、整数値０が設定されていると尺度から除かれる。一部の変数を対象とする分析例は、ウェブサイトの後半で説明している。

データ値は３行目以降に書かれる。データ値は、逆転項目が用いられている場合は、逆転項目の値を逆転して他の項目と値の方向を合わせておく必要がある。

入力データファイル名を設定してEnterキーを押すと、次に出力用テキストデータファイル名の設定が求められる（図３）。

図３

図３では、入力データファイル名として図２のファイルData10.xlsxが設定され、出力用テキストファイル名としてResults.txtが設定されている。出力用テキストファイル名は任意であるが、ファイル拡張子が「.txt」であることに注意。

出力用テキストファイル名を入力すると計算が始まる。

計算が終了すると、出力用テキストファイルが保存されたことを示すメッセージが表示されて、プログラムの実行終了となる。

図４

実行終了後、出力用テキストファイルを開くことができる。

出力ファイルをPython IDLEで開こうとすると、図５のダイアログボックスが表示されることがある。

図５

この場合は、「s-jis」を設定して（図６）、「OK」ボタンをクリックすればよい。

図６

出力用テキストファイルの内容は、以下のようになっている。説明は、このリストの後で行っている。

Input data file = Data10.xlsx

n_data =50

n_vars = 10

var_names: ['Item_1', 'Item_2', 'Item_3', 'Item_4', 'Item_5', 'Item_6', 'Item_7', 'Item_8', 'Item_9', 'Item_10']

α = 0.879

Lambda[1]**2 / Var(X[1]) = 0.620

Lambda[2]**2 / Var(X[2]) = 0.682

Lambda[3]**2 / Var(X[3]) = 0.553

Lambda[4]**2 / Var(X[4]) = 0.702

Lambda[5]**2 / Var(X[5]) = 0.664

Lambda[6]**2 / Var(X[6]) = 0.010

Lambda[7]**2 / Var(X[7]) = 0.652

Lambda[8]**2 / Var(X[8]) = 0.746

Lambda[9]**2 / Var(X[9]) = 0.011

Lambda[10]**2 / Var(X[10]) = 0.714

Item_1: Lambda[1] = 12.990

Item_2: Lambda[2] = 12.642

Item_3: Lambda[3] = 12.765

Item_4: Lambda[4] = 16.076

Item_5: Lambda[5] = 13.397

Item_6: Lambda[6] = -1.719

Item_7: Lambda[7] = 12.301

Item_8: Lambda[8] = 15.702

Item_9: Lambda[9] = 1.611

Item_10: Lambda[10] = 15.803

GFI = 0.995

ω = 0.907

出力ファイルの説明

入力データファイル名、データ数、変数名に続いて、信頼性係数アルファαが出力されている。このαの値は、McDonald (1999)がGuttman-Cronbach alphaとも呼んでいるもので、いわゆるアルファ係数である。アルファ係数については、＜岡本安晴「心理学データ分析と測定」勁草書房＞のpp. 126-128で説明している。

続いて、変数ごとに各変数における真値の分散の各変数の分散に対する比

の値が出力されている。これは、各変数に占める真値の割合と考えられる。

上の出力では、６番目と９番目の値が低い。これらは尺度から除くのが望ましい変数と考えられる。

次に、各の値が出力されている。モデル（１）が素点に対する因子モデルであるので、因子負荷量の値は１より大きくなり得る。６番目の項目のの値が負であるが、これは変数の値の動きが真値と逆であることを意味しているので、逆転項目の可能性がある。また、値の絶対値も小さいので尺度から除くことも考えられる。

次に、GFIが示されているが、これはモデル（１）の適合度（goodness of fit index）を表し（McDonald (1999)、GFIは１に近いほど適合度が高く、0に近いほど適合度が低い（岡本安晴「心理学データ分析と測定」勁草書房の「第８章　心理学的特性を測る」、pp. 124-125参照）。

最後にの値が出力されている。式（6.20a）による値である。

信頼性係数アルファが

であるのに対し、信頼性係数ωは

である。

信頼性係数αが信頼性係数ωより小さい値であることが確認できる。

入力データファイルの一部の変数を選ぶ場合

上の出力におけるの値を見ると、項目Item-6の値は負であり、Item-9の値は他と比べて小さい。これらの項目を除いて以下のように再計算を行った。

Item-6とItem-9をアルファ係数の算出から除くため、上の入力用ファイルを図A.1に示すように変更する。

図A.1　項目Item_6とItem_9を除外（ファイル名Data10_69.xlsx）

項目Item_6と項目Item_9の2行目のセルに０を設定することにより、データ識別変数DataIDの値とともに分析対象の変数から除外する。

プログラム実行時の入力データファイルとして、図A.1のファイルData10_69.xlsxを指定する（図A.2）。

図A.2

このときの出力ファイルは、以下のようになる。

Input data file = Data10_69.xlsx

n_data =50

n_vars = 8

var_names: ['Item_1', 'Item_2', 'Item_3', 'Item_4', 'Item_5', 'Item_7', 'Item_8', 'Item_10']

α = 0.940

Lambda[1]**2 / Var(X[1]) = 0.618

Lambda[2]**2 / Var(X[2]) = 0.683

Lambda[3]**2 / Var(X[3]) = 0.546

Lambda[4]**2 / Var(X[4]) = 0.709

Lambda[5]**2 / Var(X[5]) = 0.659

Lambda[6]**2 / Var(X[6]) = 0.653

Lambda[7]**2 / Var(X[7]) = 0.749

Lambda[8]**2 / Var(X[8]) = 0.715

Item_1: Lambda[1] = 12.968

Item_2: Lambda[2] = 12.648

Item_3: Lambda[3] = 12.685

Item_4: Lambda[4] = 16.145

Item_5: Lambda[5] = 13.350

Item_7: Lambda[6] = 12.306

Item_8: Lambda[7] = 15.736

Item_10: Lambda[8] = 15.817

GFI = 0.999

ω = 0.942

信頼性係数ωは

となり、すべての変数を対象とした最初の信頼積係数ω

より、大きくなっている。

また、各変数における真値の分散の変数の分散に対する比

Lambda[j]**2 / Var(X[j])

に特に小さい値も見当たらない。

直接最小２乗法を用いる場合（特異値分解の利用）

主因子法のアルゴリズムの場合、繰り返し計算が収束しない場合がある。これに対して、特異値分解は常に成功することが期待される。

上のモデル（１）および（２）

において、特異値分解に定数項の影響が負荷量の推定に影響を与えるのを避けるため、

いま、個人の得点をで、特性値をで表し、モデル（PC.1）を次式で表す。Pythonスクリプトでは、Tは紛らわしいので避けた次第である。

モデル（PC.2）において、誤差項を最小2乗の意味で最小にするとを求める。

いま、の特異値分解を次式の形で表す。

また、

と置けば、

より、

を得る。したがって、

と置けば、は平均０，分散であり、とは式（PC.2）の最小2乗解である。因子の平均が０，分散が1ということは、尺度の原点と単位を決めるもので、本質的な制約ではない。

上のモデルに基づくプログラムをリストSVD1のように用意した。プログラムファイルなどは圧縮ファイルfilesEstOmegaSVD.zipとしてまとめた。

なお、プログラムはJupyter Notebook用も作成して、圧縮ファイルfilesEstOmegaSVD_jn.zipとしてまとめた。Jupyter Notebookでの実行例の画面はここに用意した。

リストSVD1を実行すると、入力データファイル名の入力が求められる。

Input Excel data file (*.xlsx) =

入力データファイルは、Excelファイル（*.xlsx）として図SVD１のように用意する。

図SVD１

１行目に変数名を置く。２行目は変数の選択を示す。１であれば分析に含められるが、１以外の変数は分析から除外される。図SVD1では、変数DataID、Item_6、Item_9が分析から除かれている。

３行目以降にデータ値を書く。

次の例では、図SVD1のExcelファイル名が設定されている。入力データファイル名を入力すると、次に出力ファイル名の設定が求められる。出力ファイル名は任意であるが、テキストファイルなので、ファイル名拡張子を「.txt」とする。

Input Excel data file (*.xlsx) = Data10_69.xlsx

Output text data file (*.txt) =

次の例では、Results.txtと入力している。

Input Excel data file (*.xlsx) = Data10_69.xlsx

Output text data file (*.txt) = Results.txt

出力ファイル名を入力すると、計算が始まり、特異値分解における特異値のスクリープロットが表示される（図SVD2）。

図SVD2

１番目の特異値が特に大きい値であり、２番目以降はなだらかに減少している。モデル（１）、（２）の1因子モデルが成り立っていると考えられる。

図SVD２を閉じると、次に図SVD３の散布図が表示される。

図SVD３

モデルによる予測値

と、データ値

の散布図である。

図SVD３を閉じると、プログラムの終了である。

端末には、以下のようなメッセージが表示されている。

Results_Scree.png was saved.

ResultsScatter.png was saved.

Results.txt was saved.

スクリープロットがResults_Scree.pngに、散布図がResultsScatter.pn保存されていることを意味している。これらの名前は、出力ファイル名を基に作成されている。

出力ファイルは、以下のような内容である。

Input data file = Data10_69.xlsx

n_data =50

n_vars = 8

var_names: ['Item_1', 'Item_2', 'Item_3', 'Item_4', 'Item_5', 'Item_7', 'Item_8', 'Item_10']

Item_1 ==> X[1]

Item_2 ==> X[2]

Item_3 ==> X[3]

Item_4 ==> X[4]

Item_5 ==> X[5]

Item_7 ==> X[6]

Item_8 ==> X[7]

Item_10 ==> X[8]

Lambda[1]**2 / Var(X[1]) = 0.648

Lambda[2]**2 / Var(X[2]) = 0.690

Lambda[3]**2 / Var(X[3]) = 0.595

Lambda[4]**2 / Var(X[4]) = 0.747

Lambda[5]**2 / Var(X[5]) = 0.682

Lambda[6]**2 / Var(X[6]) = 0.665

Lambda[7]**2 / Var(X[7]) = 0.766

Lambda[8]**2 / Var(X[8]) = 0.746

Item_1: Lambda[1] = 13.280

Item_2: Lambda[2] = 12.717

Item_3: Lambda[3] = 13.244

Item_4: Lambda[4] = 16.572

Item_5: Lambda[5] = 13.578

Item_7: Lambda[6] = 12.422

Item_8: Lambda[7] = 15.906

Item_10: Lambda[8] = 16.154

GFI = 0.997

ω = 0.999

分析された変数名が、プログラムでの変数として何番目であるかの一覧が示されている。

最後に、ωの値も書き出されている。

リストSVD1　十九一分解によるωの計算

import numpy

import pandas as pd

import matplotlib.pyplot as plt

s = input("Input Excel data file (*.xlsx) = ")

xlsx_fl = pd.ExcelFile(s)

data_xlsx = pd.read_excel(xlsx_fl, header = None)

data = data_xlsx.values

s_out = input("Output text data file (*.txt) = ")

fout = open(s_out, "w")

fout.write("Input data file = " + s + "\n")

data = data.T

n_vars = 0

var_names = []

X = []

for v in data:

if v[1] == 1:

var_names.append(v[0])

X.append(list(v[2:]))

n_vars += 1

n_data = len(X[0])

print('n_data =', n_data)

print('n_vars =', n_vars)

print('var_names = :', var_names)

for i, nm in enumerate( var_names):

print(f'{nm:<.10s} ==> X[{i+1}]')

print()

fout.write('\nn_data ={}\n'.format(n_data))

fout.write('n_vars = {}\n'.format(n_vars))

fout.write('var_names: {}\n'.format(var_names))

for i, nm in enumerate( var_names):

fout.write(f'{nm:<.10s} ==> X[{i+1}]\n')

fout.write('\n')

X = numpy.array(X).T

RawScore = numpy.sum(X, axis=1)

meanX = numpy.mean(X, axis=0)

X = X - meanX # Centering of X

#

# SVD of the centered X

#

U, Sv, Vh = numpy.linalg.svd(X, full_matrices=False)

#

# Scree plot of sigular values

#

plt.plot(numpy.arange(1, len(Sv)+1), Sv)

plt.title('Scree Plot', fontsize= 16)

plt.ylabel('Singular value', fontsize=14)

plt.tight_layout()

fnm_scree = s_out[:-4] + '_Scree.png'

plt.savefig(fnm_scree)

plt.show()

if numpy.sum(Vh[0]) < 0.0:

U *= -1.0

Vh *= -1.0

S = numpy.cov(X, rowvar=False)

A = Vh[0] * Sv[0] / (n_data**0.5) # Vector of Coefficients, i.e. Lambda's

for i in range(n_vars):

print('Lambda[{0}]**2 / Var(X[{0}]) = {1:.3f}'.format(

i+1, A[i]**2 / S[i][i]))

fout.write('Lambda[{0}]**2 / Var(X[{0}]) = {1:.3f}\n'.format(

i+1, A[i]**2 / S[i][i]))

print()

fout.write('\n')

for i in range(n_vars):

print('{0}: Lambda[{1}] = {2:.3f}'.format(

var_names[i], i+1, A[i]))

fout.write('{0}: Lambda[{1}] = {2:.3f}\n'.format(

var_names[i], i+1, A[i]))

#

# Calculation of GFI (Goodness of Fit Index)

#

sum_num = 0.0

sum_den = 0.0

for i in range(n_vars):

for j in range(n_vars):

if i == j:

sum_den += S[i][i] ** 2

else:

sum_den += S[i][j] ** 2

sum_num += (S[i][j] - A[i] * A[j]) ** 2

GFI = 1.0 - sum_num / sum_den

print('\nGFI = {0:.3f}'.format(GFI))

fout.write('\nGFI = {0:.3f}\n'.format(GFI))

#

# Calculation of omega

#

sx2 = numpy.var(X.sum(axis=1))

omega = (A.sum() ** 2) / sx2

print('\nω = {0:.3f}'.format(omega))

fout.write('\nω = {0:.3f}\n'.format(omega))

f = U[:,0] * (n_data**0.5)

f = numpy.reshape(f, (n_data,1))

A = numpy.reshape(A, (1, len(A)))

ModelScore = numpy.sum(f @ A, axis=1) + numpy.sum(meanX)

#

# Scatter plot of (model prediction, raw score)

#

plt.scatter(ModelScore, RawScore)

plt.xlabel('Model Score', fontsize=14)

plt.ylabel('Raw Score', fontsize=14)

plt.plot([numpy.min(ModelScore), numpy.max(ModelScore)],[numpy.min(ModelScore), numpy.max(ModelScore)],

label = '$Y=X$', color='k')

plt.title(f'$\\omega={omega:.3f}$', fontsize=16)

plt.legend(fontsize=14)

plt.tight_layout()

fnm_scatter = s_out[:-4] + 'Scatter.png'

plt.savefig(fnm_scatter)

plt.show()

fout.close()

print()

print(fnm_scree + ' was saved.')

print(fnm_scatter + ' was saved.')

print(s_out + " was saved.\n")

ベイズ分析例

ベイズ分析のためのStanスクリプトをリストB1に、リストB1のスクリプトを用いてベイズ分析を行うPythonスクリプトをリストB2に示した。これらのファイルはまとめて圧縮ファイルfilesPS3.zipとした。ダウンロード・解凍すれば、自由に利用できる。ベイズ分析とPython用のStan（PyStan）については、岡本（2019）で説明している。

リストB1のStanスクリプトは、式（１）

のモデルに対するものであるが、式（１）におけるTはリストB1ではFで表されている。

パラメータFは、平均０，標準偏差１の標準正規分布に従うという仮定が置かれているが、事前分布での制約ではshrinkageが起こり、このためが過大評価されるという傾向が見られる。これを避けるために、Fの標準化を行っている。

また、MCMCサンプリングにおいてlabel switchingを避けるためにパラメータに非負の制約を置いている。

事前分布は、soft constraintを用いたが、Gelmanら（2021、p. 124）より緩いものを用いた。これは、Gelmanらでは、データ数がsparseの場合にも対応できるように設定されているが、尺度の信頼性を求めるときのデータはある程度の多さがあると考えられるからである。データ数が少なくてパラメータ推定が不安定な場合は、それに応じて事前分布を調整すればよい。PyStan 3では、事前分布として分散の大きいものを設定すると警告が出るので、一様分布に近いものとして分散の大きい分布を用いるときは、まず警告の出ない事前分布のパラメータを用意して、これを変換して分散の大きい分布としている。

リストB1およびリストB2はPyStan 3対応である。PyStan 3をWindowsで利用するときは、WindowsにUbuntuをインストールしてUbuntu上で利用するが、UbuntuのWindowsへのインストール、およびUbuntuでのPythonの準備は簡単なコマンドで出来る。このウェブサイトで説明しているので、必要ならば参照されたい。

リストB2のスクリプトを実行すると、入力データファイル名を聞いてくる。

(py311) ・・・filesPS3$ python EstOmega.py

Input Excel data file (*.xlsx) =

入力データファイルはExcelファイルとして図B1のように用意する。

図B1

上の主因子法のアルゴリズムを用いたスクリプト用データファイル（図２）のものと同様であるが、図B1の場合は、データは２行目から設定されている。１行目に変数名を設定して、次の２行目からデータを設定する。ただし、１列目はデータの識別用値を設定する。１列目の値は、分析では無視される。

入力データファイル名を次のように設定する。

Input Excel data file (*.xlsx) = Data10PS3.xlsx

ファイル名を設定すると、計算が始まる。MCMCサンプリングが終了すると、まずトレース図が表示される（図B2）。

図B2

図B2のフォームを閉じると、次にの事後分布のグラフが表示される（図B3）。

図B3

図B3のフォームを閉じると、の事後分布が表示される（図B4）。

図B4

図B4のフォームを閉じると、の事後分布のグラフが表示される（図B5）。

図B5

図B5のフォームを閉じると、次に、変数ごとに各変数における真値の分散の各変数の分散に対する比

の事後分布のグラフが表示される。この値は、各変数に占める真値の割合と考えられる。

図B6

図B6のフォームを閉じると、ωの事後分布のグラフが表示される（図B7）。

図B7

グラフのタイトル部に、事後分布の4分位数が表示されている。点推定値としては、MAP推定値（モード）がよいという報告（岡本、2022）があるが、計算時間が掛かるので、ここではGelmanら（2021）に従って、中央値を点推定値とした。

図B7のフォームを閉じると、図B8の散布図が表示される。

図B8

横軸に各の推定値（事後分布の中央値）、縦軸にテストスコアの和をとり、プロットしたものである。点がほぼ直線沿いに散らばっていることが解る。

図B8のフォームを閉じると、プログラムの実行終了となる。端末画面には、次のメッセージが表示されている。

Results.txt was saved.

(py311) ・・・/filesPS3$

出力ファイルResults.txtを開くと、以下のようになっている。

Input data file = Data10PS3.xlsx

data...

[[ 1 79 46 68 65 89 38 58 70 49 65]

[ 2 48 71 61 50 43 27 48 54 34 61]

[ 3 38 42 38 45 55 43 47 39 65 54]

　　　・

　　　・

　　　・

[48 36 46 49 54 31 75 43 40 30 49]

[49 44 24 56 51 58 34 42 52 67 48]

[50 71 75 64 84 58 25 67 65 47 56]]

Summary...

mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk \

Lmbd[0] 12.987 1.574 10.091 15.940 0.020 0.014 6152.0

Lmbd[1] 12.569 1.370 9.967 15.139 0.019 0.013 5338.0

Lmbd[2] 12.684 1.804 9.237 15.933 0.023 0.017 5933.0

　　　・

　　　・

　　　・

F[47] -0.617 0.239 -1.052 -0.170 0.003 0.002 7217.0

F[48] -0.415 0.251 -0.891 0.062 0.003 0.003 7061.0

F[49] 0.987 0.248 0.539 1.453 0.003 0.002 7256.0

ess_tail r_hat

Lmbd[0] 2580.0 1.0

Lmbd[1] 2604.0 1.0

Lmbd[2] 2870.0 1.0

　　　・

　　　・

　　　・

F[47] 3045.0 1.0

F[48] 2993.0 1.0

F[49] 3036.0 1.0

Mu...

Mu[1]: Q1 = 50.755, Med = 51.759, Q3 = 52.775

Mu[2]: Q1 = 53.383, Med = 54.236, Q3 = 55.100

Mu[3]: Q1 = 49.078, Med = 50.193, Q3 = 51.336

Mu[4]: Q1 = 51.594, Med = 52.628, Q3 = 53.611

Mu[5]: Q1 = 51.900, Med = 52.866, Q3 = 53.825

Mu[6]: Q1 = 44.403, Med = 46.027, Q3 = 47.688

Mu[7]: Q1 = 51.587, Med = 52.501, Q3 = 53.413

Mu[8]: Q1 = 51.408, Med = 52.296, Q3 = 53.215

Mu[9]: Q1 = 42.277, Med = 43.762, Q3 = 45.190

Mu[10]: Q1 = 53.024, Med = 54.027, Q3 = 55.042

Sigma...

Sigma[1]: Q1 = 9.662 Med = 10.441, Q3 = 11.333

Sigma[2]: Q1 = 8.263 Med = 8.963, Q3 = 9.735

Sigma[3]: Q1 = 11.047 Med = 11.854, Q3 = 12.786

Sigma[4]: Q1 = 9.848 Med = 10.695, Q3 = 11.594

Sigma[5]: Q1 = 9.242 Med = 9.987, Q3 = 10.835

Sigma[6]: Q1 = 16.279 Med = 17.439, Q3 = 18.742

Sigma[7]: Q1 = 8.612 Med = 9.306, Q3 = 10.086

Sigma[8]: Q1 = 8.594 Med = 9.398, Q3 = 10.244

Sigma[9]: Q1 = 14.562 Med = 15.635, Q3 = 16.777

Sigma[10]: Q1 = 9.434 Med = 10.249, Q3 = 11.132

Lambda...

Lambda[1]: Q1 = 11.953 Med = 13.019, Q3 = 14.060 (Item_1)

Lambda[2]: Q1 = 11.668 Med = 12.595, Q3 = 13.499 (Item_2)

Lambda[3]: Q1 = 11.473 Med = 12.677, Q3 = 13.932 (Item_3)

Lambda[4]: Q1 = 15.082 Med = 16.157, Q3 = 17.192 (Item_4)

Lambda[5]: Q1 = 12.246 Med = 13.262, Q3 = 14.256 (Item_5)

Lambda[6]: Q1 = 0.573 Med = 1.284, Q3 = 2.244 (Item_6)

Lambda[7]: Q1 = 11.354 Med = 12.252, Q3 = 13.180 (Item_7)

Lambda[8]: Q1 = 14.757 Med = 15.747, Q3 = 16.738 (Item_8)

Lambda[9]: Q1 = 1.162 Med = 2.267, Q3 = 3.578 (Item_9)

Lambda[10]: Q1 = 14.782 Med = 15.840, Q3 = 16.830 (Item_10)

Lambda_j**/Var(X_j)...

Lambda[1]**2 / Var(X[1]): Q1 = 0.548, Med = 0.608, Q3 = 0.661

Lambda[2]**2 / Var(X[2]): Q1 = 0.608, Med = 0.664, Q3 = 0.712

Lambda[3]**2 / Var(X[3]): Q1 = 0.468, Med = 0.534, Q3 = 0.591

Lambda[4]**2 / Var(X[4]): Q1 = 0.646, Med = 0.696, Q3 = 0.740

Lambda[5]**2 / Var(X[5]): Q1 = 0.583, Med = 0.638, Q3 = 0.689

Lambda[6]**2 / Var(X[6]): Q1 = 0.001, Med = 0.005, Q3 = 0.016

Lambda[7]**2 / Var(X[7]): Q1 = 0.578, Med = 0.634, Q3 = 0.684

Lambda[8]**2 / Var(X[8]): Q1 = 0.690, Med = 0.737, Q3 = 0.779

Lambda[9]**2 / Var(X[9]): Q1 = 0.005, Med = 0.021, Q3 = 0.051

Lambda[10]**2 / Var(X[10]): Q1 = 0.655, Med = 0.705, Q3 = 0.748

omega...

Q1 = 0.897, Med = 0.904, Q3 = 0.911

まず、入力データファイル名が出力されている。

続いて、MCMCサンプリングの統計量が出力されている。

事後分布などのグラフで表示されていた情報が、統計量の数値データとして読むことができる。

図B5（）、図B6（）における他とは著しく異なる項目は、上の統計量から６番目と９番目の項目であることが解る。

図B1のデータから、６番目と９番目の項目のデータを抜いたものを、図C1に示すように用意した。

図C1　　図B1のデータから６番目と９番目のデータを抜いたもの

図C1のデータに対して分析を行った。

(py311) ・・・/filesPS3$ python EstOmega.py

Input Excel data file (*.xlsx) = Data10PS3_69.xlsx

トレース図が、図C2のように表示された。

図C2

の事後分布は、図C3のようである。

図C3

の事後分布は図C4のようである。

図C4

の事後分布は、図Ｃ５のようである。

図C5

の事後分布は図C6のようである。

図C6

ωの事後分布は図C7のようである。図B7のものより右に寄っていることが分かる。

図C7

点の散布図は図C8のようである。図B8のものより、点が狭い帯内にまとまっていることが分かる。

図C8

出力ファイルResults.txtを開くと、以下のようになっている。

Input data file = Data10PS3_69.xlsx

data...

[[ 1 79 46 68 65 89 58 70 65]

[ 2 48 71 61 50 43 48 54 61]

[ 3 38 42 38 45 55 47 39 54]

　　　・

　　　・

　　　・

[48 36 46 49 54 31 43 40 49]

[49 44 24 56 51 58 42 52 48]

[50 71 75 64 84 58 67 65 56]]

Summary...

mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk \

Lmbd[0] 13.013 1.544 10.169 15.920 0.025 0.018 3958.0

Lmbd[1] 12.601 1.357 9.989 15.217 0.019 0.014 5049.0

Lmbd[2] 12.653 1.744 9.191 15.731 0.025 0.018 4911.0

Lmbd[3] 16.122 1.599 13.280 19.305 0.023 0.017 4707.0

Lmbd[4] 13.229 1.502 10.362 15.968 0.023 0.016 4469.0

Lmbd[5] 12.273 1.370 9.569 14.739 0.021 0.015 4444.0

Lmbd[6] 15.738 1.416 13.138 18.418 0.020 0.014 4846.0

Lmbd[7] 15.794 1.548 12.671 18.533 0.022 0.016 4902.0

mu[0] 51.790 1.500 48.967 54.636 0.020 0.014 5585.0

mu[1] 54.193 1.299 51.846 56.753 0.018 0.012 5535.0

mu[2] 50.176 1.742 46.900 53.454 0.023 0.016 5982.0

mu[3] 52.680 1.562 49.715 55.525 0.020 0.014 5853.0

mu[4] 52.885 1.444 50.044 55.507 0.019 0.014 5559.0

mu[5] 52.545 1.339 49.991 54.953 0.017 0.012 5808.0

mu[6] 52.262 1.372 49.468 54.673 0.017 0.012 6748.0

mu[7] 54.047 1.447 51.325 56.794 0.020 0.014 5425.0

sgm[0] 10.492 1.239 8.374 12.929 0.020 0.014 3904.0

sgm[1] 9.006 1.103 7.091 11.146 0.018 0.013 3777.0

sgm[2] 11.991 1.401 9.464 14.560 0.022 0.016 4343.0

sgm[3] 10.813 1.378 8.393 13.415 0.022 0.016 4145.0

sgm[4] 10.114 1.200 7.917 12.302 0.019 0.013 4231.0

sgm[5] 9.343 1.093 7.442 11.483 0.017 0.012 4417.0

sgm[6] 9.454 1.235 7.271 11.876 0.021 0.015 3445.0

sgm[7] 10.335 1.295 7.965 12.680 0.021 0.015 3747.0

omg 0.936 0.007 0.923 0.949 0.000 0.000 3502.0

F[0] 0.938 0.255 0.464 1.425 0.003 0.002 5820.0

F[1] 0.130 0.243 -0.346 0.569 0.003 0.003 6219.0

F[2] -0.517 0.238 -0.982 -0.093 0.003 0.002 5968.0

　　　・

　　　・

　　　・

F[47] -0.626 0.249 -1.092 -0.160 0.003 0.002 6323.0

F[48] -0.434 0.253 -0.925 0.025 0.003 0.003 5618.0

F[49] 0.995 0.243 0.547 1.435 0.003 0.002 6610.0

ess_tail r_hat

Lmbd[0] 2450.0 1.0

Lmbd[1] 2752.0 1.0

Lmbd[2] 2711.0 1.0

　　　・

　　　・

　　　・

F[47] 2721.0 1.0

F[48] 2869.0 1.0

F[49] 3485.0 1.0

Mu...

Mu[1]: Q1 = 50.834, Med = 51.762, Q3 = 52.781

Mu[2]: Q1 = 53.341, Med = 54.216, Q3 = 55.039

Mu[3]: Q1 = 48.997, Med = 50.174, Q3 = 51.340

Mu[4]: Q1 = 51.662, Med = 52.696, Q3 = 53.687

Mu[5]: Q1 = 51.951, Med = 52.902, Q3 = 53.827

Mu[6]: Q1 = 51.654, Med = 52.554, Q3 = 53.445

Mu[7]: Q1 = 51.374, Med = 52.262, Q3 = 53.167

Mu[8]: Q1 = 53.084, Med = 54.027, Q3 = 54.999

Sigma...

Sigma[1]: Q1 = 9.627 Med = 10.376, Q3 = 11.237

Sigma[2]: Q1 = 8.239 Med = 8.938, Q3 = 9.702

Sigma[3]: Q1 = 11.007 Med = 11.906, Q3 = 12.859

Sigma[4]: Q1 = 9.868 Med = 10.683, Q3 = 11.654

Sigma[5]: Q1 = 9.288 Med = 9.999, Q3 = 10.857

Sigma[6]: Q1 = 8.588 Med = 9.247, Q3 = 9.998

Sigma[7]: Q1 = 8.606 Med = 9.367, Q3 = 10.182

Sigma[8]: Q1 = 9.446 Med = 10.237, Q3 = 11.122

Lambda...

Lambda[1]: Q1 = 11.996 Med = 13.046, Q3 = 14.031 (Item_1)

Lambda[2]: Q1 = 11.732 Med = 12.607, Q3 = 13.518 (Item_2)

Lambda[3]: Q1 = 11.510 Med = 12.664, Q3 = 13.789 (Item_3)

Lambda[4]: Q1 = 15.088 Med = 16.150, Q3 = 17.189 (Item_4)

Lambda[5]: Q1 = 12.236 Med = 13.259, Q3 = 14.252 (Item_5)

Lambda[6]: Q1 = 11.389 Med = 12.296, Q3 = 13.175 (Item_7)

Lambda[7]: Q1 = 14.793 Med = 15.781, Q3 = 16.693 (Item_8)

Lambda[8]: Q1 = 14.775 Med = 15.811, Q3 = 16.828 (Item_10)

Lambda_j**/Var(X_j)...

Lambda[1]**2 / Var(X[1]): Q1 = 0.552, Med = 0.611, Q3 = 0.661

Lambda[2]**2 / Var(X[2]): Q1 = 0.611, Med = 0.667, Q3 = 0.713

Lambda[3]**2 / Var(X[3]): Q1 = 0.467, Med = 0.532, Q3 = 0.591

Lambda[4]**2 / Var(X[4]): Q1 = 0.643, Med = 0.695, Q3 = 0.738

Lambda[5]**2 / Var(X[5]): Q1 = 0.579, Med = 0.637, Q3 = 0.686

Lambda[6]**2 / Var(X[6]): Q1 = 0.586, Med = 0.638, Q3 = 0.685

Lambda[7]**2 / Var(X[7]): Q1 = 0.692, Med = 0.739, Q3 = 0.780

Lambda[8]**2 / Var(X[8]): Q1 = 0.655, Med = 0.705, Q3 = 0.747

omega...

Q1 = 0.931, Med = 0.936, Q3 = 0.941

他と特に異なる値の項目は認められない。ωの値も高くなっており

１０項目のときの

Q1 = 0.897, Med = 0.904, Q3 = 0.911

から、第６項目と第９項目を抜いた８項目でのωは

Q1 = 0.931, Med = 0.936, Q3 = 0.941

であり、四分位数の比較でも８項目に減らしたときのωの方が高い値であると言えよう。

リストB1　　ベイズ分析のためのStanスクリプト（est_omega.stan）

//

// Yasuharu Okamoto, 2023.08

//

data {

int N;

int M;

array[N,M] real X;

array[M] real i_mu;

array[M] real i_std;

}

parameters {

array[N] real vF;

array[M] real v_mu;

array[M] real<lower = 0.0> v_Lmbd;

array[M] real<lower = 0.0> v_sgm;

}

transformed parameters {

array[N] real F;

real sumF;

real ssF;

real meanF;

real sdF;

array[M] real mu;

array[M] real Lmbd;

array[M] real sgm;

array[N,M] real V;

//

// Standardized F values

//

sumF = 0.0;

ssF = 0.0;

for (i in 1:N) {

sumF += vF[i];

ssF += pow(vF[i], 2);

}

meanF = sumF / N;

sdF = pow((ssF - N * pow(meanF, 2)) / (N - 1), 0.5);

for (i in 1:N) {

F[i] = (vF[i] - meanF) / sdF;

}

//

// Rescaled parameters

//

for (j in 1:M) {

mu[j] = i_mu[j] + 10.0 * i_std[j] * v_mu[j];

Lmbd[j] = 10.0 * i_std[M] * v_Lmbd[j];

sgm[j] = 0.001 + 10.0 * i_std[M] * v_sgm[j];

}

//

// Regression model

//

for (i in 1:N) {

for (j in 1:M) {

V[i][j] = mu[j] + Lmbd[j] * F[i];

}

}

}

model {

for (i in 1:N) {

vF[i] ~ normal(0.0, 1.0);

}

for (j in 1:M) {

v_mu[j] ~ normal(0.0, 10.0);

v_Lmbd[j] ~ exponential(0.1);

v_sgm[j] ~ exponential(0.1);

}

for (i in 1:N) {

for (j in 1:M) {

X[i][j] ~ normal(V[i][j], sgm[j]);

}

}

}

generated quantities {

real omg; // omega

array[M] real L2V;

real v1;

real v2;

v1 = 0.0;

v2 = 0.0;

for (j in 1:M) {

L2V[j] = pow(Lmbd[j], 2) / (pow(Lmbd[j], 2) + pow(sgm[j], 2));

v1 += Lmbd[j];

v2 += pow(sgm[j], 2);

}

omg = pow(v1, 2) / (pow(v1, 2) + v2);

}

リストB2　　リストB1のStanスクリプトを用いるPythonスクリプト（EstOmega.py）

"""

Yasuharu Okamoto, 2023.08

"""

import numpy as np

import pandas as pd

import stan

import matplotlib.pyplot as plt

import seaborn as sb

import arviz as az

pd.set_option('display.max_columns', 10000)

pd.set_option('display.max_rows', 10000)

fout = open('Results.txt', 'w')

#

# Excelデータファイルの読み込み

#

s = input("Input Excel data file (*.xlsx) = ")

xlsx_fl = pd.ExcelFile(s)

data_xlsx = pd.read_excel(xlsx_fl)

var_names = [v for v in data_xlsx.keys()][1:]

print('var_names =', var_names)

data = data_xlsx.values

print('data...\n', data)

fout.write(f'Input data file = {s}\n')

fout.write(f'\ndata...\n{data}\n')

X = data[:,1:]

mus = np.mean(X, axis=0)

stds = np.std(X, axis=0)

Data = {'N': len(X), 'M': len(X[0]), 'X': X, 'i_mu': mus, 'i_std': stds}

with open('est_omega.stan', 'r') as fstan:

sm = stan.build(fstan.read(), data = Data)

fit = sm.sample()

i_data = az.from_pystan(posterior = fit, posterior_model = sm)

az.plot_trace(i_data, var_names = ['Lmbd', 'sgm', 'omg'])

plt.tight_layout()

plt.show()

smmry = az.summary(i_data, var_names = ['Lmbd', 'mu', 'sgm', 'omg', 'F'])

print('Summary...\n', smmry)

fout.write(f'\nSummary...\n{smmry}\n')

df = fit.to_frame()

print()

fout.write('\n\nMu...\n')

for j in range(len(X[0])):

q1, med, q3 = np.percentile(df[f'mu.{j+1}'], [25, 50, 75])

print(f'Mu[{j+1}]:', end = '')

print(f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}')

fout.write(f'Mu[{j+1}]: Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}\n')

for j in range(len(X[0])):

sb.kdeplot(df[f'mu.{j+1}'])

plt.title('$\mu$', fontsize = 16)

plt.xlabel('')

plt.show()

print()

fout.write('\nSigma...\n')

for j in range(len(X[0])):

q1, med, q3 = np.percentile(df[f'sgm.{j+1}'], [25, 50, 75])

print(f'Sigma[{j+1}]:', end = '')

print(f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}')

fout.write(f'Sigma[{j+1}]: Q1 = {q1:.3f}' +

f' Med = {med:.3f}, Q3 = {q3:.3f}\n')

for j in range(len(X[0])):

sb.kdeplot(df[f'sgm.{j+1}'])

plt.title('$\sigma$', fontsize = 16)

plt.xlabel('')

plt.show()

print()

fout.write('\nLambda...\n')

for j in range(len(X[0])):

q1, med, q3 = np.percentile(df[f'Lmbd.{j+1}'], [25, 50, 75])

print(f'Lambda[{j+1}]:', end = '')

print(f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}', end = '')

print(f' ({var_names[j]})')

fout.write(f'Lambda[{j+1}]: Q1 = {q1:.3f}' +

f' Med = {med:.3f}, Q3 = {q3:.3f}' +

f' ({var_names[j]})\n')

for j in range(len(X[0])):

sb.kdeplot(df[f'Lmbd.{j+1}'])

plt.title('$\lambda$', fontsize = 16)

plt.xlabel('')

plt.savefig('FigLmbd.png')

plt.show()

print()

fout.write('\nLambda_j**/Var(X_j)...\n')

for j in range(len(X[0])):

q1, med, q3 = np.percentile(df[f'L2V.{j+1}'], [25, 50, 75])

print(f'Lambda[{j+1}]**2 / Var(X[{j+1}]):', end = '')

print(f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}')

fout.write(f'Lambda[{j+1}]**2 / Var(X[{j+1}]):' +

f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}\n')

for j in range(len(X[0])):

sb.kdeplot(df[f'L2V.{j+1}'])

plt.title('$\lambda^{2}_{j}/Var(X_{j})$', fontsize = 16)

plt.xlabel('')

plt.show()

print()

q1, med, q3 = np.percentile(df[f'omg'], [25, 50, 75])

print(f'omega: Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}')

fout.write('\nomega...\n')

fout.write(f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}\n')

sb.kdeplot(df['omg'])

plt.title(f'Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}', fontsize = 14)

plt.xlabel('$\omega$', fontsize = 16)

plt.tight_layout()

plt.show()

sumXs = np.sum(X, axis = 1)

Fs = []

for i in range(len(X)):

Fs.append(np.median(df[f'F.{i+1}']))

plt.plot(Fs, sumXs, 'o')

plt.xlabel('$\hat{F}_{median}$', fontsize = 14)

plt.ylabel('$\sum{X_{j}}$', fontsize = 14)

plt.tight_layout()

plt.show()

fout.close()

print('\nResults.txt was saved.\n')

引用文献

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and Other Stories. Cambridge University Press.

McDonald, R. P. (1999). Test Theory: A unified treatment. Psychology Press.

Up