潜在クラス項目反応理論／Stan

潜在クラス項目反応理論（単調性を仮定しない場合）用プログラム

Stanスクリプト例

現代テスト理論の項目反応理論（Item Response Model: IRT）では、正答確率を特性値の関数として

（１）

と表す。ここで、特性値は連続変数である。しかし、Piagetの発達段階理論に見られるように、心理学的概念あるいは特性は心理特性の段階を表すカテゴリ変数の場合もある。このとき、正答確率は心理特性の段階の関数として

（２）

と表される。モデル（２）では、発達段階理論のように段階（Stage）に順序関係があり、確率（２）はその単調関数になっている。しかし一般的には、順序関係のある段階ではなく単に潜在クラスとして離散状態を設定するモデル

（３）　　　　

が適切な場合も考えられる。モデル（３）では反応（Response）の確率となっているが、これは項目に対する回答が「正答／誤答」だけではなく、「当てはまる／当てはまらない」の２値反応を対象としていることを明示するためである。また、モデル（３）における反応は２値（例えば、０と１）である必要はないが、ここでは簡単のため２値の場合のモデルについて考える。２値でないデータに対して２値モデルを適用するときは、データを２値に変換すればよい。例えば４カテゴリ、「当てはまらない」「どちらかと言えば当てはまらない」「どちらかと言えば当てはまる」「当てはまる」、であれば、「当てはまらない」と「どちらかと言えば当てはまらない」を反応０とし、「どちらかと言えば当てはまる」と「当てはまる」を反応１にまとめればよい。クラスは潜在クラスである必要はないが、ここでは潜在クラス推定のための分析について考えている。心理学理論において連続変量ではなく潜在クラスが想定されているものとして、統合失調症に対するtaxon（Meehl, 2004）、感情の離散モデル（Levenson, 2014）などを上げることができる。

モデル（３）の確率モデルにベイズ分析を適用すると、データに対して事後分布が与えられるが、モデル（３）のみではMCMCサンプリングにおいてラベルスイッチングの起こる危険がある。ラベルスイッチングを避けるため、以下の制約を置く。

項目のクラスにおける正反応（「正答」、「当てはまる」など）をで表すとき、モデル（３）は

（４）　　　　

と書ける。項目がクラスを表すものであるときは、他のクラスに対する確率より高いとする。すなわち、

（５）　　　　

である。

仮定（５）の下で分析を行うが、本ウェブサイトのプログラムでは、各項目がどのクラスに対応するかまとめられているものとする。関連研究からの知見によって、どのような潜在クラスが想定できるか、各潜在クラスにどのような質問項目（検査項目）が対応するかについての仮説が設定されていることを考えている。しかし、このような仮説設定が難しい段階の研究の場合は、まず収集されたデータに対して、最初は予備分析として、主成分分析あるいはクラスタ分析などで項目をグループ分けしてまとめ、本ウェブサイトのプログラムで仮分析を行い、その後、その仮分析の結果を見て、例えばMCMCのサンプリングのの値が大き過ぎるとか、式（５）における差が小さいという項目を除いて項目群を選び、再度、本ウェブサイトのプログラムで本分析を行えばよい。なお、項目の主成分分析の例を、本ウェブサイトの最後に示した。項目のグループ分けはクラスタ分析によって行うこともできる。クラスタ分析の例も主成分分析の例に続けて示した。

本ウェブサイトの分析によりの推定値が得られると、回答パターンが与えられた回答者の属するクラスは、ベイズ分析における事後分布

により推定できる。ここで、は回答者の属する集団におけるclassの事前分布である。

データ例１に、潜在クラスが２個、各カラスに属する項目が５個である例（人工データ）を示す。各行が一人分のデータである。

データ例１

[1, 1, 1, 1, 1, 0, 0, 0, 0, 0]

[0, 0, 0, 0, 0, 1, 1, 1, 1, 1]

正反応が１、負反応が０で表されている。最初の５項目がクラス１に対応するものであり、後半の５項目がクラス２に対応している。前半の１０人がクラス１に属しており、後半の１０人がクラス２に属していることを想定した人工データである。データの分析においては、各項目がどのクラスに対応するものであるのかが分かっていて、クラスごとに項目がまとめて並べられておればよい。各回答者の属するクラスは分析結果として得られるので、回答者の行の並びの順序は任意である。データ例１のデータを分析するスクリプト例をリスト１に示す。なお、発達段階に応じて正答確率が高くなるような単調関係を設定した分析例は別のウェブサイトで扱っている。

潜在クラスモデルでは、各回答者の属するクラスはカテゴリ変数であるが、Stanではパラメータとしてカテゴリ変数を設定することはできない。これに対応するため、パラメータの確率の周辺化を行い、潜在クラスを直接パラメータとして扱うことを避ける。次の周辺化を行い、回答者の正答確率を設定する。

リスト１のスクリプトでは、確率をsimplex[NCls]型のN次元配列thetaで表している。Nは回答者の数である。このとき、

である。NCls＝クラス数であることに注意。

リスト１のスクリプトは、データ例１を生成して、その分析を行うデモンストレーションプログラムであるが、より一般的に、ファイルでデータを用意して、クラス数の３以上の場合も扱えるスクリプト例は、このウェブサイトの後半で扱う。

なお、PyStan（StanのPython用ライブラリ）のWindowsでの利用法の説明を、このウェブサイトで行っている。

リスト１　モデル（５）に基づくスクリプト例

"""

Yasuharu Okamoto, 2019.04

"""

import pystan

import pickle

import matplotlib.pyplot as plt

import seaborn as sb

import numpy as np

M1 = 5

M2 = 5

X = ([[1] * M1 + [0] * M2] * 10) + ([[0] * M1 + [1] * M2] * 10)

print('X =')

for v in X:

print(v)

N = len(X)

print('N = ', N)

M = M1 + M2

print('M1 = ', M1, ' M2 = ', M2)

print('M = ', M)

NCls = 2

print('NCls = ', NCls)

stan_data = {'N': N, 'M1': M1, 'M2': M2, 'M': M, 'NCls': NCls, 'X': X}

stan_code = """

// Binary Responses for Two Latent Classes

data {

int N; // Number of subjects

int M1; // Number of items for latent class 1

int M2; // Number of items for latent class 2

int M; // Number of items M1 + M2

int NCls; // Number of classes

int X[N, M]; // Responses: 0 -> failure, 1 -> success

}

parameters {

simplex[NCls] theta[N];

real<lower = 0.0, upper = 1.0> pp[M, NCls];

}

transformed parameters {

real p[M, NCls];

real p_person[N, M];

for (i in 1:M1){

for (stg in 1:NCls){

if (stg == 1){

p[i][stg] = pp[i][stg];

}

else{

p[i][stg] = pp[i][stg] * pp[i][1];

}

for (i in (M1+1):(M1+M2)){

for (stg in 1:NCls){

if (stg == 2){

p[i][stg] = pp[i][stg];

}

else{

p[i][stg] = pp[i][stg] * pp[i][2];

}

for (j in 1:N) {

for (i in 1:M){

p_person[j][i] = 0.0;

for (stg in 1:NCls){

p_person[j][i] = p_person[j][i] + p[i][stg] * theta[j][stg];

}

model {

for (i in 1:M){

for (stg in 1:NCls){

pp[i][1] ~ uniform(0.0, 1.0);

}

for (j in 1:N) {

for (i in 1:M){

X[j][i] ~ bernoulli(p_person[j][i]);

}

"""

sm = pystan.StanModel(model_code = stan_code)

def f_init():

theta = []

for j in range(N):

theta.append([])

for i in range(NCls):

theta[j].append(1.0/NCls)

pp = []

for i in range(M):

pp.append([])

for s in range(NCls):

pp[i].append(0.8)

return dict({'theta': theta, 'pp': pp})

fit = sm.sampling(data = stan_data, n_jobs = 1, init = f_init)

print(fit)

p_samples = fit['p']

p_item_class = np.empty((M, NCls))

for i in range(M):

for s in range(NCls):

p_item_class[i][s] = p_samples[:,i,s].mean()

print('\n\nP(correct|item,class)\n')

print('{0:>7s}{1:>9s}{2:>9s}'.format(' ', 'Class-0', 'Class-1'))

for i in range(M):

print('Item-{0} '.format(i), end = '')

for s in range(NCls):

print('{0:>9.5f}'.format(p_item_class[i][s]), end = '')

print(' ')

theta = fit['theta']

p_class = np.empty((N, NCls))

for j in range(N):

for s in range(NCls):

p_class[j][s] = theta[:,j,s].mean()

print('\n\nP(class|person)\n')

print('{0:>10s}{1:>9s}{2:>9s}'.format(' ', 'class-0', 'class-1'))

for j in range(N):

print('Person-{0: <3d}'.format(j), end = '')

for s in range(NCls):

print('{0:>9.5f}'.format(p_class[j][s]), end = '')

print(' ')

リスト１のスクリプトファイルは、圧縮ファイルSimpleSample.zipとして用意した。スクリプトの実行は、Windowsの場合は、コンソール端末ソフト（例えば、Windows PowerShell）で行うと、MCMCサンプリングの経過が表示される。Pythonのインストールと使い方は、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan、PyMC）――」丸善出版＞の「序章　Pythonの準備と使い方」で説明している。

コンソール端末（Windows PowerShell）において、pythonコマンドを入力している例を図１に示す。ファイル名LClassSimpleStan.pyは、リスト１のスクリプトのファイル名である。

図１

リスト１を実行すると（PyStan 2.17.1を使用）、データ例１のデータが生成表示され、回答者数などの表示の後、Stanスクリプトのコンパイルが始まる。コード生成が終了すると、サンプリングが始まり、経過が表示される（図２）。

図２

サンプリングが終了すると、以下の結果が表示される。

P(correct|item,class)

Class-0 Class-1

Item-0 0.89943 0.09323

Item-1 0.89838 0.09316

Item-2 0.89942 0.09239

Item-3 0.89727 0.09200

Item-4 0.89788 0.09294

Item-5 0.09279 0.89942

Item-6 0.09305 0.89746

Item-7 0.09297 0.89919

Item-8 0.09182 0.89897

Item-9 0.09269 0.89889

P(class|person)

class-0 class-1

Person-0 0.90617 0.09383

Person-1 0.90550 0.09450

Person-2 0.90734 0.09266

Person-3 0.90880 0.09120

Person-4 0.90450 0.09550

Person-5 0.90750 0.09250

Person-6 0.90654 0.09346

Person-7 0.90900 0.09100

Person-8 0.90533 0.09467

Person-9 0.90736 0.09264

Person-10 0.09131 0.90869

Person-11 0.09284 0.90716

Person-12 0.09337 0.90663

Person-13 0.09325 0.90675

Person-14 0.09355 0.90645

Person-15 0.09362 0.90638

Person-16 0.09352 0.90648

Person-17 0.09327 0.90673

Person-18 0.09283 0.90717

Person-19 0.09315 0.90685

PS D:\yasuharu\www\LaCoocan\Python\misc\LClassNoOrdStan\SimpleSample>

P(correct|item,class)は、の事後分布の平均値である。項目０から項目４までのクラス０に対する確率は高いが、クラス１に対しては低い。出力では、配列の添え字が０から始まっていることに注意。

P(class|person)は、を表すパラメータthetaの事後分布の平均値である。回答者０から回答者９まではクラス０の確率が高く、回答者１０から回答者１９まではクラス１の確率が高い。

これらの結果は、スクリプト１の分析プログラムの可能性を示すものである。これを、潜在クラスが３クラス以上の場合も扱えるようにしたものを、次に説明する。このプログラムでは、データはファイルから入力する。

ファイル入力データの分析プログラム例

データをファイルから入力する場合のプログラム例を説明する。

データはCSV形式で用意する（図３）。CSV形式のファイルは、Excelで保存時のファイル拡張子を.csvとすることにより、簡単に用意できる。

図３　サンプルデータ（Data.csv）

図３では、先頭の２２人分のデータが示されている。第１欄には、回答者の識別番号が入力されている。第２欄から１５個の項目が順番に並んでいる。Excelのグリッドの第１行目は、変数の名前が入力されている。回答データは、正反応が１、負反応が０で示されている。図３のデータの場合、潜在クラス数が３個、各クラスタの項目数はすべて同じ５個として作成されている（人工データ）。最初の５項目が第１クラスタ、次の５項目が第２クラスタ、最後の５項目が第３クラスタに対応している。各項目がどのクラスに対応するものであるかは、２行目において各項目名の下に示されている。クラスの総数がK個のときは、１からKまでの整数値でクラスを表すことに注意。また、項目はクラスごとにまとめられていることにも注意。どのクラスが何番目になるかは任意である。しかし、クラスを表す整数値は、左から順番に１、・・・、Kと昇順でなければならない。これは、リスト３のプログラムが、クラスを表す整数値が左から順番に昇順であることを前提としているからである。図３のデータの回答者は、各クラスタ１０人である。項目の順序（列の順序）は、同じクラスに対応する項目ごとにまとめられている必要があるが、回答者の順序（行の順序）は任意である（所属クラスごとにまとめられている必要はない。所属クラスは、分析結果として推定される）。項目のグループを並べ替えた例は、この後で示す。

図３のCSVファイルデータを読み込み、潜在クラス項目反応モデルで分析するためのStanスクリプトをリスト２に、リスト２のStanスクリプトを使って図３のデータを分析するスクリプトをリスト３に示す。リスト２のStanスクリプトのファイルは、リスト３のファイル、図３のデータファイルとともに圧縮ファイルGeneralSample.zipとしてまとめた。

リスト２　潜在クラス項目分析Stanスクリプト（ファイル名LClass.stan）

// Binary Responses for NCls Latent Classes

// Yasuharu Okamoto, 2019.04

data {

int N; // Number of subjects

int NCls; // Number of classes

int Msub[NCls]; // Number of items in each class

int M; // Total number of items = sum(Msub)

int X[N, M]; // Responses: 0 -> failure, 1 -> success

vector[NCls] a; // Parameter of the prior distribution of the classes

}

transformed data{

int CumM[NCls];

CumM[1] = Msub[1];

for (s in 2:NCls){

CumM[s] = CumM[s-1] + Msub[s];

}

parameters {

simplex[NCls] theta[N];

real<lower = 0.0, upper = 1.0> pp[M, NCls];

}

transformed parameters {

real p[M, NCls];

real p_person[N, M];

vector[2] p_person_cat[N, M];

for (i in 1:CumM[1]){

for (s in 1:NCls){

if (s == 1){

p[i][s] = pp[i][1];

}

else{

p[i][s] = pp[i][s] * pp[i][1];

}

for (s_g in 2:NCls){

for (i in (CumM[s_g-1]+1):CumM[s_g]){

for (s in 1:NCls){

if (s == s_g){

p[i][s] = pp[i][s_g];

}

else{

p[i][s] = pp[i][s] * pp[i][s_g];

}

for (j in 1:N) {

for (i in 1:M){

p_person[j][i] = 0.0;

for (stg in 1:NCls){

p_person[j][i] = p_person[j][i] + p[i][stg] * theta[j][stg];

}

model {

for (j in 1:N){

theta[j] ~ dirichlet(a);

}

for (i in 1:M){

for (s in 1:NCls){

pp[i][s] ~ uniform(0.0, 1.0);

}

for (j in 1:N) {

for (i in 1:M){

X[j][i] ~ bernoulli(p_person[j][i]);

}

リスト２のStanスクリプトでは、クラス（thetaで表されている）の事前分布がディリクレ分布分布で与えられている。

for (j in 1:N){

theta[j] ~ dirichlet(a);

}

である。

リスト３のスクリプトでは、リスト２のStanスクリプトをファイル名LClass.stanのファイルとして読み込んでいる。

リスト３　データファイルを読み込んで分析するスクリプト例

"""

Yasuharu Okamoto, 2019.04

"""

import csv

import pystan

import pickle

import matplotlib.pyplot as plt

import seaborn as sb

import numpy as np

fin_name = input('Data file name(*.csv) = ')

with open(fin_name, 'r') as f:

data = [v for v in csv.reader(f)]

fout_name = input('Output file name(*.txt) = ')

fout = open(fout_name, 'w')

fout.write('Data file name = {}\n'.format(fin_name))

N = len(data) - 2

M = len(data[0]) - 1

ItemNames = data[0][1:]

ck_cnt = data[1][1:]

cnts = np.bincount(ck_cnt)

NCls = len(cnts) - 1

Msub = []

for s in range(NCls):

Msub.append(cnts[s+1])

X = np.empty((N, M), dtype = int)

CaseID = []

for i in range(2, N + 2):

CaseID.append(data[i][0])

for j in range(M):

X[i - 2][j] = int(data[i][j + 1])

print('X =')

fout.write('\nX =\n')

i = 0

for v in X:

print(CaseID[i], ' ',end = '')

fout.write('{} '.format(CaseID[i]))

i += 1

print(v)

fout.write('{}\n'.format(v))

N = len(X)

print('N = ', N)

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

M = len(X[0])

print('M = ', M)

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

if M != sum(Msub):

print('\nError in number of M, or M[]')

raise Exception()

print('\nNumber of items for each class')

for j in range(NCls):

print('M[{0}] = {1}'.format(j+1, Msub[j]))

fout.write('M[{0}] = {1}\n'.format(j+1, Msub[j]))

print('NCls = ', NCls)

fout.write('NCls = {}'.format(NCls))

a = [] # Parameter of the prior distribtuion of the classes

for j in range(NCls):

a.append(1.0);

stan_data = {'N': N, 'M': M, 'Msub': Msub, 'NCls': NCls, 'a': a, 'X': X}

sm = pystan.StanModel(file = 'LClass.stan')

def f_init():

pp = []

pos_pp = -1

for s in range(NCls):

for i in range(Msub[s]):

pp.append([])

pos_pp += 1

for si in range(NCls):

pp[pos_pp].append(0.8)

return dict({'pp': pp})

fit = sm.sampling(data = stan_data, pars = ['p', 'theta'],

n_jobs = 1, init = f_init)

print(fit)

p_samples = fit['p']

p_item_class = np.empty((M, NCls))

for i in range(M):

for s in range(NCls):

p_item_class[i][s] = p_samples[:,i,s].mean()

print('\n\nP(+|item,class)\n')

fout.write('\n\nP(+|item,class)\n\n')

print('{0: <9s}'.format('Item'), end = '')

for s in range(NCls):

print('{0:>9s}'.format('Class-{}'.format(s+1)), end = '')

print()

fout.write('{0: <9s}'.format('Item'))

for s in range(NCls):

fout.write('{0:>9s}'.format('Class-{}'.format(s+1)))

fout.write('\n')

for i in range(M):

print('{0: <9s}'.format(ItemNames[i]), end = '')

fout.write('{0: <9s}'.format(ItemNames[i]))

for s in range(NCls):

print('{0:>9.5f}'.format(p_item_class[i][s]), end = '')

fout.write('{0:>9.5f}'.format(p_item_class[i][s]))

print(' ')

fout.write('\n')

theta = fit['theta']

p_class = np.empty((N, NCls))

for j in range(N):

for s in range(NCls):

p_class[j][s] = theta[:,j,s].mean()

print('\n\nP(class|person)\n')

fout.write('\n\nP(class|person)\n\n')

print('{0: <10s}'.format('Person'), end = '')

fout.write('{0: <10s}'.format('Person'))

for s in range(NCls):

print('{0:>9s}'.format('class-{}'.format(s+1)), end = '')

fout.write('{0:>9s}'.format('class-{}'.format(s+1)))

print()

fout.write('\n')

for j in range(N):

print('{0: <10s}'.format(CaseID[j]), end = '')

fout.write('{0: <10s}'.format(CaseID[j]))

for s in range(NCls):

print('{0:>9.5f}'.format(p_class[j][s]), end = '')

fout.write('{0:>9.5f}'.format(p_class[j][s]))

print(' ')

fout.write('\n')

fout.close()

print('\nOutput file {} was saved.\n'.format(fout_name))

クラスの事前分布であるディリクレ分布のパラメータの値は

a = [] # Parameter of the prior distribtuion of the classes

for j in range(NCls):

a.append(1.0);

と全て１であるので、無情報事前分布である。このパラメータの値は、データに応じて設定されるものであるが、とりあえず無情報事前分布を設定している。ディリクレ分布については別のウェブサイトで説明している。

リスト３のスクリプトをWindowsパソコンで実行する場合は、コンソール端末（例えば、Windows PowerShell）から実行すると、MCMCサンプリングの途中経過が表示される。例えば、リスト３のスクリプトのファイル名がLClassFLStan.pyである場合、ファイルの置かれているフォルダに於いて

Python LClassFLStan.py

と打ち込めばよい（図４）。

図４

Pythonコマンドを打ち込んでスクリプトを実行すると、入力データファイル名と出力用テキストファイル名の設定が求められる。入力データファイルとして図３のデータのファイル名Data.csvが設定され、出力用テキストファイル名は適当なテキストファイル名が設定されている。出力ファイル名を設定すると、次に、読み込んだデータが表示さる。データの表示の最後に出力されているNとMは、回答者数と項目数である（図５）。

図５

続いて、潜在クラス数NClsとクラス別の項目数M[class]が表示される。その後、Stanスクリプトのコンパイルが始まる。コード生成終了後、サンプリングが始まり、途中経過が表示される（図６）。IDLEからの実行のときは、サンプリングの途中経過は表示されない。

図６

MCMCサンプリング終了後、以下の出力が表示される。

P(+|item,class)

Item Class-1 Class-2 Class-3

X1 0.86952 0.09489 0.09620

X2 0.87253 0.09547 0.09301

X3 0.87348 0.09447 0.09435

X4 0.87179 0.09402 0.09562

X5 0.87338 0.09461 0.09751

X6 0.09494 0.87327 0.09327

X7 0.09393 0.87164 0.09604

X8 0.09488 0.87235 0.09469

X9 0.09608 0.87150 0.09520

X10 0.09467 0.87341 0.09366

X11 0.09525 0.09372 0.87126

X12 0.09549 0.09641 0.87296

X13 0.09714 0.09498 0.87009

X14 0.09544 0.09485 0.87282

X15 0.09677 0.09655 0.87269

P(class|person)

Person class-1 class-2 class-3

1 0.81737 0.08933 0.09331

2 0.81663 0.09209 0.09128

3 0.81954 0.09049 0.08997

4 0.81623 0.09032 0.09346

5 0.81710 0.09239 0.09051

6 0.81796 0.09214 0.08990

7 0.81587 0.09173 0.09239

8 0.81728 0.09061 0.09211

9 0.81478 0.09379 0.09143

10 0.81960 0.09026 0.09015

11 0.09249 0.81435 0.09317

12 0.09304 0.81718 0.08978

13 0.09085 0.81802 0.09113

14 0.09061 0.81771 0.09167

15 0.09011 0.81844 0.09144

16 0.09192 0.81737 0.09071

17 0.09096 0.81822 0.09082

18 0.09088 0.81702 0.09210

19 0.09268 0.81772 0.08960

20 0.09245 0.81948 0.08806

21 0.08982 0.09084 0.81934

22 0.09303 0.09140 0.81557

23 0.09072 0.09118 0.81810

24 0.09077 0.09016 0.81907

25 0.09121 0.09082 0.81797

26 0.09024 0.09054 0.81922

27 0.09279 0.09125 0.81595

28 0.09316 0.09268 0.81416

29 0.09306 0.09066 0.81627

30 0.09087 0.09076 0.81838

Output file Results.txt was saved.

PS D:\yasuharu\www\LaCoocan\Python\misc\LClassNoOrdStan\GeneralSample>

まず、p[M,NCls]の事後分布の平均値がP(correct|item,class)として出力されている。X1からX5までは、クラス１での値が高く、X6からX10まではクラス２での値が高く、X11からX15まではクラス３での値が高くなっている。

次に、theta[j,s]の事後分布の平均値が、P(class|person)として出力されている。Person１からPerson10まではクラス１の値が高く、Person11からPerson20まではクラス２の値が高く、Person21からPerson30まではクラス３の値が高い。

以上の結果は、人工データに対するものであるが満足のいくものであり、実データの分析においても本ウェブサイトのプログラムが有効であることが期待される。

図３に示されているデータにおいて項目群を入れ替えた例

項目群の順序が任意であることを確認するため、図３の項目群を入れ替えたものを図３aのように用意した。

図３a

上の図３aのデータは、図３のデータにおいて項目群X1からX5を右側の最後に移動したものである。２行目のクラス番号は、１から３までの昇順に書き直されていることに注意。

図３aのデータに対してプログラムを実行すると、読み込まれたデータが以下のように表示される。

X =

1 [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]

2 [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]

3 [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]

4 [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]

5 [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]

6 [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]

7 [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]

8 [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]

9 [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]

10 [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]

11 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]

12 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]

13 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]

14 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]

15 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]

16 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]

17 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]

18 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]

19 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]

20 [1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]

21 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]

22 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]

23 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]

24 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]

25 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]

26 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]

27 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]

28 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]

29 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]

30 [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]

N = 30

M = 15

Number of items for each class

M[1] = 5

M[2] = 5

M[3] = 5

NCls = 3

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_dee74dfece458dbb4907ad69b4d4e9ad NOW.

Stanスクリプトによる分析結果は以下のようになる。

P(+|item,class)

Item Class-1 Class-2 Class-3

X6 0.87374 0.09475 0.09286

X7 0.87370 0.09330 0.09545

X8 0.87285 0.09512 0.09488

X9 0.87275 0.09633 0.09558

X10 0.87141 0.09651 0.09592

X11 0.09580 0.87050 0.09711

X12 0.09476 0.87355 0.09493

X13 0.09518 0.87203 0.09530

X14 0.09473 0.87417 0.09619

X15 0.09516 0.87168 0.09515

X1 0.09740 0.09431 0.87021

X2 0.09553 0.09566 0.87174

X3 0.09390 0.09616 0.87287

X4 0.09507 0.09535 0.87286

X5 0.09375 0.09613 0.87035

P(class|person)

Person class-1 class-2 class-3

1 0.09192 0.09188 0.81620

2 0.09182 0.09000 0.81818

3 0.09215 0.09090 0.81695

4 0.09036 0.09170 0.81795

5 0.08944 0.09298 0.81758

6 0.09280 0.09058 0.81662

7 0.09113 0.09067 0.81820

8 0.09142 0.09214 0.81644

9 0.09125 0.09201 0.81674

10 0.09130 0.09160 0.81710

11 0.81929 0.09057 0.09014

12 0.81585 0.09184 0.09231

13 0.81669 0.09157 0.09174

14 0.81813 0.09098 0.09090

15 0.81761 0.09200 0.09039

16 0.81803 0.09062 0.09135

17 0.81951 0.09112 0.08937

18 0.81571 0.09386 0.09042

19 0.81675 0.09006 0.09319

20 0.81618 0.09097 0.09286

21 0.09065 0.81800 0.09135

22 0.09305 0.81646 0.09049

23 0.09172 0.81716 0.09112

24 0.09081 0.81900 0.09019

25 0.09178 0.81588 0.09234

26 0.09115 0.81922 0.08963

27 0.09212 0.81487 0.09301

28 0.09099 0.81697 0.09204

29 0.09181 0.81500 0.09319

30 0.09030 0.81821 0.09149

図３のデータではクラス１、２、３となっているクラスが、図３aを入力データとした場合は、クラス３，１，２となっていることに注意すれば、図３のデータを入力ファイルとする場合と図３aのデータを入力ファイルとする場合が本質的に同じ結果であることが確認できる。

項目の主成分分析

上のサンプルデータにおける項目の主成分分析例を示す。プログラムは

岡本安晴「いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan、PyMC）――」丸善出版

の「第７章　主成分分析」のプログラムをダウンロードして用いた。

データ例１の項目を主成分分析すると、図７に示すパターン行列の布置を得る。主成分数を１に設定しているので、横軸と縦軸は同じ値である。

図７

出力ファイルを開くと、以下のように出力されている。

Pattern matrix for the principal components =

comp.1

X1 -1.00000

X2 -1.00000

X3 -1.00000

X4 -1.00000

X5 -1.00000

X6 1.00000

X7 1.00000

X8 1.00000

X9 1.00000

X10 1.00000

図３のデータの項目を主成分分析すると、図８のようなパターン行列の布置が得られる。

図８

出力ファイルを開くと、パターン行列は以下のようになっている。

Pattern matrix for the principal components =

comp.1 comp.2

X1 -0.00000 1.00000

X2 0.00000 1.00000

X3 -0.00000 1.00000

X4 -0.00000 1.00000

X5 -0.00000 1.00000

X6 -0.86603 -0.50000

X7 -0.86603 -0.50000

X8 -0.86603 -0.50000

X9 -0.86603 -0.50000

X10 -0.86603 -0.50000

X11 0.86603 -0.50000

X12 0.86603 -0.50000

X13 0.86603 -0.50000

X14 0.86603 -0.50000

X15 0.86603 -0.50000

図８のマップ、および上のパターン行列を見ると、項目がX1からX5まで、X6からX10まで、X11からX15までの３グループにまとまっていることが確認できる。

データ１および図３のデータは人工データであるので主成分分析の結果は明確に項目群を示すものであるが、実データの場合も同様に項目群のまとまりが得られると予想される。

項目群のクラスタ分析

　項目群のクラスタ分析によって得られた項目のクラスタが潜在クラスに対応している可能性がある。潜在クラスを見つけるために受検者のクラスタ分析を行うことも考えられるが、ここでは項目のクラスタ分析について例示する。クラスタ分析は、ライブラリscipyに用意されているものを用いる。スクリプトは以下のようである。

リストC.1　クラスタ分析スクリプト例

import csv

import numpy as np

import pandas as pd

from scipy.cluster.hierarchy import linkage

from scipy.cluster.hierarchy import dendrogram

import matplotlib.pyplot as plt

fl_nm = input('Input Data file (*.csv) = ')

with open(fl_nm, 'r') as f:

raw_data = [v for v in csv.reader(f)]

ItemLabels = []

for i in range(1, len(raw_data[0])):

ItemLabels.append(raw_data[0][i])

print(ItemLabels)

n_items = len(ItemLabels)

Person = []

Data = []

n = 0

for i in range(1, len(raw_data)):

Person.append('P-' + raw_data[i][0])

Data.append([])

for j in range(1, n_items + 1):

Data[n].append(int(raw_data[i][j]))

n += 1

df = pd.DataFrame(np.array(Data).T, columns = Person, index = ItemLabels)

print(df)

row_clusters = linkage(df.values, method = 'ward')

fig = plt.figure(figsize = (8, 8), facecolor = 'white')

row_dendr = dendrogram(row_clusters, labels = ItemLabels,

orientation = 'right')

plt.show()

データはCSV形式として用意するが、上の図３のデータの場合は、次の図C.1のように用意して、CSV形式のファイルとして保存する。

図C.1

図３のデータにおいて、所属潜在クラスに関する行を削除したものである。また、各項目は５個ずつ同じ反応パターンであるが、このとき同じ反応パターンの項目はクラスタ分析のデンドログラムにおいて分岐枝が重なり結果が見難くなるので、少しノイズを加えてある。図C.1のクラスタ分析のデンドログラムは下図のようである（図C.2）。

図C.2　Ward法による項目変数のクラスタ分析

リストC.1を実行すると図C.3のように入力データファイル名の設定が求められる。

図C.3

図C.1のデータファイル名Data_Noise.csvを設定してEnterキーを押すと、クラスタ分析が始まり、図C.２の結果が表示される。

リストC.1のスクリプトファイルと図C.1のデータファイルData_Noise.csvはまとめて圧縮ファイルClusteringItems.zipとした。ダウンロード解凍すれば、自由に利用できる。

クラスタ分析の結果は、クラスタの構成法によって分析結果が異なることが多い。異なる方法の結果を統合するために、クラスタの構成法によらず同じクラスタに分類される項目同士は関係が強いと考える。各クラスタ分析の結果において、同じクラスタに属するかどうかを表す行列、すなわち同じクラスタに属する場合は１、そうでない場合は０が置かれている行列を構成して、クラスタ分析間の平均値を求めれば、平均値が１に近いほど、そのセルに対応する項目同士は関係が強いと考えられる。この同じクラスタに属するかどうかを表す行列を求めるスクリプトも作成してみた。リストC.1のスクリプトにリストC.2のスクリプトを追加したものである。

リストC.2　項目間のクラスタ関係を表す行列の作成

ItemList = [i for i in range(n_items)]

N_Clusters = int(input('\nN_Clusters = '))

id_cluster = n_items

step = 0

while True:

c0 = row_clusters[step][0]

c1 = row_clusters[step][1]

for i in range(len(ItemList)):

if ItemList[i] == c0:

ItemList[i] = id_cluster

if ItemList[i] == c1:

ItemList[i] = id_cluster

id_cluster += 1

step += 1

if step > len(row_clusters) - N_Clusters:

break

print('\n', ItemList)

FraternityMatrix = np.zeros((n_items, n_items), dtype = int)

for i in range(n_items):

for j in range(n_items):

if i == j:

FraternityMatrix[i][i] = 1

else:

if ItemList[i] == ItemList[j]:

FraternityMatrix[i][j] = 1

FraternityMatrix[j][i] = 1

print('\n|Fraternity Matrix = ')

print('{0:>3s}'.format(' '), end = '')

for i in range(n_items):

print('{0:>3s}'.format(ItemLabels[i]), end = '')

print()

for i in range(n_items):

print('{0:>3s}'.format(ItemLabels[i]), end = '')

for j in range(n_items):

print('{0:>3d}'.format(FraternityMatrix[i][j]), end = '')

print()

リストC.1とリストC.2のスクリプトを合わせたプログラムは圧縮ファイルClusterFraternity.zipの中に用意した。

このプログラムを実行すると、リストC.1の実行と同じように進み、図C.2のグラフが表示される。

図C.2のウィンドウを閉じると、「N_Clusters =」とクラスタ数の入力が求められる。図C.2の樹状図では３つのクラスタに分かれているので、それに合わせて図C.4では３を設定している。他のクラスタ数でもよい。クラスタ数を設定すると、計算が始まり、図C.4に示すように各項目がどの項目と同じクラスタに属するかを表す行列が出力される。

図C.4

複数のクラスタ分析を行ったとき、項目のクラスタ所属関係を表す行列（図C.4）をクラスタ分析ごとに作成して、それらの平均値を求めれば、クラスタ構成法によらず関係の強い項目を選んでクラスタを構成できることが期待される。

引用文献

Levenson, R. W. (2014). What is the underlying structure of emotions: An argument for discrete emotions. In S. N.-Hoeksema, B. L. Fredrickson, G. R. Loftus, & C. Lutz (2014). Atkinson & Hilgard”s Introduction to Psychology, 16th edition. Cengage Learning. pp. 406-407.

Meehl, P. E. (2004). What’s in a taxon? Journal of Abnormal Psychology, 2004, 113, 39-43.