Up

心理尺度の信頼性係数ω

 

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

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

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

 

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

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

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

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

の平均が0、分散が1であることから式(3)は

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

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

と書ける。

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

 

 

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

 

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

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

 

        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を実行すると、入力データファイル名の入力が求められる(図1)。

図1

 

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

2.a 入力データファイル例(Data10.xlsx

図2.b 入力データファイル例(Data10.xlsx

 

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

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

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

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

図3

 

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

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

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

図4

 

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

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

図5

 

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

図6

 

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

 

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で説明している。

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

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

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

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

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

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

信頼性係数アルファが

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

である。

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

 

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

 

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

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

A.1 項目Item_6Item_9を除外(ファイル名Data10_69.xlsx

 

項目Item_6と項目Item_92行目のセルに0を設定することにより、データ識別変数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])

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

 

 

 

直接最小2乗法を用いる場合(特異値分解の利用)

 

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

上のモデル(1)および(2)

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

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

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

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

また、

と置けば、

より、

を得る。したがって、

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

 

 

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

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

 

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

 

Input Excel data file (*.xlsx) =

 

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

SVD

 

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

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

 

次の例では、図SVD1Excelファイル名が設定されている。入力データファイル名を入力すると、次に出力ファイル名の設定が求められる。出力ファイル名は任意であるが、テキストファイルなので、ファイル名拡張子を「.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番目の特異値が特に大きい値であり、2番目以降はなだらかに減少している。モデル(1)、(2)の1因子モデルが成り立っていると考えられる。

SVD2を閉じると、次に図SVD3の散布図が表示される。

SVD

 

モデルによる予測値

と、データ値

の散布図である。

 

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

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

 

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用のStanPyStan)については、岡本(2019で説明している。

リストB1Stanスクリプトは、式(1)

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

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

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

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

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

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

 

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

Input Excel data file (*.xlsx) =

 

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

B1

 

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

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

 

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)における他とは著しく異なる項目は、上の統計量から6番目と9番目の項目であることが解る。

 

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

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

 

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

 

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

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

 

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

C2

 

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

C3

 

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

C4

 

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

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

 

 

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

10項目のときの

 

         Q1 = 0.897,  Med = 0.904,  Q3 = 0.911

 

から、第6項目と第9項目を抜いた8項目でのωは

 

         Q1 = 0.931,  Med = 0.936,  Q3 = 0.941

 

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

 

 

 

リスト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  リストB1Stanスクリプトを用いる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