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")

 

 

 

 

ベイズ分析例

 

ベイズ分析例は、別のウェブサイトとして用意した。

 

 

 

 

引用文献

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

 

 

Up