列聯分析與對應分析
人們在研究某一個事物或現象的過程中，有些時候不只考察單獨某一方面的信息，即可以把幾個方面的信息聯合起來一並考察。這個過程稱為交叉分析。列聯分析和對應分析就是交叉分析的兩種典型形式，同時也是數據降維分析的一種形式。
列聯分析
對於定類或定序等定性數據的描述和分析，通常可使用列聯表進行分析。本此主要介紹基於列聯表 \(\chi^2\) 檢驗的列聯分析。
列聯表
兩個或兩個以上變量交叉形成的二維頻數分布表格，稱之為列聯表
列聯表中變量的屬性或取值通常也叫做水平，列聯表行變量的水平個數一般用 \(R\) 表示，列變量水平的個數一般用 \(C\) 表示，一個 \(R\) 行 \(C\) 列的頻數分布表叫做 \(R\times C\) 列聯表。
\(R\times C\) 列聯表中各元素 \(f_{ij}\) 就是行列變量進行交叉分類得到的觀測值個數所形成的頻數分布，行合計表示行變量每個水平在列變量不同水平交叉分類的觀測值總數；列合計表示列變量每個水平在行變量不同水平交叉分類的觀測值總數；行合計加總應當等於列合計加總，記為總計頻數。
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 各部門對新工資方案的態度
sc = pd.read_csv('./data/salary_reform.csv')
print(sc.head(5))
# department attitude ID
# 0 3 2 1
# 1 5 2 2
# 2 5 2 3
# 3 4 2 4
# 4 1 1 5
# 數據標簽
sc['department'] = sc['department'].astype('category')
sc['department'].cat.categories = ['發展戰略部', '客戶服務部', '市場部', '研發中心', '綜合部']
sc['department'].cat.set_categories = ['發展戰略部', '客戶服務部', '市場部', '研發中心', '綜合部']
sc['attitude'] = sc['attitude'].astype('category')
sc['attitude'].cat.categories = ['支持', '反對']
sc['attitude'].cat.set_categories = ['支持', '反對']
# 數據列聯表
sc_contingencytable = pd.crosstab(sc['attitude'], sc['department'], margins=True)
print(sc_contingencytable)
# department 發展戰略部 客戶服務部 市場部 研發中心 綜合部 All
# attitude
# 支持 16 21 23 22 22 104
# 反對 25 15 20 27 29 116
# All 41 36 43 49 51 220
# 各個單元格占總人數的百分比
res = sc_contingencytable / sc_contingencytable.loc['All', 'All']
print(res)
# department 發展戰略部 客戶服務部 市場部 研發中心 綜合部 All
# attitude
# 支持 0.072727 0.095455 0.104545 0.100000 0.100000 0.472727
# 反對 0.113636 0.068182 0.090909 0.122727 0.131818 0.527273
# All 0.186364 0.163636 0.195455 0.222727 0.231818 1.000000
# 行/列百分比
def percent_observed(data):
return data / data[-1]
res_r = pd.crosstab(sc['attitude'], sc['department'], margins=True).apply(percent_observed, axis=1)
print(res_r)
# department 發展戰略部 客戶服務部 市場部 研發中心 綜合部 All
# attitude
# 支持 0.153846 0.201923 0.221154 0.211538 0.211538 1.0
# 反對 0.215517 0.129310 0.172414 0.232759 0.250000 1.0
# All 0.186364 0.163636 0.195455 0.222727 0.231818 1.0
res_c = pd.crosstab(sc['attitude'], sc['department'], margins=True).apply(percent_observed, axis=0)
print(res_c)
# department 發展戰略部 客戶服務部 市場部 研發中心 綜合部 All
# attitude
# 支持 0.390244 0.583333 0.534884 0.44898 0.431373 0.472727
# 反對 0.609756 0.416667 0.465116 0.55102 0.568627 0.527273
# All 1.000000 1.000000 1.000000 1.00000 1.000000 1.000000
列聯表的分布
列聯表的分布有兩種：一種是如上表中所示，能夠直接從樣本數據中獲得的交叉分類分布，可以直接觀測得到。其行、列合計分別為行邊緣分布和列邊緣分布；另一種是期望值的分布，是不能直接觀測出來的，可以通過樣本數據和相關理論依據進行計算。
示例
# 以上表為例，如果想要了解不同部門的員工對工資改革方案的態度是否存在顯著差異，在沒有顯著差異的假定條件下，各部分員工不同態度的分布即為列聯表的理論分布。據此可以計算出各部門態度人數的理論期望頻數
# 反對的總人數為116，支持的有104，則對整個單位而言，反對率=116/220=0.5273，支持率=104/220=0.4727.
# 現假定各部門對工資改革的態度沒有差異，故各部門反對該政策的人數=該部門被調查人數*反對率，支持該政策的人數=該部門被調查人數*支持率。由此計算出來的人數便是列聯表的期望值：
# 列聯表的期望分布
from scipy.stats import contingency
res = pd.DataFrame(contingency.expected_freq(sc_contingencytable),
columns=sc_contingencytable.columns,
index=sc_contingencytable.index)
print(res)
# department 發展戰略部 客戶服務部 市場部 研發中心 綜合部 All
# attitude
# 支持 19.381818 17.018182 20.327273 23.163636 24.109091 104.0
# 反對 21.618182 18.981818 22.672727 25.836364 26.890909 116.0
# All 41.000000 36.000000 43.000000 49.000000 51.000000 220.0
\(\chi^2\) 分布與檢驗
列聯表的分布主要有觀測值分布和期望值分布，同時也計算了觀測值和期望值之間的偏差。設 \(f_{ij}^o\) 表示各交叉分類頻數的觀測值，\(f_{ij}^e\) 表示各交叉分類頻數的期望值，則各交叉分類頻數觀測值與期望值的偏差為 \(f_{ij}^o-f_{ij}^e\)，則 \(\chi^2\) 統計量為
當樣本量較大時，\(\chi^2\) 統計量近似服從自由度為 \((R-1)(C-1)\) 的 \(\chi^2\) 分布，\(\chi^2\)值與期望值、觀測值和期望值之差均有關，\(\chi^2\)值越大表明觀測值與期望值的差異越大。
在上例中，假定各部門對工資改革的態度沒有差異，即各部門對該方案的支持率或反對率 \(P\) 均相等，即員工對該方案的態度與其所在部門無關，行列變量之間是獨立的。可以提出原假設和備擇假設
\(H_0\)：部門與對改革方案的態度獨立；\(H_1\)：部門與對改革方案態度不獨立
對假設進行檢驗
# \chi^2檢驗
res = stats.chi2_contingency(sc_contingencytable.iloc[:-1, :-1])
# 參數為不含有margin即邊緣分布的列聯表數據
# 返回\chi^2統計量的值、對應的P值，自由度、列聯表的期望值分布
print(res)
# (4.013295405516321, 0.4042095025927255, 4,
# array([[19.38181818, 17.01818182, 20.32727273, 23.16363636, 24.10909091],
# [21.61818182, 18.98181818, 22.67272727, 25.83636364, 26.89090909]]))
# 統計量值不大，p值比較大，非常不顯著，故無法拒絕原假設
# Fisher精確檢驗:適用與2*2列聯表
res = stats.fisher_exact(sc_contingencytable.iloc[:-1, :-4])
# 返回先驗odds ratio(比值比或相對風險)、P值
print(res)
# (0.45714285714285713, 0.11232756314249981)
# P值較大，非常不顯著，因此，沒有充分理由拒絕發展戰略部、客戶服務部兩個部門態度之間獨立的原假設
\(\chi^2\) 分布的期望准則
\(\chi^2\)檢驗是一種近似檢驗，依據觀測值和期望值斤蒜出來的統計量在大樣本的情況下近似服從 \(\chi^2\) 分布。因此，要求在進行列聯表檢驗的過程中，樣本量應足夠大，且每個交叉分類的期望頻數不能偏小，否則檢驗可能會得出錯誤的結論。
進行 \(\chi^2\) 檢驗時，\(\chi^2\) 分布的期望值准則主要有
1.當交叉分類為兩類時，要求每一類別的期望值不少於5；
2.當交叉分類為兩個以上類別時，期望值小於5的比例不應超過20%，否則應把期望值小於5的類別與相鄰的類別合並
對應分析
\(\chi^2\) 檢驗可以對行列變量之間是否有關聯進行檢驗，但是行列變量之間的關聯性具體是如何作用的，通過列聯表的 \(\chi^2\) 檢驗很難進行判斷。
對應分析便是交叉分析進一步研究的一種方式，它利用數據降維方法直觀明了地分析行列變量之間的相互關系。對應分析是在 \(R\) 型（樣本）和 \(Q\) 型（變量）因子分析的基礎上發展起來的一種多元統計分析方法，它不僅關注行變量或列變量本身的關系，更加關注行列變量之間的相互關系。
對應分析可根據所分析變量的數目分為：簡單對應分析、多重對應分析。簡單對應分析主要用於兩個分類變量之間關系的研究，多重對應分析用於分析3個或更多變量之間的關系。
基本思想
對應分析的基本思想是將一個聯列表的行列變量的比例結構，以散點形式在較低維的空間中表示出來。對應分析省去了因子選擇和因子軸旋轉等中間運算過程，可以從因子載荷圖上對樣品進行直觀的分類，而且能夠指示分類的主因子及分類的依據。此外，對應分析最主要是要得到能夠同時反映眾多樣本和眾多變量的對應分析圖。
為了實現上述基本思想，對應分析通常先找到能夠代表行列變量的行得分與列得分。行得分與列得分互為對方的加權均值，它們之間具有相關性。行、列得分在一個數據中可以得到多組數值，可以根據各組行列得分繪制多個二維散點圖，然后把各個散點圖堆疊起來，最終形成對應分析圖。
為了直觀明了地描述行列變量之間的對應關系，通常選擇2對行列得分，通過2張散點圖疊加得到對應分析圖。在對應分析中，把衡量行列關系強度的指標稱之為慣量（inertia）,其積累所占的百分比成為選取行、列得分對的數目的主要依據。同時，慣量所占比例也成為衡量某對行、列得分在對應分析圖中重要性的重要依據。
步驟和過程
在進行對應分析之前，應當依據列聯分析的知識先判定行列變量之間是否存在相關性。通過檢驗，如果存在相關性，可進行進一步的對應分析，以找出行列變量之間的具體影響關系。對應分析最主要的內容是依據慣量的累積百分確定選取行、列得分的數目，然后繪制對應分析圖，從圖中找出行列變量之間的對應關系。
概率矩陣P
編制兩定性變量的交叉列聯表，得到一個 \(r\times c\) 的矩陣 \(X\)，即
其中，\(r\) 為行變量的分類數，\(c\) 為列變量的分類數，且 \(x_{ij}>0\).
將矩陣 \(X\) 規格化為 \(r\times c\) 的概率矩陣 \(P\)，即
其中，\(p_{ij}=x_{ij}/\sum_{i=1}^{r}\sum_{j=1}^{c}{x_{ij}}\) 為各單元頻數的總百分比。矩陣 \(P\) 表示了一組關於比例的相對數據。
數據點坐標
將 \(P\) 矩陣的 \(r\) 行看成 \(r\) 個樣本，並將這 \(r\) 個樣本看成 \(c\) 維空間中的 \(r\) 個數據點，且各數據點的坐標定義為：\(z_{i1},z_{i2},z_{i3}, \cdots,z_{ic},i=1,2,3,\cdots,r\)。
其中，
此時，各個數據點的坐標是一個相對數據，它在各單元總百分比的基礎上，將在行和列上的分布比例考慮了進來。如果某兩個數據點相距較近，則表明行變量的相應兩個類別在列變量所有類別上的頻數分布差異均不明顯；反之，則差異明顯。
\(P\) 矩陣中，\(p_{ij}\) 表示行變量第 \(i\) 個屬性與列變量第 \(j\) 個屬性同時出現的概率，相應的 \(p_{i.},p_{.j}\) 就有邊緣概率的含義，可定義行剖面：行變量取值固定為 \(i\) 時，列變量各個屬性相對出現的概率情況，即把 \(P\) 中第 \(i\) 行的每一個元素除以 \(p_{i.}\)，同理，可定義列剖面。這樣久可以把第 \(i\) 行表示成在 \(c\) 維歐氏空間中的一個點，其坐標為
其中，\(p_i^r\) 的分量 \(\frac{p_{ij}}{p_{i.}}\) 表示條件概率 \(P(B=j|A=i)\).
第 \(i\) 個行剖面 \(p_i^r\) 就是把 \(P\) 中的第 \(i\) 行剖開單獨研究第 \(i\) 行的各個取值在 \(c\) 維超平面上的分布情況。通過剖面的定義，行變量的不同取值就可用 \(c\) 維空間中的不同點來表示，各點坐標分比為 \(p_j^c\)。
可引入距離概念來分別描述行列變量各個屬性之間的接近程度。行變量第 \(k\) 屬性與第 \(l\) 屬性的歐式距離為：
但這樣定義的距離會受到列變量各屬性邊緣概率的影響，因此可用 \(\frac{1}{p_{.j}}\) 作為權重，得到加權距離公式：
因此，上式定義的距離也可看做是坐標為
的任意兩點之間的歐式距離。
同理，將 \(P\) 矩陣的 \(c\) 列看成 \(c\) 個樣本，並將這 \(c\) 個樣本看成 \(r\) 維空間中的 \(c\) 個數據點，且各數據點的坐標定義為：\(z_{1i},z_{2i},z_{3i}, \cdots,z_{ci},i=1,2,3,\cdots,c\)。
其中，
行列變量分類降維
通過上述步驟能將兩變量的各個類別看作是多維空間上的點，並通過點與點間距離的測度分析類別間的聯系。在變量的類別較多時，數據點所在空間維數必然較高。由於高維空間比較抽象，且高維空間中數據點很難直觀地表示出來，因此最直接的解決方法便是降維。對應分析采用類似因子分析的方式分別對行變量類別和列變量類別實施降維。
如，對列變量實施分類降維：
\(P\) 矩陣的\(c\) 列看做\(c\)個變量，計算 \(c\) 個變量的協方差矩陣 \(A\)。可以證明，第 \(i\) 個變量與第 \(j\) 個變量的協方差矩陣為 \(\sum = (a_{ij})\) ，其中 \(a_{ij}=\sum_{k=1}^{r}{z_{ki}z_{kj}}\)，並記為 \(A=Z^{'}Z\)。從協方差矩陣 \(A\) 出發，計算協方差矩陣 \(A\) 的特征根 \(\lambda_1>\lambda_2>\cdots>\lambda_k,0<k\leq\min\{r,c\}-1\) 以及對應的特征向量 \(\mu_1,\mu_2,\cdots,\mu_k\)，根據累計方差貢獻率確定最終提取特征根的個數 \(m\) （通常 \(m=2\)），並計算出相應的因子載荷矩陣 \(F\)，即
同理，對行變量也可實施類似的分類降維。將 \(P\) 矩陣的 \(r\) 行看做 \(r\) 個變量，計算 \(r\) 個變量的協防擦好矩陣 \(B\)。可以證明，第 \(i\) 個變量與第 \(j\) 個變量的協方差矩陣為 \(\sum = (b_{ij})\) ，其中 \(b_{ij}=\sum_{k=1}^{c}{z_{ik}z_{jk}}\)，並記為 \(B=ZZ^{'}\)。從協方差矩陣 \(B\) 出發，計算協方差矩陣 \(B\) 的特征根和特征向量。可以證明，協方差矩陣 \(A,B\) 具有相同的非零特征根。如果 \(\mu_i\) 為矩陣 \(A\) 的響應特征根 \(\lambda_k\) 的特征向量，那么 \(\nu_i=Z\mu_i\) 就是矩陣 \(B\) 的相應特征根 \(\lambda_k\) 的特征向量，根據累計方差貢獻率確定最終提取特征根的個數 \(m\) ，並計算出相應的因子載荷矩陣 \(G\)，即
對應分析圖
因子載荷矩陣 \(F,G\) 中的元素，其取值范圍時相同的，且元素數量大小的含義也是類似的，因此可以將它們分別看成 \(c\) 個二維點和 \(r\) 個二維點繪制在一個共同的坐標平面中，形成對應分布圖，各點的坐標即為相應的因子載荷。
通過以上步驟，實現了對行列變量多類別的降維，並以因子載荷為坐標，將行列變量的多個分類點直觀地表示在對應分布圖中，實現了定性變量各類別間差異的量化。通過觀察對應分布圖中各數據點的遠近就能判斷個類別之間聯系的強弱。
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 研究不同收入群體購買手機時主要考慮的影響因素
# 指定為黑體中文字體，防止中文亂碼
plt.rcParams["font.sans-serif"] = ["Heiti TC"]
# 解決保存圖像是負號'-'顯示為方塊的問題
plt.rcParams['axes.unicode_minus'] = False
def transform_crosstable(data):
"""
用於將前兩列為行變量各水平，第3列為頻數的原始數據轉換為交叉表
返回行頻數和、列頻數和、交叉表
"""
if data.shape[1] != 3:
raise Exception('Program can only do CA with 2 variables data.')
else:
data = data
# 計算樣本總數n
n = sum(data.iloc[:, 2])
p = data.iloc[:, 2] / float(n)
# 查看每個變量的水平數
v1_len = len(data.iloc[:, 0].value_counts())
v2_len = len(data.iloc[:, 1].value_counts())
v1_name = list(data.iloc[:, 0].unique())
v2_name = list(data.iloc[:, 1].unique())
cross_table = pd.DataFrame(columns=v1_name, index=v2_name)
cross_table_array = np.array(cross_table)
for i in v2_name:
i_index = v2_name.index(i)
cross_table_array[i_index] = list(data[data.iloc[:, 1] == v2_name[i_index]].iloc[:, 2])
cross_table_f = pd.DataFrame(cross_table_array, index=v2_name, columns=v1_name)
total_r = []
total_c = []
for i in range(cross_table_f.shape[0]):
total_r.append(sum(cross_table_f.loc[v2_name[i]]))
for j in range(cross_table_f.shape[1]):
total_c.append(sum(cross_table_f[v1_name[j]]))
total = sum(total_r)
total_r = pd.DataFrame(total_r)
total_c = pd.DataFrame(total_c)
return cross_table_f
if __name__ == '__main__':
cp = pd.read_csv('./data/CellPhone.csv')
print(cp.head(5))
# Element Income Count
# 0 1 1 658
# 1 1 2 665
# 2 1 3 139
# 3 1 4 26
# 4 1 5 13
# Count表示不同Element和Income出現的次數
# 數據標簽
cp['Element'] = cp['Element'].astype('category')
cp['Element'].cat.categories = ['價格', '待機時間', '外觀', '功能', 'IO接口', '網絡兼容性', '內存大小',
'品牌', '攝像頭', '質量口碑', '操縱系統']
cp['Element'].cat.set_categories = ['價格', '待機時間', '外觀', '功能', 'IO接口', '網絡兼容性', '內存大小',
'品牌', '攝像頭', '質量口碑', '操縱系統']
cp['Income'] = cp['Income'].astype('category')
cp['Income'].cat.categories = ['<1000', '1000~3000', '3000~5000', '5000~8000', '8000~10000', '>10000']
cp['Income'].cat.set_categories = ['<1000', '1000~3000', '3000~5000', '5000~8000', '8000~10000', '>10000']
cp = transform_crosstable(cp)
print(cp)
# 價格 待機時間 外觀 功能 IO接口 網絡兼容性 內存大小 品牌 攝像頭 質量口碑 操縱系統
# <1000 658 374 528 332 50 104 143 363 90 626 52
# 1000~3000 665 406 522 323 36 86 165 387 110 637 46
# 3000~5000 139 108 119 76 9 24 46 93 26 147 19
# 5000~8000 26 19 26 13 8 10 5 20 16 27 10
# 8000~10000 13 10 14 11 9 9 8 12 13 15 12
# >10000 13 14 11 6 11 15 7 11 6 18 7
# 創建類Ca的實例，對交叉表進行分析
cp_ca = Ca()
cp_ca.ca(cp)
res = cp_ca.summary()
print(res)
# Singular Value Principal Inertia Chi-square Percent Cumulative Percent
# 0 0.175049 0.030624 243.000000 0.847168 0.847168
# 1 0.058380 0.003407 27.031250 0.094238 0.941406
# 2 0.036804 0.001355 10.750000 0.037476 0.978516
# 3 0.023636 0.000559 4.433594 0.015457 0.994141
# 4 0.014465 0.000209 1.660156 0.005787 1.000000
# 5 0.000000 0.000000 0.000000 0.000000 1.000000
# 第1列為奇異值，即行列變量進行因子分析所得綜合變量的典型相關系數，數值上等於慣量的平方根
# 第2列為慣量，第3列為卡方統計量，其值與列關聯分析中計算的\chi^2值相等，
# 本例中計算出總的\chi^2統計量為286。875,遠大於\alpha=0.05條件下的臨界值，表明行列之間有較強的關聯性
# 第4、5列分別為慣量比例與累積慣量比例。從慣量比例來看，第1維度所占比例=84.72%，其重要性非常大。
# 因此在最后得到的對應分析圖中，主要考察第1維度上的變動情況
# 行列坐標
res_co = cp_ca.get_coords()
print(res_co)
# 行變量坐標： dim1 dim2
# <1000 -0.030768 -0.028231
# 1000~3000 -0.059991 0.016077
# 3000~5000 0.005742 0.034820
# 5000~8000 0.486516 0.147748
# 8000~10000 0.877744 0.198273
# >10000 0.859663 -0.338700
# 列變量坐標： dim1 dim2
# 價格 -0.090981 -0.011639
# 待機時間 -0.031835 -0.004089
# 外觀 -0.058469 0.007325
# 功能 -0.060204 0.012071
# IO接口 0.817753 -0.178094
# 網絡兼容性 0.401819 -0.175205
# 內存大小 0.021973 0.007856
# 品牌 -0.026646 0.015874
# 攝像頭 0.331351 0.199702
# 質量口碑 -0.057744 -0.016852
# 操縱系統 0.671848 0.166509
# 繪制散點圖
cp_ca.rowplot()
cp_ca.colplot()
# 繪制疊加圖
cp_ca.caplot()
# 從圖上可以看到，月收入低於5000的人群在購買手機時主要考慮待機時間、收入、功能、質量等因素，即收入月底想法越多。
# 因此針對中低收入人群的手機應考慮從上述方面來吸引客戶；收入在5000～10000在購買時更加關注攝像頭、操作系統等方面
# 對於收入>10000遠以上的高端用戶，則基本上沒有什么印象因素與之對應。
python沒有較為成熟可靠的對應分析包或模塊，因此，自定義了用於2變量對應分析的類Ca
class Ca(object):
"""
對應分析類
"""
def __init__(self, dim=2):
self.dim = dim # 指定行列變量降維降至的維度，默認2維，無特殊需求使用默認值即可
def ca(self, data):
"""
ca方法用於對應分析的計算過程並存儲該過程中的一些主要結果
"""
# 降數據框數據轉換為數組方便后續操作
data_array = np.array(data)
# 樣本總數
total = data_array.sum()
# 行列百分比
row_percent = data_array.sum(axis=1) / float(total)
row_percent = np.array(row_percent, dtype=float)
col_percent = data_array.sum(axis=0) / float(total)
col_percent = np.array(col_percent, dtype=float)
# 交叉表的概率矩陣形式
p = data_array / float(total)
# 重構
product = np.outer(row_percent.reshape(-1, 1), col_percent.reshape(-1, 1))
center = p - product
# 卡方值
# 總卡方值
chi_squared = float(total) * ((center ** 2) / product).sum()
row_sqrt_I = np.diag(1 / np.sqrt(row_percent))
col_sqrt_I = np.diag(1 / np.sqrt(col_percent))
resid = np.dot(np.dot(row_sqrt_I, center), col_sqrt_I)
resid = np.array(resid, dtype=float)
U, D_lamb, V_T = np.linalg.svd(resid, full_matrices=False)
inertias = D_lamb ** 2
D_lamb_array = np.diag(D_lamb)
# 行列變量前兩個特征值方差貢獻率
# 各主慣量貢獻率
percent = inertias / float(inertias.sum())
# 累計方差貢獻率
cumulative_percent = []
cp = 0
for i in percent:
cp = cp + i
cumulative_percent.append(cp)
# 每個主慣量卡方值
chisquare_list = chi_squared * percent
chisquare_list = [float(item) for item in chisquare_list]
# 完整因子載荷陣
row_coords_all = np.dot(np.dot(row_sqrt_I, U), D_lamb_array)
col_coords_all = np.dot(np.dot(col_sqrt_I, V_T.T), D_lamb_array.T)
# 行變量和列變量坐標(完整因子載荷陣的前2維)
row_coords = np.dot(np.dot(row_sqrt_I, U), D_lamb_array)[:, :self.dim]
row_coords.T[1] = -row_coords.T[1]
col_coords = np.dot(np.dot(col_sqrt_I, V_T.T), D_lamb_array.T)[:, :self.dim]
col_coords.T[1] = -col_coords.T[1]
# 保存類屬性
self.data = data
# 總樣本個數
self.total = total
# 交叉表概率矩陣
self.p = p
# 行列變量因子載荷陣
self.row_coords_all = row_coords_all
self.col_coords_all = col_coords_all
# 二維行列變量坐標
self.row_coords = row_coords
self.col_coords = col_coords
# 奇異值
self.singular_value = D_lamb
# 主慣量
self.principal_inertia = inertias
# 各主慣量方差貢獻率
self.percent = percent
# 各主慣量累積方差貢獻率
self.cumulative_percent = cumulative_percent
# 總卡方值
self.chi_squared = chi_squared
# 各卡方值
self.chisquare_list = chisquare_list
def get_coords(self):
"""
用於輸出行列變量坐標
"""
row_coords_df = pd.DataFrame(data=self.row_coords, index=self.data.index, columns=['dim1', 'dim2'])
col_coords_df = pd.DataFrame(data=self.col_coords, index=self.data.columns, columns=['dim1', 'dim2'])
print('行變量坐標：', row_coords_df)
print('列變量坐標：', col_coords_df)
def summary(self):
"""
用於匯總輸出主要結果
"""
sv = pd.DataFrame(self.singular_value).astype('float16')
pi = pd.DataFrame(self.principal_inertia).astype('float16')
chi = pd.DataFrame(self.chisquare_list).astype('float16')
per = pd.DataFrame(self.percent).astype('float16')
cuper = pd.DataFrame(self.cumulative_percent).astype('float16')
summary_df = pd.concat([sv, pi, chi, per, cuper], axis=1)
summary_df.columns = ['Singular Value', 'Principal Inertia', 'Chi-square', 'Percent', 'Cumulative Percent']
return summary_df
# 考慮到圖像可視化效果，行列變量，雙變量散點圖與類參數dim無關
# 取載荷系數的前兩個維度作為坐標，輸出二維散點圖
def rowplot(self):
"""
行變量圖
"""
x = self.row_coords.T[0]
y = self.row_coords.T[1]
fig1 = plt.figure()
plt.scatter(x, y, c='r', marker='o')
plt.title('Plot of Rows variables')
plt.xlabel('dim1')
plt.ylabel('dim2')
# 為點加標簽
zipxy_r = zip(x, y)
for i, (x, y) in enumerate(list(zipxy_r)):
plt.annotate((list(self.data.index))[i],
xy=(x, y), xytext=(x, y), xycoords="data",
textcoords='offset points', ha='center', va='top')
plt.show()
def colplot(self):
"""
列變量圖
"""
x = self.col_coords.T[0]
y = self.col_coords.T[1]
fig2 = plt.figure()
plt.scatter(x, y, c='b', marker='+', s=50)
plt.title('Plot of Columns variables')
plt.xlabel('dim1')
plt.ylabel('dim2')
# 為點添加標簽
zipxy_c = zip(x, y)
for i, (x, y) in enumerate(zipxy_c):
plt.annotate((list(self.data.columns))[i],
xy=(x, y), xytext=(x, y), xycoords="data",
textcoords="offset points", ha='center', va='top')
plt.show()
def caplot(self):
"""
行列變量圖疊加的對應分析圖
"""
x_c = self.col_coords.T[0]
y_c = self.col_coords.T[1]
x_r = self.row_coords.T[0]
y_r = self.row_coords.T[1]
fig = plt.figure()
plt.scatter(x_r, y_r, c='r', marker='o')
plt.scatter(x_c, y_c, c='b', marker='+', s=50)
plt.title('Plot of two variables')
plt.xlabel('dim1 ({}%)'.format(self.percent[0] * 100))
plt.ylabel('dim2 ({}%)'.format(self.percent[1] * 100))
# 為點添加標簽
zipxy_c = zip(x_c, y_c)
for i, (x, y) in enumerate(zipxy_c):
plt.annotate((list(self.data.columns))[i],
xy=(x, y), xytext=(x, y), xycoords="data",
textcoords="offset points", ha='center', va='top')
zipxy_r = zip(x_r, y_r)
for i, (x, y) in enumerate(zipxy_r):
plt.annotate((list(self.data.index))[i],
xy=(x, y), xytext=(x, y), xycoords="data",
textcoords="offset points", ha='center', va='top')
plt.show()