Correspondence Analysis CASE

案例：
 将1660个人组成的样本按照心理健康状况[good->slight->moderate->damaged]和社会经济状况[Level_A(高)-->LevelE(低)]进行交叉分类如下表所示
 原案例传送门

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

# 频数矩阵
X = np.array([[121,57 ,72 ,36 ,21 ],
              [188,105,141,97 ,71 ],
              [112,65 ,77 ,54 ,54 ],
              [86 ,60 ,94 ,78 ,71 ]
                                   ])

df = pd.DataFrame(X,
                  columns=['Level_A','Level_B','Level_C','Level_D','Level_E'],
                  index=['good','slight','moderate','damaged'])
print(df)

            Level_A  Level_B  Level_C  Level_D  Level_E
good          121       57       72       36       21
slight        188      105      141       97       71
moderate      112       65       77       54       54
damaged        86       60       94       78       71

# 频数总和
X_Sum = np.sum(df.values)
# 频率矩阵
P = df.values/X_Sum
# 行边缘频率
Pr_margin = np.sum(P,1)
# 列边缘频率
Pc_margin = np.sum(P,0)
# 行边缘频率-1/2次方对角方阵
Dr_1_2 = np.diag(np.power(Pr_margin,-0.5))
# 列边缘频率-1/2次方对角方阵
Dc_1_2 = np.diag(np.power(Pc_margin,-0.5))
# 过渡矩阵
Z = np.zeros((P.shape))
for i in range(P.shape[0]):
    for j in range(P.shape[1]):
        Z[i][j] = (P[i][j] - Pr_margin[i]*Pc_margin[j])/np.power(Pr_margin[i]*Pc_margin[j],0.5)
print(Z)

[[ 0.06903397  0.01321385  0.00286337 -0.04560926 -0.07412414]
 [ 0.00748675  0.0022116   0.00362348  0.00224728 -0.02129074]
 [ 0.00335509  0.007487   -0.01807692 -0.01223391  0.02382773]
 [-0.07387793 -0.02171255  0.01038691  0.04952401  0.06934967]]

卡方检验：
 检验行变量和列变量之间(心理健康状况和社会经济状况)是否相互独立，若两者之间相互独立则不用进行对应分析说明两变量之间无关
 .原假设H0：行变量与列变量相互独立

# 检验统计量chi2
chi2 = X.sum() * np.sum(np.power(Z,2))
print('chi2:{}'.format(chi2))
# 显著性水平alpha
alpha = 0.05
# 自由度freedom
freedoom = (Z.shape[0]-1) * (Z.shape[1]-1)
# 显著性水平alpha的临界值
chi2_alpha = stats.chi2.ppf(1-alpha*0.5,freedoom)
print('chi2_alpha:{}'.format(chi2_alpha))

chi2:45.594052238116525
chi2_alpha:23.33666415864534

 因为chi2>chi_alpha，所以拒绝原假设，行和列变量之间不相互独立，可以进行对应分析

# 奇异值分解
U,Sigma,VT = np.linalg.svd(Z)
# 主惯量
inertia = np.power(Sigma,2)
# 主惯量贡献
inertia_contri = np.power(Sigma,2)/np.sum(inertia)
# 打印主惯量贡献表
table1 = pd.DataFrame(np.c_[Sigma,inertia,inertia_contri].T,
                      columns=['dim1','dim2','dim3','dim4'],
                      index=['Sigma','inertia','inertia_contri'])
print(table1)
# 主惯量贡献碎石图
plt.plot([1,2,3,4],inertia_contri,'r*-',linewidth = 2)
plt.xlabel('Strnge Num')
plt.ylabel('Inertia Contribute')

                    dim1      dim2      dim3          dim4
Sigma           0.161318  0.037088  0.008195  9.435155e-17
inertia         0.026024  0.001376  0.000067  8.902214e-33
inertia_contri  0.947475  0.050080  0.002445  3.241141e-31

 从运碎石图可以看出第一维(dim1)贡献了94%的主惯量，第二维(dim2)贡献了5%的主惯量，使用两维作对应分析已经包含了原有数据99%的信息量，且第一维(dim1)占主要地位，对应分析时主要关注dim1。

# 心理健康状况荷载矩阵(取两维)
F = Dr_1_2.dot(U[:,0:2]).dot(np.diag(Sigma[0:2]))
# 心理健康状况标签
mental = ['good','slight','moderate','damaged']
# 社会经济水平荷载矩阵(取两维)
G = Dc_1_2.dot(VT.T[:,0:2]).dot(np.diag(Sigma[0:2]))
# 社会经济水平标签
economic = ['Level_A','Level_B','Level_C','Level_D','Level_E']
# 绘制对应分析图
plt.scatter(F[:,0],F[:,1],c='r',marker='*',linewidths=3)
for i in range(F.shape[0]):
    plt.text(F[i,0],F[i,1],mental[i])
    pass
plt.scatter(G[:,0],G[:,1],c='b',marker='o',linewidths=3)
for j in range(G.shape[0]):
    plt.text(G[j,0],G[j,1],economic[j])
    pass
plt.xlabel('dim1')
plt.ylabel('dim2')
print(F)
print(G)

[[-0.25966303  0.01327174]
 [-0.02945593  0.02257353]
 [ 0.01420067 -0.06995721]
 [ 0.2372966   0.01969363]]
[[-0.1828982  -0.01552113]
 [-0.05902637 -0.02244335]
 [ 0.00888363  0.04233544]
 [ 0.16540302  0.04332697]
 [ 0.28768131 -0.06188021]]

 细心的小伙伴可能已经看出来我们所求得荷载矩阵F等价于原案例中的-X1，我们所求的荷载矩阵G等价于原案例中的-Y1;
 这是因为原案例中的U做了一个变换使U的每一维度(列)的数值之和>1(若U的某一列数值之和小于1则这一维度的数值取负，对应的V的列也取负，以使其满足Z=U.SIGMA.VT)
 U是否做变换不影响对应分析，因为我们只关注相对位置的大小，等比例的放缩不影响相对关系

# 对U和VT进行变换
U_SUM = U.sum(axis=0)
for k in range(U_SUM.shape[0]):
    if U_SUM[k]  < 0:
        U[:,k] = -U[:,k]
        VT.T[:,k] = -VT.T[:,k]
# 心理健康状况荷载矩阵(取两维)
F = Dr_1_2.dot(U[:,0:2]).dot(np.diag(Sigma[0:2]))
# 心理健康状况标签
mental = ['good','slight','moderate','damaged']
# 社会经济水平荷载矩阵(取两维)
G = Dc_1_2.dot(VT.T[:,0:2]).dot(np.diag(Sigma[0:2]))
# 社会经济水平标签
economic = ['Level_A','Level_B','Level_C','Level_D','Level_E']
# 绘制对应分析图
plt.scatter(F[:,0],F[:,1],c='r',marker='*',linewidths=3)
for i in range(F.shape[0]):
    plt.text(F[i,0],F[i,1],mental[i])
    pass
plt.scatter(G[:,0],G[:,1],c='b',marker='o',linewidths=3)
for j in range(G.shape[0]):
    plt.text(G[j,0],G[j,1],economic[j])
    pass
plt.xlabel('dim1')
plt.ylabel('dim2')
print(F)
print(G)

[[ 0.25966303 -0.01327174]
 [ 0.02945593 -0.02257353]
 [-0.01420067  0.06995721]
 [-0.2372966  -0.01969363]]
[[ 0.1828982   0.01552113]
 [ 0.05902637  0.02244335]
 [-0.00888363 -0.04233544]
 [-0.16540302 -0.04332697]
 [-0.28768131  0.06188021]]

简书的markdown语法支持真是差劲！

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

推荐阅读更多精彩内容

  • bash
    官网 中文版本 好的网站 Content-type: text/htmlBASH Section: User ...
    不排版阅读 4,433评论 0赞 5
  • 2017-05-03
    1. 简述相关分析和回归分析的区别和联系。 回归分析和相关分析都是研究两个或两个以上变量之间关系的方法。 广义上说...
    安也也阅读 8,746评论 0赞 3
  • 我只有七块钱了
    我只有七块钱了 它被我牢牢的攥在手里 因为我害怕我把它花完 街边这么多都是我爱吃的好吃的 它们好坏 使劲的飘着香气...
    李子默_阅读 155评论 0赞 0
  • 开局先挂！大巴黎锋线遭遇重大伤病，这下内马尔走不了了？
    大巴黎与巴萨这个夏天因为内马尔转会已经扯皮了两个月时间了，他们目前为止还没有任何达成协议的迹象，然...
    7390b43169fd阅读 150评论 0赞 0
  • 竞争最好的方式是提高自己
    大人谈起竞争，一定有些火气，想着如何压倒对方。而面对所谓的竞争，没有其他方式去应对吗？突然就想到了很多...
    布兰妮田田阅读 434评论 0赞 0