Zenn
📙

一対比較法、サーストン法における内的整合性の検定

2021/01/17に公開
Python
統計学
心理統計
サーストン法
tech

一対比較法、サーストン法における内的整合性の検定

Python3を利用して一対比較法シリーズ最終回です。

前回は 一対比較法、サーストン法を用いて尺度化,プロットまでを行った。
https://zenn.dev/_kazuya/articles/a1179ed3f6e027

今回は算出したZ、尺度値を利用して、サーストン法における内的整合性の検定を行う。

Githubにも上げています↓
https://github.com/Iovesophy/paired_comparison_method.git

内的整合性の検定

前回、刺激に対して、尺度値を与えたならば、一体どのような判断の比率が期待されるか？
その比率は実験によって得られた値と十分に一致するか？

これらは私たちの仮説が正しければ十分一致するはずである。
このことから、尺度化の場合にとった逆の手順を辿ることによって、期待比率を求めることができる。

内的整合性の検定では、このようにして得られた各期待比率を、それに対応する実測比率と比較し、検定を用いて調べることである。[1]

手法について

まず、期待される尺度距離に基づく行列\hat Zを作る。
この期待されるという意味は、刺激について得られた尺度値に基づき、期待されるという意味である。

ここで、前回生成した比率行列から考えていく。
このサンプルデータはn=100で筆者が生成したものである。

表２.1

i>j A_j B_j C_j D_j E_j F_j G_j H_j
A_i 0.5 0.07 0.06 0.09 0.03 0.05 0.06 0.08
B_i 0.93 0.5 0.08 0.05 0.08 0.08 0.06 0.13
C_i 0.94 0.92 0.5 0.03 0.09 0.07 0.11 0.17
D_i 0.91 0.95 0.97 0.5 0.15 0.08 0.24 0.28
E_i 0.97 0.92 0.91 0.85 0.5 0.17 0.31 0.33
F_i 0.95 0.92 0.93 0.92 0.83 0.5 0.27 0.34
G_i 0.94 0.94 0.89 0.76 0.69 0.73 0.5 0.42
H_i 0.92 0.87 0.83 0.72 0.67 0.66 0.58 0.5

表２からZを算出

表3

i>j A_j B_j C_j D_j E_j F_j G_j H_j
A_i 0. -1.47579103 -1.55477359 -1.34075503 -1.88079361 -1.64485363 -1.55477359 -1.40507156
B_i 1.47579103 0. -1.40507156 -1.64485363 -1.40507156 -1.40507156 -1.55477359 -1.12639113
C_i 1.55477359 1.40507156 0. -1.88079361 -1.34075503 -1.47579103 -1.22652812 -0.95416525
D_i 1.34075503 1.64485363 1.88079361 0. -1.03643339 -1.40507156 -0.70630256 -0.58284151
E_i 1.88079361 1.40507156 1.34075503 1.03643339 0. -0.95416525 -0.49585035 -0.43991317
F_i 1.64485363 1.40507156 1.47579103 1.40507156 0.95416525 0. -0.61281299 -0.41246313
G_i 1.55477359 1.55477359 1.22652812 0.70630256 0.49585035 0.61281299 0. -0.20189348
H_i 1.40507156 1.12639113 0.95416525 0.58284151 0.43991317 0.41246313 0.20189348 0.

表3からZの合計を算出(ソート)

表4

i>j SUM
A_i -10.856812
B_i -7.065442
C_i -3.91818789
D_i 1.13575325
E_i 3.77312482
H_i 5.12273923
F_i 5.85967691
G_i 5.94914772

表4からZの平均値を算出（尺度値、ソート)

表5

i>j MEAN
A_i -1.35710151
B_i -0.88318025
C_i -0.48977349
D_i 0.14196916
E_i 0.4716406
H_i 0.6403424
F_i 0.73245961
G_i 0.74364347

合計\sum z_{ij}から、平均値M_{zij}を算出したら、この平均値を、刺激S_iの尺度値とみなすが、特に刺激S_0(刺激S_{ij}はここではA_{ij} ... H_{ij}のことを表す)を原点とすれば、S_iの尺度値は次式のようになる。(R_1 = 0)

R_i　＝　M_{zij} - M_{z0j}

ここは少しわかりにくいので、具体的に言えば、サンプルデータを用いると以下のようになる。

M_{zij} = [-1.35710151 , -0.88318025 , -0.48977349 , 0.14196916 , 0.4716406 , 0.6403424 , 0.73245961 , 0.74364347]
M_{z1j} = -1.35710151
R_1 = -1.35710151 - (-1.35710151) = 0
R_2 = -0.88318025 - (-1.35710151) = 0.47392
R_3 = -0.48977349 - (-1.35710151) = 0.86733
R_4 = 0.14196916 - (-1.35710151) = 1.49907067
R_5 = 0.4716406 - (-1.35710151) = 1.82874211
R_6 = 0.6403424 - (-1.35710151) = 1.99744391
R_7 = 0.73245961 - (-1.35710151) = 2.08956112
R_8 = 0.74364347 - (-1.35710151) = 2.10074498

ここから、尺度値R_iを使って刺激データに対する期待される尺度距離\hat Zを導き出す。

計算例
表6.1

i>j A_j B_j C_j D_j E_j H_j F_j
B_i 0.47392
C_i 0.86733 0.86733 - 0.47392 = 0.39340676
D_i 1.49907067 1.49907067 - 0.47392 = 1.02514941 1.02514941 - 0.39340676 = 0.63174265
E_i 1.82874211 1.82874211 - 0.47392 = 1.35482085 1.35482085 - 0.39340676 = 0.96141409 0.96141409 - 0.63174265 = 0.32967144
H_i 1.99744391 1.99744391 - 0.47392 = 1.52352265 1.52352265 - 0.39340676 = 1.13011589 1.13011589 - 0.63174265 = 0.49837324 0.49837324 - 0.32967144 = 0.1687018
F_i 2.08956112 2.08956112 - 0.47392 = 1.61563986 1.61563986 - 0.39340676 = 1.2222331 1.2222331 - 0.63174265 = 0.59049045 0.59049045 - 0.32967144 = 0.26081901 0.26081901 - 0.1687018 = 0.09211721
G_i 2.10074498 2.10074498 - 0.47392 = 1.62682372 1.62682372 - 0.39340676 = 1.23341696 1.23341696 - 0.63174265 = 0.60167431 0.60167431 - 0.32967144 = 0.27200287 0.27200287 - 0.1687018 = 0.10330107 0.10330107 - 0.09211721 = 0.01118386

求めた期待尺度距離\hat Z (行列\hat Z_5)

表6.2

j<i A_i B_i C_i D_i E_i H_i F_i
B_j 0.47392
C_j 0.86733 0.39340676
D_j 1.49907067 1.02514941 0.63174265
E_j 1.82874211 1.35482085 0.96141409 0.32967144
H_j 1.99744391 1.52352265 1.13011589 0.49837324 0.1687018
F_j 2.08956112 1.61563986 1.2222331 0.59049045 0.26081901 0.09211721
G_j 2.10074498 1.62682372 1.23341696 0.60167431 0.27200287 0.10330107 0.01118386

表6.2より期待比率\hat pをケースVの過程もとで、正規曲線表を使用して求める。(行列\hat P_5)

表7

i>j A_j B_j C_j D_j E_j H_j F_j
B_i 0.682221971
C_i 0.807118846 0.652990461
D_i 0.93307235 0.847353653 0.736222473
E_i 0.966280875 0.912262679 0.831827999 0.629175882
H_i 0.977111509 0.936186001 0.870786303 0.690889503 0.566984399
F_i 0.981671379 0.9469139 0.889190255 0.722569057 0.602883958 0.536697543
G_i 0.982168321 0.948112729 0.89128988 0.726304523 0.607190094 0.541137987 0.504461622

また表2.1のP行列を以下のように取り出す。

表2.2

i>j A_j B_j C_j D_j E_j H_j F_j
B_i 0.93
C_i 0.94 0.92
D_i 0.91 0.95 0.97
E_i 0.97 0.92 0.91 0.85
H_i 0.92 0.87 0.83 0.72 0.67
F_i 0.95 0.92 0.93 0.92 0.83 0.34
G_i 0.94 0.94 0.89 0.76 0.69 0.42 0.73

また、非常に大きな標本だった場合でも、同様に比率の標本の分散が正規の型では無いような範囲である場合は、これらの比率のほとんどは1.00に近い値ということになる。
標本分布が正規型をなす、統計量\theta , \hat\thetaに、各比率を変換する必要が出てくる。
つまり比率の実測値と期待値のばらつきを見るため、実測値(P)と期待値(\hat P_5)の比率を逆正弦変換した値\thetaと\hat\thetaにする必要があるということ。
変換式は、下記の通り。[1:1][2]

\theta = \frac{180 \cdot{} \arcsin \sqrt{P}}{\pi}
\hat\theta = \frac{180 \cdot{} \arcsin \sqrt{\hat P_5}}{\pi}

行列\hat P_5と行列Pを用いて以下の表8,9のように変換する。

\hat\thetaの表8

表8

i>j A_j B_j C_j D_j E_j H_j F_j
B_i 55.68668298
C_i 63.94827524 53.90859414
D_i 75.00683661 67.00194891 59.09656034
E_i 79.41886624 72.77019969 65.78964981 52.48614205
H_i 81.29834183 75.3677312 68.93277512 56.22196098 48.84949733
F_i 82.21921875 76.67912952 70.55626733 58.21608778 50.93722791 47.10450664
G_i 82.32606745 76.83313619 70.7486932 58.45560064 51.18958341 47.35970041 45.25563548

\thetaの表9

表9

i>j A_j B_j C_j D_j E_j H_j F_j
B_i 74.65829145
C_i 75.82118171 73.57005981
D_i 72.54239688 77.07903362 80.02577821
E_i 80.02577821 73.57005981 72.54239688 67.213502
H_i 73.57005981 68.86570779 65.6499364 58.05194057 54.93843704
F_i 77.07903362 73.57005981 74.65829145 73.57005981 65.6499364 35.66853756
G_i 75.82118171 75.82118171 70.6302877 60.66612575 56.16684133 40.39655189 58.69355375

\theta_{ij}と\hat \theta_{ij} の差の平方和

表10

i>j A_j B_j C_j D_j E_j H_j F_j
B_j 359.9219278
C_j 140.9659082 386.5732324
D_j 6.073463195 101.5476362 438.0321606
E_j 0.368342132 0.639776213 45.59959296 216.895131
H_j 59.72634302 42.27630845 10.77703009 3.348825285 37.07518673
F_j 26.42150316 9.666314452 16.82660188 235.744457 216.463791 130.7813888
G_j 42.3135385 1.024051865 0.014019862 4.886421275 24.77309642 48.48543735 180.5776475

カイ二乗値を求める。[1:2][2:1]
Nは刺激対ごとの判断数で、各対を判断した評価者の人数ということになる、今回は100。

\chi^2 = \frac{N \sum (\theta - \hat \theta)^2}{821}

n は、試料の数

f = \frac{(n - 1)(n - 2)}{2}

カイ二乗分布表[3]より、実測値と期待値の間には、統計的に有意な差がない場合に内的整合性が確認されたといえる。

Python3を用いて内的整合性の検定を行う。

環境

  • Python 3.9.1 (default, Dec 24 2020, 16:23:16)
  • macOS Big Sur ver11.1（このOSで動作を確認）

流れ

  • まず一意性の検定、一致性の検定を必要に応じてあらかじめ終えておく
  • 一対比較法、サーストン法を用いて尺度値を算出
    (任意の数の試料についての尺度距離行列を作成[各要素の確率に対応する正規偏差zの値を求める])
  • Zの平均値を算出（尺度値、ソート)
  • 尺度値R_iを使って刺激データに対する期待される尺度距離\hat Zを算出
  • 期待比率\hat pをケースVの過程もとで、正規曲線表を使用して行列\hat P_5を求める
  • 行列P（行列\hat P_5ではない）を整形して取り出す
  • 統計量\theta , \hat\thetaに、各比率を変換
  • 自由度fを求める
  • カイ二乗値\chi^2を求める

一意性の検定は以前の記事vol1へ
一致性の検定、有効被験者の判定結果を集計の部分は以前の記事Vol2へ
一対比較法、サーストン法を用いて尺度値を算出する部分は以前の記事vol3へ

事前準備

まず、welcomeメッセージを出力

def welcome_mes(): # 起動時のメッセージ
    print("Welcome to 一対比較法内的整合性検定システムver3:2021,1,16")
    print("paired comparison method internal consistency software, PCIC")
    print("made by kazuya yuda.")

次にz_calculation.pyで生成した解析データを読み込む。

def import_npy():
    b = ['','','']
    return_code = subprocess.check_output(['ls','./internal_consistency_data/'])
    code = return_code.split(b"\n")
    for i in range(len(code)):
        stdout_txt = str(code[i]).replace("b","").replace('\'',"")
        if re.search("npy",stdout_txt):
            print(stdout_txt)
    print("一括で解析済みデータを読み込みますか？y,n:",end = " ")
    all_import_ans = input()
    if all_import_ans == "y":
        for i in range(len(code)):
            stdout_txt = str(code[i]).replace("b","").replace('\'',"")
            if re.search("npy",stdout_txt):
                filename = stdout_txt
                if filename == "npy":
                    print("npyファイルを指定してください。")
                    sys.exit()
                elif re.search("npy",filename):
                    pass
                else:
                    print("npyファイルを指定してください。")
                    sys.exit()
                print(filename,end=" ")
                print("を読み込みました.")
                b[i] = np.load('./internal_consistency_data/'+filename)
                print(b[i])

    return b

試料数k,被験者数nを読み込む

def get_k(data):
    print(data.shape[1])
    return data.shape[1]

def get_n():
    path = './n.txt'
    with open(path) as f:
        s = f.read()
    if s == "":
        print("一致性の検定が完了していません。")
        sys.exit()
    else:
        return int(s)

尺度値R_iを使って刺激データに対する期待される尺度距離\hat Zを算出

まず、尺度値R_iを算出する。

def calculation_R(mean,k):
    print(mean[1])
    M_z1j = float(mean[1][0])
    R = [0] * k
    for i in range(k):
        M_zij = float(mean[1][i])
        R[i] = M_zij-(M_z1j)
    print(R)
    return R

次に刺激データに対する期待される尺度距離\hat Zを算出

def calculation_hatZ(mean,R,k):
    print(R)
    mapping = k-1
    hatZ_array = np.zeros((mapping,mapping))
    np.place(hatZ_array, hatZ_array == 0, 0.5)
    np.place(hatZ_array, hatZ_array == 1, 0.5)

    print(hatZ_array)
    M_z1j = float(mean[1][0])
    rock = 0
    for row in range(mapping):
        col=rock
        hatZ_array[row,col] = R[row+1]
    mapping=mapping-1
    for a in range(mapping):
        for b in range(mapping-a):
            col=a+1;row=a+b+1
            hatZ_array[row,col] = hatZ_array[row,col-1] - hatZ_array[a,a]
    print("hatZ")
    print(hatZ_array)
    return hatZ_array

期待比率\hat pをケースVの過程もとで、正規曲線表を使用して行列\hat P_5を求める

ここでは、NormalDistライブラリを使って行列を求める。

def calculation_hatP(p_array,k):
    k = k-1
    print(k)
    for a in range(k):
        for b in range(k):
            # パラメータ指定がないと、標準正規分布(mu = 0およびsigma = 1)に自動補完される。
            p_array[a,b] = NormalDist().cdf(float(p_array[a,b]))
    print(p_array)
    return p_array

行列P（行列\hat P_5ではない）を整形して取り出す

並び変えなければ、期待される尺度距離\hat Zを算出する際、誤差が生じてしまうので注意。
基本的に最後に求められる尺度値が小から大へと順番に並ぶように、刺激を配列するのが望ましい[1:3]

def swap_p(p,k,mean):
    k=k-1;a_del=np.delete(p,0,0)
    p=np.delete(a_del,0,1)
    print("origin")
    print(p)

    swap_mapdata=[]
    mapdata_origin=get_info()
    mapdata_sorted=list(mean[0])
    for i in range(k+1):
        swap_mapdata.append(mapdata_origin.index(mapdata_sorted[i]))

    col_swap=p[:,swap_mapdata]
    col_swap=col_swap[swap_mapdata,:]
    print("swaped")
    print(col_swap)

    b_del=np.delete(col_swap,0,0)
    col_swap=np.delete(b_del,k,1)
    for row in range(k):
        for col in range(k-1-row):
            col=row+col+1
            col_swap[row,col]=0
    print("result")
    print(col_swap)
    return col_swap

統計量\theta , \hat\thetaに、各比率を変換

np.nditerを使いイテレーションしている。mathライブラリを使って一度に行列計算は厳しいので、要素を取り出す必要がある。

def calculation_arcsin_P(array):
    print("def calculation_arcsin_P(array):")
    it = np.nditer(array, flags=['multi_index'])
    for x in it:
        array[it.multi_index[0],it.multi_index[1]] = math.sqrt(array[it.multi_index[0],it.multi_index[1]])
        array[it.multi_index[0],it.multi_index[1]] = math.asin(array[it.multi_index[0],it.multi_index[1]])
        array[it.multi_index[0],it.multi_index[1]] = math.degrees(array[it.multi_index[0],it.multi_index[1]])
    print(array)
    return array

def calculation_arcsin_hatP(array,k):
    print("def calculation_arcsin_hatP(array):")
    it = np.nditer(array, flags=['multi_index'])
    for x in it:
        array[it.multi_index[0],it.multi_index[1]] = math.sqrt(array[it.multi_index[0],it.multi_index[1]])
        array[it.multi_index[0],it.multi_index[1]] = math.asin(array[it.multi_index[0],it.multi_index[1]])
        array[it.multi_index[0],it.multi_index[1]] = math.degrees(array[it.multi_index[0],it.multi_index[1]])
    k=k-1
    for row in range(k):
        for col in range(k-1-row):
            col=row+col+1
            array[row,col]=0
    print(array)
    return array

統計量\theta , \hat\thetaの差の平方和を算出

def calculation_RSS(p_array,hatp_array):
    array =(p_array-hatp_array)**2
    print(array)
    return np.sum((p_array - hatp_array)**2)

自由度f,カイ二乗値\chi^2を求める

mainで各種サブルーチンを呼び出し、main内で値を算出

def main():
    welcome_mes()
    data=import_npy()
    mean=data[2]
    k=get_k(mean)
    n=get_n()
    R=calculation_R(mean,k)
    hatZ=calculation_hatZ(mean,R,k)
    p=swap_p(data[0],k,mean)
    hatP=calculation_hatP(hatZ,k)
    rss = calculation_RSS(calculation_arcsin_P(p),calculation_arcsin_hatP(hatP,k))
    print("カイ二乗値",end=" ")
    chi_2 = n * rss / 821
    print(chi_2)
    print("自由度f",end=" ")
    f = float((k-1) * (k-2)/2)
    print(f)

カイ二乗分布表[3:1]より、実測値と期待値の間には、統計的に有意な差がなければ、内的整合性が確認されたといえる。

以上で一対比較法シリーズを一旦完結とする。
余裕があれば他の解法ケースについても書き加えていきたいと思う。

Githubに上がっているものはまだまだ完成度が低いですが、少しずつ改善していく予定です。
https://github.com/Iovesophy/paired_comparison_method
(お気づきのことがあれば、pullrequest,issueなどいただけると幸いです。)

脚注

  1. 「精神測定法」 (1959年) J.P.ギルホード (著), 秋重 義治 (翻訳) ↩︎ ↩︎ ↩︎ ↩︎

  2. 「製品開発に役立つ感性・官能評価データ解析　ーRを利用してー」　,市原茂,梶谷哲也,小松原良平　（株式会社メディア・アイ） ↩︎ ↩︎

  3. 「統計の基礎と実際ー保健・臨床・家政・栄養学のためのー」 水野哲夫 （光生館） ↩︎ ↩︎ ↩︎

GitHubで編集を提案

Discussion