多元统计分析-基于R（第3版） 课件 费宇 第12、13章-典型相关分析、多维标度分析

2024/9/5主编：费宇2第12章典型相关分析相关系数可以衡量两个变量间的相关关系，但两组变量之间的相关关系如何来度量呢？本章讨论的典型相关分析(canonicalcorrelationanalysis)就是研究两组变量之间相关关系的一种多元统计分析方法，它利用主成分的思想来讨论两组随机变量的相关性问题，分别对两组变量提取主成分，通过它们的相关性来度量两组变量整体的线性相关关系.典型相关分析的思想首先由Hotelling于1936年提出，现在其已经成为一种常用的分析两组变量相关性的多元分析方法，在实际中应用广泛.2024/9/5主编：费宇3第12章典型相关分析12.1

典型相关分析基本理论12.2

总体典型相关变量的概念及其解法12.3典型相关变量的性质12.4原始变量与典型相关变量的相关系数12.5简单相关、复相关和典型相关之间的关系12.6分量的标准化处理2024/9/5主编：费宇4第12章典型相关分析12.7样本典型相关系数及其对应典型相关变量的计算12.8典型相关系数的显著性检验12.9被解释样本方差的比例12.10案例分析及R操作2024/9/5主编：费宇512.1

典型相关分析基本理论典型相关分析是研究两组变量之间整体的相关关系,它将每一组变量作为一个整体来进行研究,所研究的两组变量可以是一组变量是自变量,另一组变量是因变量;当然,也可以两组变量处于同等地位.2024/9/5主编：费宇612.1

典型相关分析基本理论典型相关分析的基本原理借助主成分分析的思想,在每组变量中找出变量的线性组合即新的综合变量,使生成的综合变量能代表原始变量的主要信息,同时,与由另一组变量生成的新的综合变量的相关程度最大,这样得到的一组新变量称为第一对典型相关变量;同样的方法可以找到第二对典型相关变量,第三对典型相关变量…，要求各对典型相关变量之间互不相关.典型相关变量间的相关系数称为典型相关系数,它度量了这两组变量之间关系的强度.此项最大化技术是努力将两组变量间的一个高维关系浓缩到用少数几个典型变量来表现．2024/9/5主编：费宇712.2总体典型相关变量的概念及其解法1.总体典型相关变量假设有两组变量,一组变量为x=(x1,x2,…,xp

)T

,另一组变量为y=(y1,y2,…,yq

)T,且p≤q

,变量x与变量y的协方差阵为:2024/9/5主编：费宇812.2总体典型相关变量的概念及其解法为研究变量x与变量y之间的线性相关关系,我们考虑它们之间的线性组合u和v的方差和协方差分别为2024/9/5主编：费宇912.2总体典型相关变量的概念及其解法两个新变量u和v之间的相关系数（即典型相关系数）为由于变量u和v乘以不为零常数不改变它们之间的相关性,即对任意常数c≠0,d≠0,有Corr(cu,cv)=Corr(u,v),所以通常需对a和b附加约束条件,使变量u和v不必要的重复,最好的约束条件是2024/9/5主编：费宇1012.2总体典型相关变量的概念及其解法我们的问题就变成在上述约束条件下求a和b,使得达到最大，于是有以下定义．定义12.1设，，维随机向量的均值向量为0，协方差矩阵

（不妨设）.如果存在2024/9/5主编：费宇1112.2总体典型相关变量的概念及其解法和，令，，使得这样得出的和称为x，y的第一对（组）典型相关变量，称为第一个典型相关系数；如果存在和使得1），和前面i-1对典型相关变量不相关；2）,;3）与的相关系数最大；2024/9/5主编：费宇1212.2总体典型相关变量的概念及其解法则称，是x，y的第i对（组）典型相关变量，它们之间的相关系数称为第i个典型相关系数（）．由拉格朗日乘数法,这个问题等价于求a和b使达到最大,其中μ1和μ2是拉格朗日乘数.2024/9/5主编：费宇1312.2总体典型相关变量的概念及其解法将(12.7)两边分别对向量a和b求导,并令其为0,得方程组以aT和bT分别左乘(12.8)两式得但,所以μ1=μ2=ρ,即μ1恰好就是u和v的相关系数.两边左乘以

得,同理可得,记,则得2024/9/5主编：费宇1412.2总体典型相关变量的概念及其解法由方程组(12.8)的第二式得,将其代入(12.8)的第一式得

M1和M2的非零特征值皆为正数，，既是M1的特征根又是M2的特征根,a和b分别是M1和M2相相应的特征向量.于是求

和a,b的问题就转化为求矩阵M1和M2的特征根和特征向量的问题.设ai是M1的属于的特征向量，令2024/9/5主编：费宇1512.2总体典型相关变量的概念及其解法有则是M2的属于的特征向量.2024/9/5主编：费宇1612.2总体典型相关变量的概念及其解法设M1的m个正特征根为,相对应的特征向量分别为a1,a2,…,am由式(12.10)得出且正交化，b1,b2,…,bm由式(12.11)得出且正交化,

从而可得m对线性组合每一对变量称为一对典型变量,其中u1和v1称为第一对典型变量,它们之间的相关系数ρ1即为第一典型相关系数.

ui和vi称为第i对典型变量,它们之间的相关系数ρi即为第i典型相关系数.2024/9/5主编：费宇1712.3典型相关变量的性质我们给出典型变量以下四个性质（证明见本章附录12）:(1)每一对典型变量ui及vi(i=1,2,…,m)的标准差为1.(2)同一组的任意两个典型变量ui(i=1,2,…,m)彼此不相关,典型相关变量vi(i=1,2,…,m)彼此不相关,,，

.(3)不同组的任意两个典型变量ui,vj(i=1,2,…,m;2024/9/5主编：费宇1812.3典型相关变量的性质j=1,2,…,m)的关系为：(4)典型变量ui及vi的相关系数为ρi(i=1,2,…,m)，典型相关系数满足关系式．2024/9/5主编：费宇1912.3典型相关变量的性质在理论上,典型变量的对数和相对应的典型相关系数的个数可以等于两组变量中数目较少的那一组变量的个数,其中,u1及v1的相关系数ρ1反映的相关成分最多,所以称为u1,v1第一对典型变量;u2及v2的相关系数ρ2反映的相关成分次之,所以称u2,v2为第二对典型变量;以此类推.2024/9/5主编：费宇2012.4原始变量与典型相关变量的相关系数记2024/9/5主编：费宇2112.4原始变量与典型相关变量的相关系数2024/9/5主编：费宇2212.4原始变量与典型相关变量的相关系数则上面四个等式可以表示为其中.

2024/9/5主编：费宇2312.4原始变量与典型相关变量的相关系数其中.

其中.

2024/9/5主编：费宇2412.4原始变量与典型相关变量的相关系数其中.

其中.

2024/9/5主编：费宇2512.4原始变量与典型相关变量的相关系数其中.

其中.

2024/9/5主编：费宇2612.5简单相关、复相关和典型相关之间的关系当p=q=1时，x与y之间的（唯一）典型相关就是它们之间的简单相关；当p=1或q=1时，x与y之间的（唯一）典型相关就是它们之间的复相关；可见，复相关是典型相关的一个特例，而简单相关是复相关的一个特例.第一个典型相关系数至少同x（或y）的任一分量与y（或x）的复相关系数一样大，即使所有这些2024/9/5主编：费宇2712.5简单相关、复相关和典型相关之间的关系复相关系数都较小，第一个典型相关系数仍可能很大；同样，从复相关的定义也可以看出，当p=1（或q=1）时，x（或y）与y（或x）之间的复相关系数也不会小于x（或y）与y（或x）的任一分量之间的相关系数，即使所有这些相关系数都较小，复相关系数仍可能很大.2024/9/5主编：费宇2812.6分量的标准化处理一般来说，典型变量是人为定义的，也就是说它没有实质意义.如果使用原始变量，那么典型系数a，b的单位与x和y的单位成比例．而x和y的各分量的单位往往不全相同.我们希望在对各分量作标准化变换之后再作典型相关分析，这样原始变量就有零均值和单位方差，典型变量就没有测量值单位．2024/9/5主编：费宇2912.6分量的标准化处理记，，

，.为的相关矩阵.对x和y的各分量作标准化变换，即令.现在来求和的典型相关变量：2024/9/5主编：费宇3012.6分量的标准化处理于是2024/9/5主编：费宇3112.6分量的标准化处理因为，

，所以其中，.同理2024/9/5主编：费宇3212.6分量的标准化处理其中，.由此可见，为和的第i对典型系数，其第i个典型相关系数仍为ρi，在标准化变换下具有不变性，这一点与主成分分析有所不同．和的第i对典型变量具有零均值，且与x和y的第i对典型变量只相差一个常数．2024/9/5主编：费宇3312.7样本典型相关系数及其

对应典型相关变量的计算前面我们是从变量x与变量y的协方差阵Σ出发考虑x与y的典型相关变量,这称为总体典型相关变量,但在实际例子中一般并不知道Σ,因此通常采用样本协方差阵S代替Σ.由12.6节的分析可知，在大多数情况下，我们在进行典型相关分析时，需将数据标准化，这时样本协方差矩阵S即为样本相关阵.根据样本相关阵计算得到的典型相关变量，称为样本典型相关变量，具体计算过程如下．设容量为n的样本来自正态总体,两组变量的观测值分别记为x=(x1,x2,…,xp

)T和y=(y1,y2,…,yq

)T,不妨设p≤q,则样本数据矩阵为2024/9/5主编：费宇3412.7样本典型相关系数及其

对应典型相关变量的计算2024/9/5主编：费宇35（1）计算样本相关系数（1）计算样本相关系数阵

,并将

剖分为其中,

是第一组变量x的关系数阵,

是第二组变量y的相关系数阵,而

、

(=)为变量x与变量y的相关系数阵.（2）计算典型相关系数及典型变量（2）计算典型相关系数及典型变量设首先求

的特征根，，…，（），并求对应的特征向量

,它是a1,a2,…,am的估计值;再求的特征根对应的特征向量,它是b1,b2,…,bm的估计值.这里

称为样本典型相关系数,而

称为样本典型相关变量.2024/9/5主编：费宇3612.7样本典型相关系数及其

对应典型相关变量的计算（3）记，由式(12.12第一等式)和(12.13第二等式)得2024/9/5主编：费宇3712.8典型相关系数的显著性检验典型相关系数是否显著的不为零,可以通过Bartlett大样本卡方检验来完成.设的m个特征根为,则典型相关系数λ1的显著性检验等价于以下检验H0:λ1=0,H1:λ1≠0.检验统计量为:式中，2024/9/5主编：费宇38式中，12.8典型相关系数的显著性检验在检验水平α下,如果，则拒绝原假设,认为第一对典型变量显著相关.一般,若前j-1个典型相关系数在水平α下是显著的,则当检验第j个典型相关系数的显著性时,检验统计量为

2024/9/539主编：费宇12.8典型相关系数的显著性检验需要指出的是,在实际应用上,通常通过典型相关系数的显著性检验以及典型变量和典型相关系数的实际解释,来确定究竟保留几对典型变量.所求得的典型变量的对数愈少愈容易解释,最好是第一对典型变量就能反映足够多的相关成分,这样只保留一对典型变量便比较理想.2024/9/540主编：费宇12.9被解释样本方差的比例在进行样本典型相关分析时，我们也想了解每组变量提取出的典型变量所能解释的该组样本总方差的比例，由此定量出典型变量所包含的原始信息量的大小．对于经标准化变换后的样本数据，第一组变量的样本总方差为，第二组变量的样本总方差为．称为样本典型相关变量，其中分别是原始变量x,y的标准化结果．2024/9/541主编：费宇12.9被解释样本方差的比例前r对典型相关变量对样本总方差的贡献为：2024/9/542主编：费宇其中，可依据式（12.19）计算．则第一组变量样本方差由前r个典型变量解释的比例为：12.9被解释样本方差的比例同理，第二组变量样本方差由前r个典型变量解释的比例为：2024/9/543主编：费宇其中，可依据式(12.22)计算．例12.1数据文件为exam12.1例12.1(数据文件为exam12.1)康复俱乐部对20名中年人测量了体重(x1)、腰围(x2)、脉搏(x3)三个生理指标和引体向上次数(y1)、仰卧起坐次数(y2)、跳高(y3)三个训练指标，数据详见表12-1，分析生理指标与训练指标的相关性．解:先读取数据，求样本相关系数矩阵．R程序和运行结果如下：9/5/2024主编：费宇44例12.1数据文件为exam12.1表12-1康复俱乐部数据9/5/2024主编：费宇45例12.1数据文件为exam12.19/5/2024主编：费宇46#exam12.1康复俱乐部数据的典型相关分析#打开数据文件exam12.1.xls，选取B1:G21区域，然后复制data12.1<-read.table("clipboard",header=T)#将exam12.1.xls数据读入到data12.1中R=round(cor(data12.1),3);R#求样本相关系数矩阵，保留三位小数x1x2x3y1y2y3x11.0000.870-0.366-0.390-0.493-0.226x20.8701.000-0.353-0.552-0.646-0.191x3-0.366-0.3531.0000.1510.2250.035y1-0.390-0.5520.1511.0000.6960.496y2-0.493-0.6460.2250.6961.0000.669y3-0.226-0.1910.0350.4960.6691.000例12.1数据文件为exam12.1生理指标和训练指标之间的相关性强度中等，其中腰围和仰卧起坐次数的相关系数最大为-0.646；组内较大的是体重和腰围的相关系数为0.87；引体向上次数和起坐次数的相关系数为0.696，仰卧起坐次数和跳高的相关系数为0.669．作典型相关分析，求典型相关系数和对应的典型变量的系数，R程序和运行结果如下：9/5/2024主编：费宇47例12.1数据文件为exam12.19/5/2024主编：费宇48X=scale(data12.1)#对数据进行标准化处理x=X[,1:3]#指定一组变量数据y=X[,4:6]#指定另一组变量数据library(CCA)#载入典型相关分析所用CCA包CCA=cc(x,y)#进行典型相关分析CCA$cor#输出典型相关系数[1]0.79560.20060.0726CCA$xcoef#输出x的典型载荷

[,1][,2][,3]x10.77541.8844-0.1910x2-1.5793-1.18060.5060x30.05910.23111.0508例12.1数据文件为exam12.19/5/2024主编：费宇49CCA$ycoef#输出y的典型载荷

[,1][,2][,3]y10.34950.3755-1.2966y21.0540-0.12351.2368y3-0.7164-1.0621-0.4188因六个变量没有用相同单位测量，这里用标准化后的系数进行分析．第一典型相关系数为0.796，它比生理指标和训练指标两组间的任一其他对的典型相关系数都大．调用相关系数检验脚本进行典型相关系数检验，确定典型变量对数，R程序和运行结果如下：例12.1数据文件为eg12.19/5/2024主编：费宇50source('corcoef_test.R')#调用典型相关系数检验脚本,若该脚本不在R的当前工作路径下,则要将路径设置清晰,如source('C:/ProgramFiles/corcoef_test.R')corcoef_test(r=CCA$cor,n=nrow(x),p=ncol(x),q=ncol(y))#进行典型相关系数检验

rQP[1,]0.795616.25500.0617[2,]0.20060.74500.9457[3,]0.07260.21090.6461检验总体中所有典型相关系数均为0的零假设时概率水平为0.062，故在（或）的显著性水平下，拒绝所有典型相关系数均为0的假设.也就是至少有一对典型相关是显著的.从后面的检验结果例12.1数据文件为exam12.19/5/2024主编：费宇51可知，只有一对典型相关是显著的．结合前面输出的典型相关载荷结果可知，生理指标的第一典型变量为：它近似地是腰围和体重的加权和，在腰围上的权数更大些，在脉搏上的权数近似为0.来自训练指标的第一典型变量为：其在仰卧起坐次数上的权数最大．这对典型变量主要是反映腰围和仰卧起坐的负相关关系．例12.1数据文件为exam12.19/5/2024主编：费宇52输出原始变量和典型变量的相关系数.R程序和运行结果如下：CCA$scores$corr.X.xscores#输出第一组典型变量与X组原始变量之间的相关系数

[,1][,2][,3]x1-0.62060.7724-0.1350x2-0.92540.3777-0.0310x30.3328-0.04150.9421CCA$scores$corr.Y.xscores#输出第一组典型变量与Y组原始变量之间的相关系数

[,1][,2][,3]y10.5789-0.0475-0.0467y20.6506-0.11490.0040例12.1数据文件为exam12.19/5/2024主编：费宇53y30.1290-0.1923-0.0170CCA$scores$corr.X.yscores#输出第二组典型变量与X组原始变量之间的相关系数

[,1][,2][,3]x1-0.49380.15498-0.0098x2-0.73630.07578-0.0022x30.2648-0.00830.0684CCA$scores$corr.Y.yscores#输出第二组典型变量与Y组原始变量之间的相关系数

[,1][,2][,3]y10.7276-0.2370-0.6438y20.8177-0.57300.0544y30.1622-0.9586-0.2339例12.1数据文件为exam12.19/5/2024主编：费宇54整理后得表12-2．表12-2原始变量与第一对典型变量的相关系数由表12-2可知来自生理指标的第一典型变量u1与腰围的相关系数为-0.925，与体重的相关系数为-0.621，它们都是负的．但在典型变量中体重的载荷为正（0.775），即体重在中的载荷和它与的相关系数反号．来自例12.1数据文件为exam12.19/5/2024主编：费宇55训练指标的第一典型变量与三个训练指标的相关系数都是正数，其中跳高在中的载荷（-0.716）和它与的相关系数（0.1622）反号；因此，体重和跳高在这组变量中分别是一个校正（或抑制）变量．一个变量具有同典型变量的相关系数相反符号的载荷似乎是矛盾的.为了理解这是怎样发生的，考虑简单的情况：用多元回归方法由腰围和体重来预测仰卧起坐次数.一般来说，胖的人比瘦的人仰卧起坐次数少，这似乎是有道理的.假定这组样本中没有非常高的人，于是腰围和体重之间的相关系数（0.87）是很大的．检验肥胖同自变量之间的相关性：例12.1数据文件为exam12.19/5/2024主编：费宇56腰围大的人倾向于比腰围小的人胖，因此腰围与仰卧起坐次数之间的相关为负相关．体重大的人倾向于比体重小的人胖，于是体重与起坐起坐次数之间的相关为负相关．固定体重的值，腰围大的人倾向于较强壮和较胖，于是腰围的多元回归系数应是负的．固定腰围的值，体重大的人倾向于比较高和比较瘦，因此体重的多元回归系数应为正的，它与体重和仰卧起坐次数间的相关反号．因此，第一典型相关一般解释为以体重和跳高作为例12.1数据文件为exam12.19/5/2024主编：费宇57抑制变量来提高腰围和起坐次数之间的相关性，但样本的大小对于得出确定的结论还不够大．计算典型变量解释原变量方差的比例，R程序和运行结果如下：apply(CCA$scores$corr.X.xscores,2,function(x){mean(x^2)})#第一组典型变量解释原第一组变量方差的比例[1]0.45080.24700.3022apply(CCA$scores$corr.Y.xscores,2,function(x){mean(x^2)})#第一组典型变量解释原第二组变量方差的比例[1]0.25840.01750.0008例12.1数据文件为exam12.19/5/2024主编：费宇58第一对典型变量中u1解释生理指标的标准方差的比例为0.451，第一对典型变量中v1解释训练指标的标准方差的比例为0.408，但两者都不能很好地全面预测对应的那组变量.因为来自生理指标的标准方差被对方第一个典型变量v1解释的方差比例为0.285，而来自训练指标的标准方差被对方第一典型变量u1解apply(CCA$scores$corr.X.yscores,2,function(x){mean(x^2)})#第二组典型变量解释原第一组变量方差的比例[1]0.28540.00990.0016apply(CCA$scores$corr.Y.yscores,2,function(x){mean(x^2)})#第二组典型变量解释原第二组变量方差的比例[1]0.40810.43450.1574例12.1数据文件为exam12.19/5/2024主编：费宇59释的方差比例为0.258．计算得分，并绘制得分等值平面图.R程序如下：u<-as.matrix(x)%*%CCA$xcoef#计算得分v<-as.matrix(y)%*%CCA$ycoef#计算得分plot(u[,1],v[,1],xlab="u1",ylab="v1")#绘制第一对典型变量得分的散点图,x轴名称为u1,y轴名称为v1，见图12-1abline(0,1)#在散点图上添加一条y等于x的线,以查看散点分布情况图12-1康复俱乐部数据第一对典型相关得分等值平面图例12.1数据文件为exam12.19/5/2024主编：费宇60例12.1数据文件为exam12.19/5/2024主编：费宇61通过作第一对典型相关变量得分等值平面图可以看出，散点在一条近似直线上，虽然有偏离情况发生，但总体还是呈现出了线性相关关系.综合来看，生理指标与训练指标之间的关系虽有波动，但从整体来看较为明显．12.10案例分析及R操作案例12.1(数据文件为case12.1)表12-3给出了2008-2016年我国科技活动和经济发展的部分代表指标.其中，科技活动指标：x1为R&D人员全时当量（单位：万人年），x2为R&D经费支出（单位：亿元），x3为R&D项目(课题)数（单位：项），x4为发表科技论文数（单位：篇），x5为专利申请授权数（单位：件）；经济发展指标：y1为国内生产总值（单位：亿元），y2为城镇居民家庭人均可支配收入（单位：元），y3为农村居民家2024/9/562主编：费宇案例12.1数据文件为case12.1庭人均纯收入（单位：元）；利用这些数据进行典型相关分析来分析我国科技活动和经济发展的关系．2024/9/563主编：费宇案例12.1数据文件为case12.1解:先读取数据，求样本相关系数矩阵．R程序和运行结果如下：2024/9/564主编：费宇#case12.1我国科技活动和经济发展的典型相关分析#打开数据文件case12.1.xls,选取B1:I10区域,然后复制case12.1<-read.table("clipboard",header=T)#将case12.1.xls数据读入到case12.1中R=round(cor(case12.1),3);R#求样本相关系数矩阵,保留三位小数案例12.1数据文件为case12.12024/9/565主编：费宇x1x2x3x4x5y1y2y3x11.0000.9880.9870.9950.9690.9900.9840.979x20.9881.0000.9990.9840.9920.9950.9970.995x30.9870.9991.0000.9810.9920.9930.9950.993x40.9950.9840.9811.0000.9690.9850.9820.979x50.9690.9920.9920.9691.0000.9870.9950.997y10.9900.9950.9930.9850.9871.0000.9980.993y20.9840.9970.9950.9820.9950.9981.0000.997y30.9790.9950.9930.9790.9970.9930.9971.000科技活动指标和经济发展指标之间的相关性很强，组内相关性也很强.作典型相关分析，求典型相关系数和对应的典型变量的系数，R程序和运行结果如下：案例12.1数据文件为case12.12024/9/566主编：费宇X=scale(case12.1)#对数据进行标准化处理x=X[,1:5]#指定一组变量数据y=X[,6:8]#指定另一组变量数据library(CCA)#载入作典型相关分析所用CCA包CCAc12.2=cc(x,y)#进行典型相关分析CCAc12.2$cor#输出典型相关系数[1]0.999980.867060.29489CCAc12.2$xcoef#输出x的典型载荷

[,1][,2][,3]x10.1188-7.5900.1671x2-0.7634-1.44013.9738x30.65382.655-19.2193x4-0.16523.6411.6708x5-0.84572.6453.3987案例12.1数据文件为case12.12024/9/567主编：费宇CCAc12.2$ycoef#输出y的典型载荷

[,1][,2][,3]y10.6016-13.1187.092y2-0.89369.993-19.181y3-0.70583.03612.106因六个变量没有用相同单位测量，这里用标准化后的系数进行分析.第一典型相关系数为0.99998，它比科技活动指标和经济发展指标间的任一相关系数都大．调用相关系数检验脚本进行典型相关系数检验，确定典型变量对数，R程序和运行结果如下：案例12.1数据文件为case12.12024/9/568主编：费宇source('corcoef_test.R')#调用典型相关系数检验脚本,若该脚本不在R的当前工作路径下,则要将路径设置清晰,如source('C:/ProgramFiles/corcoef_test.R')corcoef_test(r=CCAc12.2$cor,n=nrow(x),p=ncol(x),q=ncol(y))#进行典型相关系数检验

rQP[1,]0.9999840.491410.00038[2,]0.867065.195750.73646[3,]0.294890.348460.95067检验总体中所有典型相关系数均为0的零假设时概率水平远小于,否定所有典型相关系数均为0的假设，也就是至少有一对典型相关是显著案例12.1数据文件为case12.12024/9/569主编：费宇的；典型相关系数检验p值的第二个值为0.736、第三个值为0.951，因此在显著性水平为0.05的情况下只有一对典型相关是显著的．结合前面输出的典型相关载荷结果来看，科技活动指标的第一典型变量为：它近似地是专利申请授权数、R&D经费支出和R&D项目（课题）数的加权和.在专利申请授权数上的权数最大，其次是R&D经费支出，在R&D项目（课题）数上的权数也较大．案例12.1数据文件为case12.12024/9/570主编：费宇来自经济发展指标的第一典型变量为：它在城镇居民家庭人均可支配收入上的权数最大，其次为农村居民家庭人均纯收入．输出原始变量和典型变量的相关系数，R程序和运行结果如下：CCAc12.2$scores$corr.X.xscores#输出第一组典型变量与X组原始变量之间的相关系数

[,1][,2][,3]x1-0.97482-0.20944-0.02677x2-0.99413-0.08663-0.04131案例12.1数据文件为case12.12024/9/571主编：费宇x3-0.99275-0.07465-0.08732x4-0.97608-0.157720.01888x5-0.998780.02189-0.03766CCAc12.2$scores$corr.Y.xscores#输出第一组典型变量与y组原始变量之间的相关系数

[,1][,2][,3]y1-0.99059-0.11673-0.00708y2-0.99682-0.05889-0.01212y3-0.99907-0.024930.00931CCAc12.2$scores$corr.X.yscores#输出第二组典型变量与X组原始变量之间的相关系数

[,1][,2][,3]x1-0.97480-0.18160-0.00790案例12.1数据文件为case12.12024/9/572主编：费宇x2-0.99411-0.07511-0.01218x3-0.99273-0.06472-0.02575x4-0.97606-0.136750.00557x5-0.998750.01898-0.01111CCAc12.2$scores$corr.Y.yscores#输出第二组典型变量与y组原始变量之间的相关系数

[,1][,2][,3]y1-0.99061-0.13463-0.02399y2-0.99684-0.06792-0.04108y3-0.99909-0.028750.03156整理后得到表12-4.案例12.1数据文件为case12.1表12-4原始变量与第一对典型变量的相关系数2024/9/573主编：费宇来自科技活动指标的第一典型变量u1与R&D经费支出、发表科技论文数、专利申请授权数的相关系数分别为-0.97482、-0.97608、-0.99878，u1与R&D人员全时当量、R&D项目（课题）数的相关系数分别为-0.99413、案例12.1数据文件为case12.12024/9/574主编：费宇-0.99275，它们都是负的，因此R&D人员全时当量、R&D项目（课题）数是抑制变量，其含义是它们在中的载荷（0.119，0.654）和它们与的相关系数（-0.97482，-0.99275）反号．来自经济发展指标的第一典型变量v1与三个经济发展指标的相关系数是负值，因国内生产总值在中的系数和它与的相关系数是反号，故国内生产总值也是一个抑制变量．案例12.1数据文件为case12.12024/9/575主编：费宇计算典型变量解释原变量方差的比例，第一对典型变量能很好地全面预测对应的那组变量，来自科技活动指标的标准方差被第一个典型变量u1解释的方差比例为0.9749，第一对典型变量v1解释经济发展指标的标准方差的比例为0.9911；而来自经济发展指标的标准方差被第一典型变量v1解释的方差比例为0.9748，经济发展指标的标准方差被对方第一典型变量u1解释的方差比例为0.991．R程序和运行结果如下：案例12.1数据文件为case12.12024/9/576主编：费宇apply(CCAc12.2$scores$corr.X.xscores,2,function(x){mean(x^2)})#第一组典型变量解释原第一组变量方差的比例[1]0.97490.01650.0024apply(CCAc12.2$scores$corr.Y.xscores,2,function(x){mean(x^2)})#第一组典型变量解释原第二组变量方差的比例[1]0.99100.00590.000009apply(CCAc12.2$scores$corr.X.yscores,2,function(x){mean(x^2)})#第二组典型变量解释原第一组变量方差的比例[1]0.97480.01240.0002apply(CCAc12.2$scores$corr.Y.yscores,2,function(x){mean(x^2)})#第二组典型变量解释原第二组变量方差的比例[1]0.99110.00790.0011案例12.1数据文件为case12.12024/9/577主编：费宇计算得分，并绘制得分等值平面图.R程序如下：u<-as.matrix(x)%*%CCAc12.2$xcoef#计算得分v<-as.matrix(y)%*%CCAc12.2$ycoef#计算得分plot(u[,1],v[,1],xlab="u1",ylab="v1")#绘制第一对典型变量得分的散点图,x轴名称为u1,y轴名称为v1，见图12-2abline(0,1)#在散点图上添加一条y等于x的线,以查看散点分布情况图12-2科技活动和经济发展数据第一对典型相关得分等值平面图案例12.1数据文件为case12.19/5/2024主编：费宇78案例12.1数据文件为case12.12024/9/579主编：费宇由得分等值平面图图12-2可以看出，第一对典型相关变量得分散点在一条直线上分布，两者之间呈高度线性相关关系，散点图上没有离开群体的差异点.这说明我国科技活动与经济发展之间的关系很稳定，整体波动平稳．2024/9/5主编：费宇主编：费宇，鲁筠中国人民大学出版社，2024年第13章多维标度分析2024/9/5

82多维标度法（MultidimensionalScaling，MDS）是一种维数缩减方法，把高维的数据点映射到一个低维空间上，然后在低维空间中以较高的相似度重新展示这些点的数据结构，并由此对原始高维数据进行统计分析．多维标度法起源于上世纪40年代的心理测度.1958年

Torgerson

正式提出了这一方法．目前应用广泛，内容丰富，方法较多．两类方法：度量分析法，非度量分析法.第13章多维标度分析13.1

多维标度法的基本思想13.2

古典多维标度法

多维标度法的几个基本概念；

已知距离矩阵时CMDS解；

已知相似系数矩阵时CMDS解；13.3

案例分析与R实现2024/9/5

8313.1

多维标度法的基本思想2024/9/5

84多维标度法是用低维空间

Rk

(k<p)中的n个点去重新标度和展示高维空间Rp中相应的n个点，将高维空间中的研究对象(样本或变量)简化到低维空间中进行定位、归类和分析且有效保留研究对象间原始关系的多元数据分析方法．是一种维数缩减方法，主要思想是在降维的同时让新得到的n个低维点与原来的n高维个点保持较高的相似度(如位置关系、距离、类别等)，使降维后重新标度的数据结构发生的“形变”尽量小．2024/9/5

85多维标度分析与主成分分析的异同：

共同点：均先降维，再做简明有效的分析；

均依据特定矩阵的特征值的累积贡献率来

决定降维后的维数．

不同点：主成分法按包含信息大小选取主成分；

多维标度法按标度前后距离矩阵尽量接近

或尽量相似的原则来构造拟合构造点．13.2.1

多维标度法的几个基本概念2024/9/5

8613.2

古典多维标度法古典多维标度(CMDS)解的定义对于距离阵，多维标度法的目的是要寻找较小的

k和Rk中的n个点,记，表示

xi

与xj

在Rk中的欧氏距离，使得与在某种意义下尽量接近,记称为的一个古典多维标度(CMDS)解,称xi

为的一个拟合构造点，亦称为的拟合构图，称为的拟合距离阵．2024/9/587特别，当时，称xi

为的一个构造点，称为的构图，注意的构图不唯一2024/9/5

88注：欧式距离阵也称为欧式型距离阵或Euclid距离阵

如何判断一个距离阵

是否为欧氏型距离阵呢？几个相关矩阵的构造：2024/9/5

89

距离阵为欧氏型距离阵的充要条件：

从欧氏距离阵D出发得到构图

X的步骤：见下面例13.1．2024/9/5

90例13.1(数据文件exam13.1)给定距离矩阵D如下，问它是否是欧式型矩阵？并求其的古典多维标度解（CMDS）

X及拟合构造点．2024/9/5

91由B的特征值进行判断并求X

解：在Excel中由矩阵D计算出B很方便：将D的上三角部分补齐输入单元格区域A1:D4（如图13-1

）；然后在A6单元格内输入“=-(A1^2)/2”后拖放填充至D9得到矩阵A；接着如图13-1在A的左方、下方输入公式并拖放填充计算出行均值、列均值和总均值；最后在单元格A12内输入公式“=A6-$E6-A$10+$E$10”后拖放填充至D15得到矩阵B．图13-1在Excel中由矩阵D算出矩阵B（A12:D15）的简单过程2024/9/5

92

最后在R中计算B的所有特征值来判断D是否是欧氏型距离阵，并求出CMDS解

X及拟合构造点．2024/9/5

93>B=read.table("clipboard",header=F)

#读入矩阵B（先复制单元格区域A12:D15）>eig<-eigen(B);eig

#求B的特征值特征向量并显示eigen()decomposition$values[1]70.6013138636.6152425180.795952899-0.002509279

B的四个特征值均大于等于0（第四个特征值很小可视为0），由定理13.1知D为欧氏型距离阵．为求D的CMDS解

X及拟合构造点（实际上可近似看作构造点），在R中输入如下命令:2024/9/5

94>

D=read.table("clipboard",header=F)

#读入矩阵D(先复制对应单元格区域A1:D4)>D13.1=cmdscale(D,k=2,eig=T)；D13.1#使用stats中的cmdscale函数，k取为2$points

[,1]

[,2][1,]2.526493-0.1531671[2,]

-7.110732-0.4721962[3,]3.545859-1.4437017[4,]1.038381

2.0690650$eig[1]7.059702e+016.611734e+007.937158e-012.428205e-15>sum(abs(D13.1$eig[1:2]))/sum(abs(D13.1$eig))

#计算a1.2[1]0.9898245>sum((D13.1$eig[1:2])^2)/sum((D13.1$eig)^2)

#计算a2.2[1]0.9998747D的CMDS解：2024/9/595前两个特征值的累积绝对贡献率和累积平方贡献率均超过98%，说明k取为2是适当的．四个构造点在坐标分别为(2.526,-0.153),(-7.111,-0.472),(3.546,-1.444),(1.038,2.069)．拟合距离阵（可与原始距离阵D对比）和拟合构图如下：图13-2

距离阵D的拟合构图13.2.2

已知距离矩阵时CMDS解上面计算CMDS解的过程在R软件中可使用stats包中的cmdscale()函数来实现，也可以使用MASS包中处理非度量MDS问题的isoMDS()函数来实现，但cmdscale()函数的好处是可以同时计算出的特征值和特征向量以及两个累积贡献率a1.k和a2.k的值．

2024/9/5

96

例13.2对表13-1给出的我国六个城市间的距离矩阵D，利用R软件stats包中的cmdscale函数求

的CMDS解，给出拟合构图

及拟合构造点．2024/9/597

北京济南青岛郑州上海南京北京0

济南4390

青岛6683620

郑州7144437720

上海12598867769840

南京10656266177103220表13-1我国六个城市间的道路距离（单位:公里）2024/9/5

98解在R中的程序为:>setwd("C:/data")#设定工作路径>exam13.2<-read.csv("exam13.2.csv",header=T)#将exam13.2数据读入>d13.2=exam13.2[,-1]#先去掉exam13.2第一列样本名称>rownames(d13.2)=exam13.2[,1]#为d13.2的行重新命名>D13.2=cmdscale(d13.2,k=2,eig=T);D13.2#使用基本包stats中的cmdscale函数,k取为2由R计算结果可见，矩阵的六个特征值分别为：1051894，311141，59859，10289，0，-11999最后一个特征值为负，表明距离矩阵不是欧氏型距离阵．a1.k=94.3%，a2.k=99.7%，故k=2就可以了．由前两个特征向量可得六个拟合构造点分别为：(612.2,119.4)，(218.2,11.8)，(38.0,319.8)(193.7,430.2)，(646.1,57.4)，(416.0,78.3)2024/9/5

99>x=D13.2$points[,1];>y=D13.2$points[,2]>plot(x,y,xlim=c(-700,800),ylim=c(-300,600))#根据两个特征向量的分量大小绘散点图>text(x,y,labels=s(d13.2),adj=c(0,-0.5),cex=0.8)#将拟合点用行名标出再画六个城市距离矩阵的拟合构图，并用中文标明（注意：拟合构图主要表示六个城市间的相对距离，和各城市在地图上的实际位置可能不一致）图13-3

我国六城市距离阵的拟合构图2024/9/5

1002024/9/5101

易计算出六个拟合构造点在R2中的欧氏距离阵，将它们与表13-1中城市间的原始距离数据进行对比，可见大多数距离数据拟合较好，有6个城市距离相差公里数为个位数，只有两个城市距离相差在50~60公里．>D1<-dist(D13.2$points[,1:2],method="euclidean",diag=T,p=2)>D1

北京

济南

青岛

郑州

上海

南京北京

0.0000济南

408.48960.0000青岛

608.1445356.81290.0000郑州

690.8417442.6457765.94660.0000上海1259.8374865.4658732.6946971.06630.0000南京1047.0937640.5733603.8211704.0036267.07110.000013.2.3

已知相似系数矩阵时的CMDS解2024/9/5102

例13.3对55个国家和地区的男子径赛记录作统计，每位运动员记录8项指标：100米、200米、400米、800米、1500米、5000米、10000米和马拉松．这8项指标的相关矩阵C如表13-2所示．求C的CMDS解，并给出拟合构图及拟合构造点

(数据文件为exam13.3)．

2024/9/5

1032024/9/5

104

100米200米400米800米1500米5000米10000米马拉松100米10.9230.8410.7560.7000.6190.6330.520200米0.92310.8510.8070.7750.6950.6970.596400米0.8410.85110.8700.8350.7790.7870.705800米0.7560.8070.87010.9180.8640.8690.8061500米0.7000.7750.8350.91810.9280.9350.8665000米0.6190.6950.7790.8640.92810.9750.93210000米0.6330.6970.7870.8690.9350.97510.943马拉松0.5200.5960.7050.8060.8660.9320.9431表13－2运动员径赛8项指标的相关系数矩阵2024/9/5

105>setwd("C:/data")#设定工作路径>eg13.3<-read.csv("example13.3.csv",header=T)#将数据读入>c13.3=eg13.3[,-1]#exam13.3的第一列为样本名称不是数值，先去掉>d13.3=round(sqrt(2-2*c13.3),3)#相似阵转换成广义距离阵,取三位小数>rownames(d13.3)=exam13.3[,1]>d13.3

X100米X200米X400米

X800米X1500米X5000米X10000米马拉松100米0.0000.3920.5640.6990.7750.8730.8570.980200米0.3920.0000.5460.6210.6710.7810.7780.899400米0.5640.5460.000

0.5100.5740.6650.6530.768800米0.6990.6210.5100.0000.4050.522

0.5120.6231500米0.7750.6710.574

0.4050.0000.3790.3610.5185000米0.8730.7810.665

0.5220.3790.0000.2240.36910000米0.8570.7780.6530.5120.3610.2240.0000.338马拉松0.9800.8990.7680.6230.5180.3690.3380.0002024/9/5

106>D13.3=cmdscale(d13.3,k=2,eig=T),D13.3

#k取为2，给出特征向量和特征值$points[,1][,2]100米-0.5273079050-0.4241746300.12620103400米-0.219741318-0.21411107800米0.003249949-0.227485801500米0.146859676-0.084591405000米0.3046913650.0667815610000米0.2994955630.06799955马拉松0.4169272990.12579130$eig[1]8.842051e-011.650194e-011.258240e-018.178272e-026.736816e-02[6]4.726735e-022.266967e-023.045637e-18……>sum(abs(D13.3$eig[1:2]))/sum(abs(D13.3$eig))#计算a1.2[1]