列联分析与对应分析
人们在研究某一个事物或现象的过程中，有些时候不只考察单独某一方面的信息，即可以把几个方面的信息联合起来一并考察。这个过程称为交叉分析。列联分析和对应分析就是交叉分析的两种典型形式，同时也是数据降维分析的一种形式。
列联分析
对于定类或定序等定性数据的描述和分析，通常可使用列联表进行分析。本此主要介绍基于列联表 \(\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()