Python进行对应分析

从网上找对应分析的资料的时候，发现很多都是理论性的文章，有案例的文章很少。这篇博文主要是利用参考资料中2中的数据，复现一下论文中的整个实验过程。参考资料2中关于对应分析的介绍有错误，所以对应分析的计算过程来源于资料3

1 对应分析

对应分析(Correspondence Analysis)也称关联分析、R-Q型因子分析，是近年新发展起来的一种多元相依变量统计分析技术，通过分析由定性变量构成的交互汇总表来揭示变量间的联系。该技术可以揭示同一变量中各个类别之间的差异，以及不同变量各个类别之间的对应关系。对应分析的基本思想是将一个联列表中的行和列中的各元素的比例结构以点的形式在较低维的空间中表示出来。

2 对应分析过程

2.1 原始矩阵标准化

假设原始矩阵$X$
$$X=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \dots & x_{1p}\\ x_{21}& x_{22}& \dots & x_{2p}\\ ⋮& ⋮& & ⋮\\ x_{n1}& x_{n2}& \dots & x_{np}\end{array}\right]$$其中，行表示样品(比如网站)，列表示属性(比如浏览量)。
(1) 将数据矩阵$X$转化为概率矩阵$P$,其中矩阵$P$中每个元素的定义如下：
$$P=(p_{ij}{)}_{n\times p},p_{ij}=x_{ij}/T,T=\sum _{i=1}^n\sum _{j=1}^p{x}_{ij}$$变换后，概率矩阵$P$的元素之和为1。
(2)求边缘概率$p_{.j}$和$p_{i.}$的边缘概率，其计算公式如下：
$$p_{i.}=\sum _{j=1}^p{p}_{ij},p_{.j}=\sum _{i=1}^n{p}_{ij}$$(3)对概率矩阵$P$做中心化及标准化变换得到矩阵$Z$,其中矩阵$Z$中的元素定义如下：
$$Z=(z_{ij}{)}_{n\times p}=\frac{p_{ij}-p_{i.}{p}_{.j}}{\sqrt{p_{i.}{p}_{.j}}}$$

2.2 计算$R$型因子分析载荷矩阵

(1) 计算协方差矩阵$S_R=Z^{T}Z$的特征值和标准化的特征向量。设特征值${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_m>0$,相应标准化特征向量为$v_1,v_2,\dots ,v_m$。在实际应用中常按累计贡献率
$$\frac{{\lambda }_1+{\lambda }_2+\cdots +{\lambda }_l}{{\lambda }_1+\cdots +{\lambda }_l+\cdots +{\lambda }_m}\ge 0.85$$确定所取公共因子$l\le m$。注意$d_i=\sqrt{{\lambda }_i}，i=(1,2,\dots ,m)$称为奇异值。
(2) $R$型因子的“因子载荷矩阵”为
$$R=\left[\begin{array}{cccc}\frac{\sqrt{{\lambda }_1}{v}_{11}}{\sqrt{p_{.1}}}& \frac{\sqrt{{\lambda }_2}{v}_{12}}{\sqrt{p_{.1}}}& \cdots & \frac{\sqrt{{\lambda }_l}{v}_{1l}}{\sqrt{p_{.1}}}\\ \frac{\sqrt{{\lambda }_1}{v}_{21}}{\sqrt{p_{.2}}}& \frac{\sqrt{{\lambda }_2}{v}_{22}}{\sqrt{p_{.2}}}& \cdots & \frac{\sqrt{{\lambda }_2}{v}_{2l}}{\sqrt{p_{.2}}}\\ ⋮& ⋮& & ⋮\\ \frac{\sqrt{{\lambda }_1}{v}_{p1}}{\sqrt{p_{.p}}}& \frac{\sqrt{{\lambda }_2}{v}_{p2}}{\sqrt{p_{.p}}}& \cdots & \frac{\sqrt{{\lambda }_l}{v}_{pl}}{\sqrt{p_{.p}}}\end{array}\right]$$

2.3 计算$Q$型因子分析载荷矩阵

(1) 计算协方差矩阵$S_Q=Z{Z}^T$的特征值和标准化的特征向量。$S_R$和$S_Q$两个矩阵矩阵大于0的特征值相同的。设$S_Q$的特征值${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_l$其对应的标准化特征向量为$u_1,u_2,\dots ,u_m$。
(2) $Q$型因子的“因子载荷矩阵”为
$$Q=\left[\begin{array}{cccc}\frac{\sqrt{{\lambda }_1}{u}_{11}}{\sqrt{p_{1.}}}& \frac{\sqrt{{\lambda }_2}{u}_{12}}{\sqrt{p_{1.}}}& \cdots & \frac{\sqrt{{\lambda }_l}{u}_{1l}}{\sqrt{p_{1.}}}\\ \frac{\sqrt{{\lambda }_1}{u}_{21}}{\sqrt{p_{2.}}}& \frac{\sqrt{{\lambda }_2}{u}_{22}}{\sqrt{p_{2.}}}& \cdots & \frac{\sqrt{{\lambda }_2}{u}_{2l}}{\sqrt{p_{2.}}}\\ ⋮& ⋮& & ⋮\\ \frac{\sqrt{{\lambda }_1}{u}_{n1}}{\sqrt{p_{n.}}}& \frac{\sqrt{{\lambda }_2}{u}_{n2}}{\sqrt{p_{n.}}}& \cdots & \frac{\sqrt{{\lambda }_l}{u}_{nl}}{\sqrt{p_{n.}}}\end{array}\right]$$

2.4 计算${\chi }^2$统计量和总惯量

(1) 计算总惯量$Q$，总惯量表示$n$个样品到中心$c$的加权平方距离的总和，其最终的计算公式如下：
$$Q=\sum _{i=1}^n\sum _{j=1}^p{z}_{ij}^2=\sum _{i=1}^l{\lambda }_i$$(2) 计算${\chi }^2$统计量，其计算公式如下：
$${\chi }^2=T\sum _{i=1}^n\sum _{j=1}^p{z}_{ij}^2=TQ$$

3 实验过程

本文使用的是参考文献2中的数据。

import pandas as pd
import numpy as np
from math import sqrt
from matplotlib import pyplot as plt
import seaborn as sns
from scipy import stats

path='Documents/data.xlsx'
data=pd.read_excel(path,index_col=0,sheet_name='网站评价分析数据')

#1 卡方检验
# 判断这种分析方法是否适用
data_j=data.sum()/data.sum().sum()
data_i=data.sum(axis=1)
data_exp=np.dot(data_i.values.reshape(-1,1),data_j.values.reshape(1,-1)) 
data_exp=pd.DataFrame(data_exp,index=data_i.index,columns=data_j.index) #期望值
data_obj=data.values.flatten()
data_exp=data_exp.values.flatten()
chi_val,p_val=stats.chisquare(f_obs=data_obj,f_exp=data_exp)
#这里只展示了person卡方分析结果。
#在参考资料4中找到了一个似然卡方公式，但是计算出来的结果和论文中的结果相差太大，所以这里就忽略了

#2 将data转化为概率矩阵,然后进行标准化
data_p=data/data.sum().sum() #概率矩阵
data_j=data_p.sum() #列边缘概率
data_i=data_p.sum(axis=1) #行边缘概率

#概率矩阵标准化
tmp=np.dot(data_i.values.reshape(-1,1),data_j.values.reshape(1,-1))
data_z=(data_p-tmp)/np.sqrt(tmp)
data_z=pd.DataFrame(data_z,index=data.index,columns=data.columns)

#3 求解奇异值和惯量，确定公共因子数量
data_z_np=data_z.to_numpy()
S_R=np.dot(data_z_np.T,data_z_np)
eig_val_R,eig_fea_R=np.linalg.eig(S_R) #返回特征值和特征向量

#注意eig_val_R没有排序,并且只去eig_val_R中大于0的特征值
#返回维度惯量 惯量其实就是特征值
dim_matrix=pd.DataFrame(sorted([i for i in eig_val_R if i>0],reverse=True),columns=['惯量'])
dim_matrix['奇异值']=np.sqrt(dim_matrix['惯量'])
dim_matrix['对应部分']=dim_matrix['惯量']/dim_matrix['惯量'].sum()
dim_matrix['累计']=dim_matrix['对应部分'].cumsum()

#R型因子载荷矩阵
#由dim_matrix['累计']可以得出，我们选择3个公共因子
com_fea_index=[x[0] for x in sorted(enumerate(eig_val_R),reverse=True,key=lambda x:x[1])][:3]
cols=['c'+str(i+1) for i in range(len(com_fea_index))]
eig_fea_R=np.multiply(eig_fea_R[:,com_fea_index],np.sqrt(eig_val_R[com_fea_index]))
R_matrix=pd.DataFrame(eig_fea_R,index=data_j.index,columns=cols)
R_matrix['tmp']=np.sqrt(data_j)
for col in cols:
    R_matrix[col]=R_matrix[col]/R_matrix['tmp']
R_matrix.drop('tmp',axis=1,inplace=True)

#Q型因子载荷矩阵
S_Q=np.dot(data_z_np,data_z_np.T)
eig_val_Q,eig_fea_Q=np.linalg.eig(S_Q)
com_fea_index=[x[0] for x in sorted(enumerate(eig_val_Q),reverse=True,key=lambda x:x[1])][:3]
cols=['c'+str(i+1) for i in range(len(com_fea_index))]
eig_fea_Q=np.multiply(eig_fea_Q[:,com_fea_index],np.sqrt(eig_val_Q[com_fea_index]))
Q_matrix=pd.DataFrame(eig_fea_Q,index=data_i.index,columns=cols)
Q_matrix=Q_matrix.astype(float,copy=False) #Q_matrix中的数据类型为复数类型，但是虚部都为0，所以这里转化成float
Q_matrix['tmp']=np.sqrt(data_i)
for col in cols:
    Q_matrix[col]=Q_matrix[col]/Q_matrix['tmp']
Q_matrix.drop('tmp',axis=1,inplace=True)


#4 选取公共因子c0和c1画定位图
plot_data=pd.concat([Q_matrix[['c1','c2']],R_matrix[['c1','c2']]],axis=0)
plot_data.index=list(data_i.index)+list(data_j.index)
plot_data['style']=['地区']*data_i.shape[0]+['指标']*data_j.shape[0]

#画图
plt.rcParams["font.family"] = 'Arial Unicode MS'
marks={'地区':'o','指标':'s'}
ax=sns.scatterplot(x='c2',y='c1',hue='style',style='style',markers=marks,data=plot_data)
ax.set_xlim(left=-1,right=1)
ax.set_ylim(bottom=-1,top=1)
ax.set_xticks([-1,-0.5,0,0.5,1])
ax.set_yticks([-1,-0.5,0,0.5,1])
ax.axhline(0,color='k',lw=0.5)
ax.axvline(0,color='k',lw=0.5)
for idx in plot_data.index:
    ax.text(plot_data.loc[idx,'c2']+0.005,plot_data.loc[idx,'c1']+0.005,idx)
plt.show()

最后的结果和论文中的结果有些差别，需要对此进行说明：

• 原始数据中有一些缺失值，我在原始数据中用0进行了填充。
• 计算出来的R型因子载荷矩阵和Q型因子载荷矩阵和论文中展示出的结果可能符号相反，误差在百分位上。
  

参考资料

1. 百度百科: 对应分析
2. 基于对应分析法的省级政府门户网站评价研究
3. MATLAB对应分析
4. Pearson and Log-likelihood Chi-square Test of Fit for Latent Class Analysis Estimated with Complex Samples