対応分析
Correspondence Analysis
クロス表の対応分析を、Greenacre (2021)に基づいて行うPythonスクリプトを書いてみた。
度数のクロス表が表１のように与えられたとする。
|
表１ 度数のクロス表 |
|
|
|
||
|
|
列１ |
|
列ｊ |
|
列Ｊ |
|
行１ |
|
|
|
|
|
|
|
|
|
|
|
|
|
行i |
|
|
|
|
|
|
|
|
|
|
|
|
|
行Ｉ |
|
|
|
|
|
ここで、
とおく。
比率で表すと、表２のようになる。
|
表２ 比率のクロス表 |
|
|
|
|
||
|
|
列１ |
|
列ｊ |
|
列Ｊ |
行和 |
|
行１ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
行i |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
行Ｉ |
|
|
|
|
|
|
|
列和 |
|
|
|
|
|
1 |
表１の度数のクロス表を
表２の比率のクロス表を
と行列で表せば、
である。
行項目
のプロファイルは
で表されるが、これは
次元空間の点である。行項目の関係は、次元空間の布置として表されている。これを低次元空間に射影して、例えば２次元空間に写せば、平面上に布置を描くことができる。空間における点の布置を低次元空間への射影で近似する方法として特異値分解がある。
まず、
の表す布置の中心が原点になるように次式による変換を行う。
ここで、
である。上式における
は、１のみからなる縦ベクトルであるが、要素の数はそれぞれＪとＩで異なることに注意。
このを次式によりstandardizeする。
ここで、
である。
布置
に対して低次元への射影を特異値分解（SVD）により求める。特異値分解を
と表す。ここで、
である。１次元空間への射影は、｛
、
、｝によって与えられ、２次元空間への射影は｛、
、、
、、｝によって与えられ、３次元空間への射影は｛、、
、、、
、、、｝によって与えられ、・・・というように部分空間への射影は入れ子（nest）になっている。
行項目の点を射影した座標は
で与えられるが、標準座標（standard coordinate）は
で与えられる。これに対して、特異値でスケールした
は主座標（principal coordinate）と呼ばれている。
次式が成り立つ。
ここで、
とおいた。
列項目も同様に考える。
列項目
のプロファイルは
で表されるが、これは
次元空間の点である。列項目の関係は、次元空間の布置として表されている。これを低次元空間に射影して、例えば２次元空間に写せば、平面上に布置を描くことができる。空間における点の布置を低次元空間への射影で近似する方法として特異値分解がある。
まず、の表す布置の中心が原点になるように次式による変換を行う。
このを次式によりstandardizeする。
布置に対して低次元への射影を特異値分解（SVD）により求める。特異値分解を
と表す。
列項目の点を射影した座標は
で与えられるが、標準座標（standard coordinate）は
で与えられる。これに対して、特異値でスケールした
は主座標（principal coordinate）と呼ばれている。
次式が成り立つ。
次式の値はtotal inertiaと呼ばれている。
いま、
とおく。
であるので、
より、
すなわち、
式（１）の値は、the k-th principal inertiaと呼ばれている。
Inertiaについて行項目の観点から説明したが、列項目の観点からも同様に論じられる。
上の説明に基づいて対応分析のPythonスクリプトをリスト１のように作成した。スクリプトファイルとサンプルデータファイルは、まとめて圧縮ファイルcorresanalgrnacrfiles.zipとした。ダウンロード展開すれば自由に利用できる。
リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。
= RESTART:
D:\・・・\corresanalgrnacrfiles\CorresAnal.py
Data File
Name (*.xlsx) =
入力データファイルはExcelファイルの形式で用意する。図１は、政府統計ポータルサイトe-Statからダウンロードした「学校における教育の情報化の実態等に関する調査（令和３年度）」のファイルをリスト１のスクリプト用に編集したものである。
図１ 「学校における教育の情報化の実態等に関する調査（令和３年度）」（e-Statより）を編集
１行目に列項目名、１列目に行項目名が設定されている。セル内の数値はPC（デバイス）の設置台数である。
図１のファイル名情報調査結果.xlsxを設定する。
Data File
Name (*.xlsx) = 情報調査結果.xlsx
ファイル名を設定して、Enterキーを押すと計算が始まる。
計算が進むと、まずscree plot（図２）が表示される。
図２
第１軸でデータのかなりの部分が説明され、２軸もデータの傾向を若干表していることが読み取れる。
図２のフォームを閉じると、図３のフォームが提示される。
図３
行項目をprincipal coordinaesで、列項目をstandard coordinatesで表示したasymmetric mapである。
図３のフォームを閉じると、端末画面に次のメッセージが表示される。
Row items/principal coordinates:
Column items/standard coordinates
x-axis =
次に描く図３のようなマップで用いる横軸の次元の設定が求められている。横軸と縦軸の組み合わせを変えて表示したいときはそれぞれに次元を順番に設定していけばよいが、０以下の値（整数値）を設定すると、図３の形式のマップの表示は終わり、次の形式のマップが表示される（図４）。
図４
列項目がprincipal coordinaesで、行項目をstandard coordinatesで表示されたasymmetric mapである。
図４のフォームを閉じると、端末画面に次のメッセージが表示される。
Row items/standard coordinates: Column
items/principal coordinates
x-axis =
横軸と縦軸の次元を変えて図４の形式のマップを表示するときは、それぞれの次元を設定すればよい。次の形式のマップ（図５）に進むときは、０以下の整数値を設定する。
図５
図５は、列項目、行項目ともにprincipal coordinatesのsymmetric mapである。図１のデータに認められる「特別支援学校ではiOSまたはiPadの度数が多いこと」、「小学校・中学校ではChrome OSが多いこと」、「高等学校・高等学校（学科設置校）ではWindowsが主流であること」に対応した布置が見られる。
図５のフォームを閉じると、次のメッセージが表示される。
Row and Column items/principal
coordinates
x-axis =
次元の組み合わせを変えて図５の形式のマップを表示するときは、次元を設定する。図５の形式のマップの表示を終了するときは、０以下の整数を設定する。
０以下の整数値を設定すると、プログラムの実行終了となり、端末には以下のメッセージが表示される。
Row and Column items/principal
coordinates
x-axis = 0
Results.txt was saved.
出力ファイルResults.txtの内容は、以下の通りである。
Input data file = 情報調査結果.xlsx
Data =
Unnamed: 0 Windows 11 Windows 10 ...
Chrome OS Android その他のOS（Linux等）
0
小学校 105148 2306588 ... 2560813 4656
1780
1
中学校 56976 1271639 ... 1404207 3385
529
2 義務教育学校
527
28248 ... 18223 432
1
3
高等学校 58870 702666 ... 424957 10796
2629
4
高等学校（学科設置校） 40230 447571 ... 228314 5386
1934
5 中等教育学校
369
12311 ... 6250 74
20
6 特別支援学校
2705
48263 ... 4158 496
604
[7 rows x 10 columns]
Principal inertias
Axis Principal inertia % inertia Cumulative %
1
0.05869
77.7
77.7
2
0.01603
21.2
98.9
3
0.00041
0.5
99.5
4
0.00028
0.4
99.8
5 0.00012
0.2
100.0
6
0.00001
0.0
100.0
Total
0.07555
「Data =」に続く、読み込まれたデータの内容の表示のフォーマットに乱れが見られるのは、Pythonの日本語を含む文字列の書式制御がうまくいかないためである。内容に問題は無い。
Greenacre, M. (1984). Theory and Applications of Correspondence Analysis. Academic Press
Greenacre, M. (2021). Correspondence Analysis in Practice, third edition (paperback). CRC Press.
岡本安晴（2019）いまさら聞けないPythonでデータ分析：多変量解析、ベイズ統計分析（PyStan, PyMC）．丸善出版
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
data = pd.read_excel(xlsx)
fout = open('Results.txt', 'w')
fout.write(f'Input data file =
{flnm}\n')
print(data)
fout.write(f'\nData =\n{data}\n')
N = data.values
row_names = N[:,0]
col_names = data.keys()[1:]
N = N[:,1:]
sum_row = np.sum(N, axis=1)
sum_col = np.sum(N, axis=0)
n = np.sum(N)
P = N/n
I, J = np.shape(P)
print('\nP =')
for i in range(I):
for j in range(J):
print(f'{P[i][j]:>10.4f}', end = '')
print()
print()
r = P @ np.ones((J, 1))
c = P.T @ np.ones((I, 1))
print('r =')
for v in r:
print(f'{v[0]:>10.4f}')
print('\nc =')
for v in c:
print(f'{v[0]:>10.4f}')
Dr = np.diag(r[:,0])
Dc = np.diag(c[:,0])
Dr_n5 = np.diag(r[:,0]**-0.5)
Dc_n5 = np.diag(c[:,0]**-0.5)
S = Dr_n5 @ (P - r @ c.T) @ Dc_n5
S = np.array(S, dtype='float')
U, Da, Vh = np.linalg.svd(S,
full_matrices = False)
ck = [1 if v > 1.0e-6 else 0 for v
in Da]
K = sum(ck)
print('\nK =', K)
U = U[:,:K]
Da = np.diag(Da[:K])
V = Vh[:K].T
Phi = Dr_n5 @ U
Gmm = Dc_n5 @ V
print('\nPhi =')
for i, v in zip(row_names, Phi):
print(f'{i:>20s}',
end = '')
for v1 in v:
print(f'{v1:>10.4f}', end = '')
print()
print('\nGmm =\n')
for i, v in zip(col_names, Gmm):
print(f'{i:>20s}', end =
'')
for v1 in v:
print(f'{v1:>10.4f}', end = '')
print()
F = Phi @ Da
G = Gmm @ Da
print('\nF =')
for i, v in zip(row_names, F):
print(f'{i:>20s}', end =
'')
for v1 in v:
print(f'{v1:>10.4f}', end = '')
print()
print('\nG =')
for i, v in zip(col_names, G):
print(f'{i:>20s}', end =
'')
for v1 in v:
print(f'{v1:>10.4f}', end = '')
print()
Dlmbd = Da @ Da
alphas = np.diag(Da)
lambdas = np.diag(Dlmbd)
cum_lambdas = lambdas.cumsum()
pcum_lambdas = cum_lambdas /
cum_lambdas[-1]
p_lambdas = lambdas / cum_lambdas[-1]
print('\nPrincipal inertias')
print(f'{"Axis
":>10s}{"Principal inertia":>20s}{"%
inertia":>12s}{"Cumulative %":>15s}')
for i in range(K):
print(f'{i+1:>10d}', end
= '')
print(f'{lambdas[i]:>20.5f}', end = '')
print(f'{p_lambdas[i]*100:>12.1f}', end = '')
print(f'{pcum_lambdas[i]*100:>15.1f}')
print(f'{"Total":>10s}{cum_lambdas[-1]:>20.5f}'
)
print()
fout.write('\nPrincipal inertias\n')
fout.write(f'{"Axis
":>10s}{"Principal inertia":>20s}' +\
f'{"% inertia":>12s}{"Cumulative %":>15s}\n')
for i in range(K):
fout.write(f'{i+1:>10d}')
fout.write(f'{lambdas[i]:>20.5f}')
fout.write(f'{p_lambdas[i]*100:>12.1f}')
fout.write(f'{pcum_lambdas[i]*100:>15.1f}\n')
fout.write(f'{"Total":>10s}{cum_lambdas[-1]:>20.5f}\n'
)
plt.bar(np.arange(K),lambdas)
plt.xticks(np.arange(K),
np.arange(K)+1)
plt.xlabel('Axis', fontsize = 18)
plt.ylabel('Principal inertia',
fontsize = 18)
plt.tight_layout()
plt.show()
#
#
Row items/principal coordinates
#
Column items/standard coordinates
#
dim1 = 1
dim2 = 2 if K > 1 else 1
while True:
plt.figure(figsize = (7,7))
plt.title('Row items:
Principal coordinates\n' +
'Column items: Standard coordinates')
for i in range(I):
plt.plot(F[i][dim1-1],
F[i][dim2-1], 'o', c = 'g')
plt.text(F[i][dim1-1], F[i][dim2-1], row_names[i])
for j in range(J):
plt.plot(Gmm[j][dim1-1], Gmm[j][dim2-1], '^', c = 'orange')
plt.text(Gmm[j][dim1-1], Gmm[j][dim2-1], col_names[j])
v1 = np.max([F[:,dim1-1],
F[:,dim2-1]])
v2 = np.max([-F[:,dim1-1],
-F[:,dim2-1]])
v3 = np.max([Gmm[:,dim1-1],
Gmm[:,dim2-1]])
v4 = np.max([-Gmm[:,dim1-1],
-Gmm[:,dim2-1]])
max_v = np.max([v1, v2, v3,
v4])
plt.plot([-max_v, max_v],
[0, 0], ls = ':', c = 'y')
plt.plot([0, 0], [-max_v,
max_v], ls = ':', c = 'y')
plt.xlabel(f'Axis {dim1}:
{lambdas[dim1-1]:.4f} ({100*p_lambdas[dim1-1]:.1f}%)')
plt.ylabel(f'Axis {dim2}:
{lambdas[dim2-1]:.4f} ({100*p_lambdas[dim2-1]:.1f}%)')
plt.show()
print('Row items/principal
coordinates: Column items/standard coordinates')
dim1 = int(input('x-axis =
'))
if dim1 < 1:
break
if dim1 > K:
dim1
= K
dim2 = int(input('y-axis =
'))
if dim2 < 1:
dim2 = 1
if dim2 > K:
dim2
= K
#
# Row
items/standard coordinates
#
Column items/principal coordinates
#
dim1 = 1
dim2 = 2 if K > 1 else 1
while True:
plt.figure(figsize = (7,7))
plt.title('Row items:
Standard coordinates\n' +
'Column items: Principal coordinates')
for i in range(I):
plt.plot(Phi[i][dim1-1], Phi[i][dim2-1], 'o', c = 'g')
plt.text(Phi[i][dim1-1], Phi[i][dim2-1], row_names[i])
for j in range(J):
plt.plot(G[j][dim1-1], G[j][dim2-1], '^', c = 'orange')
plt.text(G[j][dim1-1], G[j][dim2-1], col_names[j])
v1 = np.max([Phi[:,dim1-1],
Phi[:,dim2-1]])
v2 = np.max([-Phi[:,dim1-1],
-Phi[:,dim2-1]])
v3 = np.max([G[:,dim1-1],
G[:,dim2-1]])
v4 = np.max([-G[:,dim1-1],
-G[:,dim2-1]])
#print('v4 =', v4)
#print(-Gmm[:,0:1])
max_v = np.max([v1, v2, v3,
v4])
plt.plot([-max_v, max_v],
[0, 0], ls = ':', c = 'y')
plt.plot([0, 0], [-max_v,
max_v], ls = ':', c = 'y')
plt.xlabel(f'Axis {dim1}:
{lambdas[dim1-1]:.4f} ({100*p_lambdas[dim1-1]:.1f}%)')
plt.ylabel(f'Axis {dim2}:
{lambdas[dim2-1]:.4f} ({100*p_lambdas[dim2-1]:.1f}%)')
plt.show()
print('Row items/standard
coordinates: Column items/principal coordinates')
dim1 = int(input('x-axis =
'))
if dim1 < 1:
break
if dim1 > K:
dim1
= K
dim2 = int(input('y-axis =
'))
if dim2 < 1:
dim2
= 1
if dim2 > K:
dim2
= K
#
#
Symmetric Map
#
dim1 = 1
dim2 = 2 if K > 1 else 1
while True:
plt.figure(figsize = (7,7))
plt.title('Row items:
Principal coordinates\n' +
'Column items: Principal coordinates')
for i in range(I):
plt.plot(F[i][dim1-1], F[i][dim2-1], 'o', c = 'g')
plt.text(F[i][dim1-1],
F[i][dim2-1], row_names[i])
for j in range(J):
plt.plot(G[j][dim1-1], G[j][dim2-1], '^', c = 'orange')
plt.text(G[j][dim1-1], G[j][dim2-1], col_names[j])
v1 = np.max([F[:,dim1-1],
F[:,dim2-1]])
v2 = np.max([-F[:,dim1-1],
-F[:,dim2-1]])
v3 = np.max([G[:,dim1-1],
G[:,dim2-1]])
v4 = np.max([-G[:,dim1-1],
-G[:,dim2-1]])
max_v = np.max([v1, v2, v3,
v4])
plt.plot([-max_v, max_v],
[0, 0], ls = ':', c = 'y')
plt.plot([0, 0], [-max_v, max_v],
ls = ':', c = 'y')
plt.xlabel(f'Axis {dim1}:
{lambdas[dim1-1]:.4f} ({100*p_lambdas[dim1-1]:.1f}%)')
plt.ylabel(f'Axis {dim2}:
{lambdas[dim2-1]:.4f} ({100*p_lambdas[dim2-1]:.1f}%)')
plt.show()
print('Row and Column
items/principal coordinates')
dim1 = int(input('x-axis =
'))
if dim1 < 1:
break
if dim1 > K:
dim1
= K
dim2 = int(input('y-axis =
'))
if dim2 < 1:
dim2
= 1
if dim2 > K:
dim2
= K
fout.close()
print('Results.txt was saved.')