一対比較法、サーストン法を用いて尺度化,プロットまで（Python3+Matplotlib）

一対比較法、サーストン法を用いて尺度化,プロットまで

Python3を利用して一対比較法シリーズ第三弾です。
前回は一対比較法における一致性の検定を行いました。

https://zenn.dev/_kazuya/articles/66ef65022407ef

今回は一致性の検定を終え、サーストン法による各刺激の尺度次元上への位置付けを行う。

Github↓
https://github.com/Iovesophy/paired_comparison_method/blob/main/README.md

(Plotの様子)

サーストン法とは

※サーストンは人名(Louis Leon Thurstone.1887-1955)^[1]
一対比較法において、各刺激の尺度次元上への位置付けを行うことができる。
一般的にサーストン法を用いることが多い。
すなわち、それぞれの対象に対してより程度が大きいと判断した観察者数を全観察者数で割った値を標準得点化する方法である。この方法を用いて標準得点化することで、尺度上での刺激間の相対的な関係がわかりやすくなる。^[2]

※上図はイメージです

手法について
Louis Leon Thurstoneは、人の判断の内容、反応強度には確率的変動があるという考えに基づき、一対比較法を用いた測定の結果から判断の対象に対する尺度を与える方法を提案している。^{[1:1][2:1][3]}
ここで、一例として、試料$i$と試料$j$を一対比較する場合を考える。
特定の対象試料$i$によって感じる評価者の感覚の主観的大きさは、ばらつきを持つ。
つまり、判断を繰り返すたびに、それに対応している、反応強度$q_i$は
確率変動し(ばらつき)、期待値$\mu_i$,分散$\sigma_i^2$の正規分布になると考える。

同様に、試料$j$によって感じる評価者の感覚の主観的大きさは、反応強度$q_j$
は確率変動し、期待値$\mu_j$,分散$\sigma_j^2$の正規分布に従うと考えられる。
仮にそれぞれの反応強度を$q_i,q_j$とし,$q_i-q_j > 0$の時には、

$$試料_i > 試料_j$$

$$確率 p_{ij}=Pr(q_i-q_j > 0)$$

このように記述できる。

※上図はイメージです
試料iの分布は、試料jの分布より、評価者の感覚の主観的大きさの尺度上で、右に位置する。
しかし、試料iの分布、試料jの分布が完全には離れていなくて、重なっていた場合に問題があり、
期待値$\mu_i$の方が高くても、iとjを比較した時に、
常時、試料$i$の方が評価者の感覚の主観的大きさが大きく感じるということはないということがある。
つまり、それぞれ（$\mu_i,\mu_j$）は、正規分布にしたがい、ばらつきが有る場合と無い場合の2パターンが存在する。

両者を比較する時に、評価者がiの方がjよりも大きいと感じるのは、両者の評価者の感覚の主観的大きさの差異、($q_i-q_j$)が0よりも大きい値を取った時ということ。
さて、試料iと試料jに対する主観的な大きさの差の値を考えた時、その値は、すでに述べたように正の値を取ることもあれば、負の値を取ることもあり、その主観的な大きさの差の分布は、
つまり、$q_i - q_j$の分布は期待値、

$$\mu_i - \mu_j$$

それから、分散、

$$\sigma_i^2 + \sigma_j^2 - 2\gamma_{ij} \sigma_i \sigma_j$$

の正規分布となる。($\gamma_{ij}$は$q_iとq_j$との相関係数)
ちなみに、2つの質的変数の間にどのような関係があるかを調べたい場合は、共分散や相関係数といった指標が利用できる。
共分散は、平均$\overline{i},\overline{j}$からの偏差の積の平均であり、
iとjの共分散$S_{ij}$は次の式で求まる。

$$S_{ij} = \frac{1}{n} \sum_{a=1}^{n}(i_a - \overline{i})(j_a - \overline{j})$$

標準偏差算出における$\gamma_{ij}$は相関係数であり、共分散を２変数各々の標準偏差で割ることで求められる。

$$S_i = \sqrt{\frac{1}{n} \sum_{a=1}^{n}(i_a - \overline{i})^2}$$

$$S_j = \sqrt{\frac{1}{n} \sum_{a=1}^{n}(j_a - \overline{j})^2}$$

$$\gamma_{ij} = \frac{S_{ij}}{S_i \cdot{} S_j}$$

ここで、標準得点$z_{ij}$については以下の式で成り立つ。

$$z_{ij} = \frac{\mu_i-\mu_j}{\sqrt{\sigma_i^2 + \sigma_j^2 - 2\gamma_{ij} \sigma_i \sigma_j}}$$

※ただし、確率$p_{ij}$
が1,0の場合は求められない点に注意
Louis Leon Thurstoneは、5つのケースを考えており、^[2:2][4]

case$I$

• 一人の観察者の反復判断

case$I\hspace{-.1em}I$

• 多数の観察者が各刺激対に一回ずつ判断

case$I\hspace{-.1em}I\hspace{-.1em}I$

• $\gamma_{ij}$が0でどの刺激ついに対する回答にも相関がない場合

case$I\hspace{-.1em}V$

• $\gamma_{ij}$が0で、$\sigma_i,\sigma_j$がほぼ等しい場合

case$V$

• 最も実用的で一般的に用いられているのは、case$V$
• case$V$では、$\sigma_i$と$\sigma_j$は同じ値を取り、$\gamma_{ij}=0$と仮定する。

今回はcase$V$を使い尺度値を算出する。

case$V$を用いると以下のように簡略化される。

$$z_{ij} = \frac{\mu_i-\mu_j}{\sqrt{\sigma_i^2 + \sigma_j^2 - 2\gamma_{ij} \sigma_i \sigma_j}}$$

$$\sigma = \sigma_i = \sigma_j$$

$$\gamma_{ij}=0$$

$$z_{ij} = \frac{\mu_i-\mu_j}{\sqrt{2\sigma^2}}$$

$$\mu_i-\mu_j = z_{ij}\sigma\sqrt{2}$$

ここで、$\sigma\sqrt{2}$を尺度単位とすると以下のようになる。

$$\mu_i - \mu_j = z_{ij}$$

比較の結果から得られる各対象の尺度値$\hat{\mu_i}$は以下の式で求める。($k$は試料数)

$$\hat{\mu_i} = \frac{\sum_{i=1,i\not=j}^k z_{ij}}{k}$$

尺度値$\hat{\mu_i}$は標準化されているので、全対象の平均尺度値は0となるはずで、0から正側に離れれば離れるほど比較の際に選ばれやすいことを意味する。

サンプルデータを用いてみる

得点表（表１）から考えていく。
サンプルとしてn=100で考えている。

表1

i>j	$A_j$	$B_j$	$C_j$	$D_j$	$E_j$	$F_j$	$G_j$	$H_j$	合計
$A_i$	-	7	6	9	3	5	6	8	44
$B_i$	93	-	8	5	8	8	6	13	141
$C_i$	94	92	-	3	9	7	11	17	233
$D_i$	91	95	97	-	15	8	24	28	358
$E_i$	97	92	91	85	-	17	31	33	446
$F_i$	95	92	93	92	83	-	27	34	516
$G_i$	94	94	89	76	69	73	-	42	537
$H_i$	92	87	83	72	67	66	58	-	525

表1から比率行列Pを求める

表２

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
$F_i$	5.85967691
$G_i$	5.94914772
$H_i$	5.12273923

表4からZの平均値を算出（尺度値）

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

Python3を用いてサーストン法caseVを用いて尺度化,プロットを行う

環境

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

流れ

• まず一意性の検定、一致性の検定を必要に応じてあらかじめ終えておく
• 有効被験者の判定結果を集計
• 任意の数の試料に対する判断の確率行列を作成（比率行列、P行列、Proportion matrix)
• 任意の数の試料についての尺度距離行列を作成（各要素の確率に対応する正規偏差zの値を求める）
• 尺度値$\hat{\mu_i}$をプロットする

一意性の検定は以前の記事vol1へ
一致性の検定、有効被験者の判定結果を集計の部分は以前の記事Vol2へ

今回使用するライブラリのインポート

import pandas as pd
import csv
import subprocess
import re
import numpy as np
import sys
import datetime
from statistics import NormalDist
import matplotlib.pyplot as plt
import pylab
from IPython import get_ipython
ipy = get_ipython()
if ipy is not None:
    ipy.run_line_magic('matplotlib', 'inline')

任意の数の試料に対する判断の確率行列を作成（比率行列、P行列、Proportion matrix)
記事Vol2のコードを少し改変し、一致性が確認できた被験者のデータの行列の和をcsvに保存
データの行列の和はサブルーチン（import_csv）のmaterial[0]に格納されている。

consistency_vote_aggregatae_and_calculation.py
    now = datetime.datetime.now()
    print("プロット用ファイルを作成しますか？y,n:",end=" ")
    ans = input()
    if ans == "y":
        np.savetxt('./plot_data/' + "plotdata" + "_calculation_" + now.strftime('%Y%m%d_%H%M%S') + '.csv',material[0], delimiter=',', fmt='%d')
        print("解析用ファイルが保存されました。")
    else:
        print("解析用ファイルは保存されませんでした。")

    return material
    

ここで生成された、plotdata_calculation~.csvを元に任意の数の試料に対する判断の確率行列を作成する（比率行列、P行列、Proportion matrix)
といっても、有効被験者数nでplotdata_calculation~.csvの比率を計算する。
それぞれの得票は最大で有効被験者数nとなるので、
もし任意の試料に対する得票数が3であれば、

$$P = \frac{3}{n}$$

となる。

※標準得点$z_{ij}$について確率$p_{ij}$が1,0の場合は求められない点に注意すること。

welcomeメッセージ
def welcome_mes(): # 起動時のメッセージ
    print("Welcome to 一対比較法z計算,尺度plotシステムver3:2021,1,15")
    print("paired comparison method z calculation and plot software, PCZPloter")
    print("made by kazuya yuda.")

また、一致性の検定時に有効被験者数nを計上しているが、その値をn.txtから読み取る。

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

試料数ロード

def get_k(info):
    k = int(str(info[2]).replace("[","").replace("]","").replace("\'",""))
    return k
    

比率行列Pに変換、この際、0は0.5に1も念の為0.5に置き換える。
ただし、注意しなければならないのが、観測値が1もしくは0の場合は本来Z算出ができないことに注意すること。

def conversion_p(array,n):
    conversion_array = array/n
    np.place(conversion_array, conversion_array == 0, 0.5)
    np.place(conversion_array, conversion_array == 1, 0.5)
    return conversion_array

任意の数の試料についての尺度距離行列を作成（各要素の確率に対応する正規偏差zの値を求める）
ここは、Z算出の部分、正規曲線表から算出するのが一般的で、ここではライブラリを利用して算出（NormalDist）

def conversion_z(array):
    it = np.nditer(array, flags=['multi_index'])
    for x in it:
        # パラメータ指定がないと、標準正規分布(mu = 0およびsigma = 1)に自動補完される。
        #print(array[it.multi_index[0],it.multi_index[1]])
        #print(NormalDist().inv_cdf(array[it.multi_index[0],it.multi_index[1]]))
        array[it.multi_index[0],it.multi_index[1]] = NormalDist().inv_cdf(array[it.multi_index[0],it.multi_index[1]])
    print(array)
    return array 
    

ここは、行和を算出している。

def sum_z(conversion_z_array):
    # 行和を計上し、平均を求め尺度値を算出する
    print(conversion_z_array)
    print(np.sum(conversion_z_array, axis=1))
    return np.sum(conversion_z_array, axis=1)
    

ここは、平均を算出している。

def mean_z(conversion_z_array):
    a_del = np.delete(conversion_z_array, 0, 0)
    b_del = np.delete(a_del, 0, 1)
    # 行和を計上し、平均を求め尺度値を算出する
    print(np.mean(b_del, axis=1))
    return np.mean(b_del, axis=1)
 

最後に尺度値$\hat{\mu_i}$を使ってプロットする。

def plot_scale(k,labellist,mean_z_val):
    #配列を生成
    labellist = labellist[0].replace("[","").replace("]","").replace("\'","").replace(" ","").split(",")
    combine = sorted(zip(mean_z_val,labellist))
    print(combine)
    p,labellist = zip(*combine)
    p = list(p)
    labellist = list(labellist)

    y = [0]*k #y=0
    print(mean_z_val,labellist,p,y,k)

    #数直線
    fig,ax=plt.subplots(figsize=(15,15)) #画像サイズ
    fig.set_figheight(2.5) #高さ調整
    ax.tick_params(labelbottom=True, bottom=False) #x
    ax.tick_params(labelleft=False, left=False) #y
    ax.set_title('Psychological Scale') #タイトル

    #数値表示
    for i in range(int(k/2)):
        print(i)
        view_v = '{:.2f}'.format(p[2*i])
        ax.annotate(view_v,xy=(p[2*i],y[2*i]),xytext=(10, 20),textcoords='offset points',arrowprops=dict(arrowstyle="->"))

    for i in range(int(k/2)):
        print(i)
        view_v = '{:.2f}'.format(p[2*i+1])
        ax.annotate(view_v,xy=(p[2*i+1],y[2*i+1]),xytext=(10,-40),textcoords='offset points',arrowprops=dict(arrowstyle="->"))

    xmin, xmax= round(p[0])-.5,round(p[k-1])+.5 #数直線の最小値・最大値
    plt.tight_layout() #グラフの自動調整

    plt.subplots_adjust(left=0, right=1, bottom=0.2, top=0.8) #微調整

    markerlist = ['.','o','d','X','s','v','*','x','p','P','2']
    colorlist = ["r", "g", "b", "c", "m", "y", "#984ea3", "orange"]
    c_set = 0;m_set = 0
    for i_v in range(k):
        if m_set >= 10:
            m_set = 0
        elif c_set >= 8:
            c_set = 0
        else:
            pass
        plt.scatter(p[i_v],y[i_v], c=colorlist[c_set], s=100,marker=markerlist[m_set], label=labellist[i_v]) #size100
        c_set += 1;m_set += 1

    plt.hlines(y=0,xmin=xmin,xmax=xmax,colors='k',lw=0.8) #横軸
    #plt.vlines(x=np.arange(xmin,xmax+1,1),ymin=-0.2,ymax=0.2,colors='b') #目盛り線大
    #x=np.arange(xmin,xmax+1,1)
    #print(x)
    plt.vlines(x=np.arange(xmin,xmax,0.1),ymin=-0.1,ymax=0.1,colors='k',lw=0.3) #目盛り線小
    #x=np.arange(xmin,xmax,0.1)
    #print(x)

    plt.legend(loc='lower left') #ラベル
    line_width=0.1 #目盛り数値の刻み幅
    plt.xticks(np.arange(xmin,xmax+line_width,line_width)) #目盛り数値
    pylab.box(False) #枠を消す


    plt.show() #表示
    

mainで各種サブルーチンの呼び出し

def main():
    welcome_mes()
    material_import = import_csv(get_n())
    k = get_k(material_import[1])
    labellist = material_import[1][3]
    conversion_array = material_import[0]
    conversion_z_array = conversion_z(conversion_array)
    sum_z_val = sum_z(conversion_z_array)
    mean_z_val = mean_z(conversion_z_array)
    plot_scale(k,labellist,mean_z_val)
    

ちなみに、以前数直線グラフのプロット方法の記事を上げているので、参考までに上げておく。
https://qiita.com/_kazuya/items/9910f85909f46f319689

実行

このような形でプロットできる。

レイアウトは適宜見やすいように調節しましょう。
ここまで、できたらあと少しです！
次回は求めた尺度値から内的整合性について検定します。
https://zenn.dev/_kazuya/articles/0ad3de87944785

※間違い、誤植等、発見されましたらご指摘、ご指南いただけると幸いです。

脚注

1. 「Louis Leon Thurstone」 , (wikipedia) https://en.wikipedia.org/wiki/Louis_Leon_Thurstone

2. 「基礎心理学実験法ハンドブック」 ,日本基礎心理学会 （責任編集：坂上貴之,河原純一郎,木村英司,三浦佳世,行場次郎,石金浩史）

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

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