《多元统计分析-基于R（第3版）》课件 第13章-多维标度分析

主编：费宇中国人民大学出版社第13章多维标度分析2024/7/31

2多维标度法（MultidimensionalScaling，MDS）是一种维数缩减方法，把高维的数据点映射到一个低维空间上，然后在低维空间中以较高的相似度重新展示这些点的数据结构，并由此对原始高维数据进行统计分析．多维标度法起源于上世纪40年代的心理测度.1958年Torgerson正式提出了这一方法．目前应用广泛，内容丰富，方法较多．两类方法：度量分析法，非度量分析法.第13章多维标度分析13.1

多维标度法的基本思想13.2

古典多维标度法

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

已知距离矩阵时CMDS解；

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

案例分析与R实现2024/7/31

313.1

多维标度法的基本思想2024/7/31

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

Rk

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

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

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

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

决定降维后的维数．

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

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

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

多维标度法的几个基本概念2024/7/31

613.2

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

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

xi

与xj

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

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

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

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

如何判断一个距离阵

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

9

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

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

X的步骤：见下面例13.1．2024/7/31

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

X及拟合构造点．2024/7/31

11由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/7/31

12

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

X及拟合构造点．2024/7/31

13>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/7/31

14>

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/7/3115前两个特征值的累积绝对贡献率和累积平方贡献率均超过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/7/31

16

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

的CMDS解，给出拟合构图

及拟合构造点．2024/7/3117

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

济南4390

青岛6683620

郑州7144437720

上海12598867769840

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

18解

在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/7/31

19>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/7/31

202024/7/3121

易计算出六个拟合构造点在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/7/3122

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

(数据文件为exam13.3)．

2024/7/31

232024/7/31

24

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/7/31

25>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/7/31

26>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]0.7525982>sum((D13.3$eig[1:2])^2)/sum((D13.3$eig)^2)#计算a2.2[1]0.9644675由R输出结果知B的八个特征值均非负，说明所得的广义距离阵为欧氏型距离阵．两个累积贡献率分别为75.26%和96.45%，权衡后认为取k=2较适当．2024/7/31

27将这八个拟合构造点画出并用中文命名，得到八项指标相似系数矩阵的古典拟合构图>x=D13.3$points[,1];y=D13.3$points[,2]>plot(x,y,xlim=c(-0.6,0.7),ylim=c(-0.25,0.2))#绘制拟合图>text(x,y,labels=s(d13.3),adj=c(-0.2,-0.2),cex=0.8)#将拟合点用行名标出>abline(h=0,v=0,lty=3)#用虚线划分四个象限

从图13－4可以直观地看出，“100米”和“200米”两项指标靠得很近，“5000米”和“10000米”两项指标几乎重叠在一起，说明这两对指标各自相似性都很强，这与表13－2给出的相关系数大小是一致的，它们分别表明爆发力和耐力的相似性．说明重新降维标度后指标仍保持强相似性．2024/7/3128图13-4运动员八项径赛指标相似系数矩阵的古典拟合构图13.3案例分析与R实现

案例13.1表13-3给出了2020年我国第七次人口普查数据集中我国各地城市3岁及以上人口中接受过各类教育人群的分布情况数据．试用多维标度法对其进行统计分析，并对分析结果的实际意义进行解释(数据文件为case13.1)．2024/7/31

29

未上过学学前教育小学初中高中大专本科硕士博士

北京199060496110166071634858003135889249246544167451206254219589

天津1554313117831455412314089421544171384464182209722367829222

河北2570628644063475010664064746314793019961236582318972617925

山西1316204643861795519399297027419211799176170797215374316737

内蒙古173408284087141756127877361813584137671312229271067038628

辽宁2753135933763120335955255748671583075020323440229979532491

吉林1162582313711313922345714225288751097458122887310416015901…………………………………………

青海8449271159423099560095345708281655275581163781717

宁夏89249109750497759841008525966395997390576292452676

新疆14741137387817114632494243170027714531081097813660525816表13－3我国第七次人口普查各类受教育人群的分布情况数据本案例我们采用MASS包中的isoMDS函数来实现2024/7/31

30>setwd("C:/data")#设定工作路径>case13.1<-read.csv("case13.1.csv",header=T)#将case13.1.csv数据读入>c13.1=case13.1[,-1]#case13.1的第一列不是数值先去掉，命名为c13.1>rowna