クロス表の視覚化とクラスター分析
特異値分解(SVD)と内積モデル
クロス表(対応表)の行情報と列情報を視覚化する。視覚化することにより、クロス表に含まれる情報が直感的に読み取れることが期待される。まず、モデルについて説明を行い、次にクロス表の分析例としてCATAデータの分析例を示す。分析例では、特異値分解で得られた行情報と列情報のBiPlotと内積モデルに基づくクラスター分析(座標値によるものとベクトル間のなす角によるもの)を示した。
一般に、配列の行情報と列情報を取り出す方法として特異値分解(singular
value decomposition: SVD)がある(Brunton
& Kutz, 2019)。特異値分解SVDでは、行列を次式で展開する(岡本、2019、p.
67;岡本、2009、p.47)。
(1)
クロス表の場合、各ベクトルとをの平方根で重み付けたもの、と、を座標値として重ねて図示すると行および列に含まれる関係が視覚的に読み取れる。行の座標値と列の座標値を重ねて図示すると行と列の関係が1つの図に表され、関係を読み取るのに便利である。
岡太彬訓
立教大学名誉教授より、式(1)を行列の成分の内積による表現として見るという観点のコメントを頂いた(私信、2022.11.02)。この観点では、式(1)を以下のように展開する。
特異値分解を次式で表す。
ここで、
, ,
である。
いま、
,
と置くと、
である。すなわち、
を得る。
クロス表の各セルの値が内積で与えられる。これらの内積を与えるベクトルとは、同じ空間内のベクトルであり、所与のデータを、内積で表現するモデルである。
式(2)を変形して、
,
を得る。内積モデルにおいては、座標の直交回転ができる。
また、内積モデルで考えると、クラスター分析が可能である。
内積を与える座標ベクトルとは、同じ空間内のベクトルであるので、これに基づいてクラスター分析を行うことができる。また、2つのベクトルのなす角に基づいたクラスター分析を行うこともできる。
クロス表の例としてCATAデータを用いてみる(岡本、2022)。島村ら(2017)のTable 2のデータをExcelに設定したものを図1に示す。
図1
図1のデータ(ExcelファイルCATAData.xlsx)をSVDで分析するスクリプトCATADataSVD.pyをリスト1のように作成した。スクリプトファイルCATADataSVD.pyとサンプルデータファイルCATAData.xlsxは圧縮ファイルsamplefiles.zipにまとめた。ダウンロード・展開して自由に利用されたい。
リスト1のスクリプトを実行すると、まず入力データファイル名の設定が求められる。
Excel Data File(*.xlsx) =
入力データファイルは、図1の形式のExcelファイルの名前を設定する。図1のデータファイルであれば、CATAData.xlsxである。
Excel Data File(*.xlsx) = CATAData.xlsx
入力データファイルは、スクリプトファイルと同じフォルダに置いておくとよい。
入力データファイル名を設定して、Enterキーを押すと計算が始まる。
計算後、まず、特異値分解による行情報と列情報をBiPlotとして重ねたグラフ(図2)が描かれる。
図2
第1軸(横軸)の値で見ると、米は、コシヒカリ、ミルキークイーン、ササニシキのグループと低たんぱく米とに分かれている。評価語は、「好きな味」と「程良い粘り」がまとまっており、コシヒカリ、ミルキークイーン、ササニシキのグループと近い位置である。これは、図1における「好きな味」と「程良い粘り」の頻度がコシヒカリ、ミルキークイーン、ササニシキにおいて高いことに対応している。第2軸(縦軸)の座標について見てみると、値が0から離れているものとしては、米ではササニシキが正の側、低たんぱく米が負の側にある。評価語を見ると、正の側に「好きな味」があり、ササニシキに近い位置である。負の側には、「嫌いな味」、「ねっとりした」、「渋味」があり、低たんぱく米の近くに位置している。図3におけるこれらの組み合わせに対する分割表の頻度の値が高いことに対応している。
なお、図2において横軸の値が全て負であるが、特異値分解の場合、ベクトルとが特異ベクトルの組であれば、それらの符号を同時に反転したベクトルとも特異ベクトルの組であるので、第1軸において、米と評価語の値を同時に反転すれば、第1軸の値が正である布置を得る。
内積を与える座標ベクトルとは、同じ空間内のベクトルであるので、これに基づいてクラスター分析を行うことができる。また、2つのベクトルの関係をそのなす角で表し、角度に基づいてクラスター分析を行うこともできる。リスト1のスクリプトでは、クラスター分析も行っている。図2のグラフが表示されているウィンドウを閉じると、座標ベクトルとに基づいてWard法によるクラスター分析を行った結果が図3のように表示される。
図3
図3のウィンドウを閉じると、ベクトル間の角によるクラスター分析をMaximum(Complete, Farthest Point or Voor
Hees)アルゴリズムにより行った結果が示される(図4)。
図4
図1のデータの場合、図3のクラスター分析では、ミルキークイーン、コシヒカリ、ササニシキがまとまり、低たんぱく米が離れたクラスターに属している。「程良い粘り、好きな味」と「嫌いな味、渋味、ねっとりした、酸味」に分かれていると解釈できる。
これに対して、図4のクラスター分析では、4種の米は、それぞれのクラスターを形成していて、コシヒカリは甘い味、ミルキークイーンは程良い粘り、ササニシキはかたい、低たんぱく米は嫌いな味と特徴付けられていると解釈できる。
得られた内積モデルによる分析結果に対して、クラスター分析で注目する特徴(ベクトルの座標、ベクトル間の角度)によりデータの異なった側面が表される。
図4のグラフが表示されているウィンドウを閉じると、スクリプトの実行終了となる。
図2、図3、図4のグラフは、それぞれファイルFigMap.png、FigCluster.png、FigClusterAngles.pngに保存されている。
Brunton, S. L., & Kutz, J. N. (2019). Data-driven science and engineering. New York: Cambridge University Press.
岡太彬訓(私信、2022.11.02)
岡本安晴(2009) 統計学を学ぶための数学入門[下].培風館
岡本安晴(2019) いまさら聞けないPythonでデータ分析.丸善出版
岡本安晴(2022) 官能評価に使える統計解析ツール Python編、日本官能評価学会誌、第26巻第2号、102-106.
島村 綾、小泉昌子、峯木眞知子、市原 茂(2017) 飯の官能評価の時系列変化.日本家政学会誌、68、9、478-485.
リスト1 CATAデータのSVD(CATADataSVD.py)
import pandas as pd
import numpy as np
import
matplotlib.pyplot as plt
import
japanize_matplotlib
from
scipy.cluster.hierarchy import linkage
from
scipy.cluster.hierarchy import dendrogram
# ExcelファイルをPandasのDataFrame型として読み込む
fl_nm = input('Excel
Data File(*.xlsx) = ')
xlfile =
pd.ExcelFile(fl_nm)
df =
pd.read_excel(xlfile)
# 1行目の値(米の種類)
col_names =
df.keys()
# Excelの1行目
col_names =
col_names[1:]
# 2番目以降を取り出す
# '評価用語'の列を取り出す
row_names = df['評価用語']
M = len(row_names)
N = len(col_names)
# セル内の数値
df_data =
df.values
# 2行目以降が取り出される
df_data =
df_data[:,1:] # 取り出した各行の2列目以降を取り出す
# セル内の数値を実数型に変換
data =
np.asfarray(df_data)
# dataの特異値分解SVD
u, s, vh =
np.linalg.svd(data, full_matrices = False)
print('\ns =', s)
# 特異値の平方根による重み付け
col_coord =
np.diag(s**0.5) @ vh
row_coord = u @
np.diag(s**0.5)
#
# 重み付けられた特異ベクトルによる描画
#
plt.plot(row_coord.T[0],row_coord.T[1],'^',c
= 'k')
plt.plot([], '^',c =
'k',label = '評価項目')
plt.plot(col_coord[0],col_coord[1],'v',c
= 'k')
plt.plot([], 'v', c =
'k', label = '米')
for i in range(len(col_names)):
plt.text(col_coord[0][i],
col_coord[1][i],
col_names[i])
for i in
range(len(row_names)):
plt.text(row_coord.T[0][i],
row_coord.T[1][i],
row_names[i],rotation = -20,
rotation_mode = 'anchor')
plt.legend()
plt.savefig('FigMap.png')
plt.show()
#
# クラスター分析のためにクロス表の行と列の座標を統合
#
mrg_coord =
np.concatenate((row_coord, col_coord.T), axis = 0)
mrg_label =
np.concatenate((row_names, col_names), axis = 0)
#
# 座標値によるクラスター分析
#
row_clusters =
linkage(mrg_coord, method = 'ward')
fig =
plt.figure(figsize = (10, 7), facecolor = 'white')
row_dendr =
dendrogram(row_clusters, labels = mrg_label,
orientation = 'right')
plt.title('座標値によるクラスター分析',
fontsize = 14)
plt.tight_layout()
plt.savefig('FigCluster.png')
plt.show()
#
#
ベクトルxとyのなす角
#
def calc_angle(x, y):
#print(x,y)
return np.arccos(np.dot(x,
y) / ((np.dot(x,x))*np.dot(y,y))**0.5)
Angles =
np.empty((M+N, M+N))
#print(Angles)
for i in range(M+N):
for j in range(M+N):
if i
> j:
Angles[i][j] = Angles[j][i]
elif
i == j:
Angles[i][j] = 0.0
else:
Angles[i][j] = calc_angle(mrg_coord[i], mrg_coord[j])
#
# 角度によるクラスター分析
#
dists = []
for i in
range(len(Angles)-1):
for j in range(i+1,
len(Angles)):
dists.append(Angles[i][j])
row_clusters =
linkage(dists, method = 'complete')
fig =
plt.figure(figsize = (7, 5), facecolor = 'white')
row_dendr =
dendrogram(row_clusters, labels = mrg_label,
orientation = 'right', leaf_font_size = 14) #, link_color_func = lambda
k: 'k')
plt.title('角度によるクラスター分析',
fontsize = 14)
plt.tight_layout()
plt.savefig('FigClusterAngles.png')
plt.show()