第十二讲 双标图biplot在同一个二维坐标系中同时标识样本个体和变量
第十二讲 双标图biplot 在同一个二维坐标系中同 时标识样本个体和变量
recap假设协方差矩阵var(x）=Z的谱分解Z=VAVT·=diag()..≥为的特征根V.V,为相应的单位长度相互正交的特征向量；·V=(vi.,)载荷矩阵， VTV=VT=I,总体PCA：·主成分：y=VTxeRP,var(y)=△, =vfx-2vx, var(y,)=,j= ,., p..载荷v, =cov(x,y,)/ a,Vy =cov(xi,y,)/ j·取前k个主成分使得方差累计贡献率（+..+)/+...+a)>CR命令：>princomp(covmat=sigma）#对方差矩阵sigma进行总体PCA分析2
2 recap ( ,., ) . ,., ; ( ,., ) . , var( ) 1 1 1 1 2 p p p p p V V V VV I diag V V T T T 载荷矩阵， 为相应的单位长度相互正交的特征向量 ： 为 的特征根 假设协方差矩阵 的谱分解 v v v v x . )/ . . cov( , )/ cov( , )/ var( ) , 1,., . , var( ) . PCA : 1 1 1 k C y v x y y v x y j p V R k p j j j ij i j j j j p i j j ij i p 取前 个主成分使得方差累计贡献率（ （ ） 载荷 ， ， 主成分： 总体 v x v x y x y T T R命令：> princomp(covmat=sigma) # 对方差矩阵sigma进行总 体PCA分析
样本PCA：样本x..X,ERP,样本协方差矩阵S的谱分解：S=VAVT,A= diag(...a,), a...≥ap, VTV=VVT =Ip, Vpxp =(Vi..,V,).X,ERP的主成分变换：y,=VT(x,-x),(y,=VTx,也可以定义为主成分，与y,=VT(x,-x)相差一个平移常数)主成分矩阵：Yxp=XV或Y=XV取前k个主成分使得(+..+)(+..+,)>0.8PC散点图：在二维散点图上画出前k个主成分的两两散点图。R命令：>princomp(x）#对数据x进行PCA分析> prcomp(x)#不充许总体PCA，充许p>nm
3 diag( ,., ) . ( ,., ). ,., , PCA 1 1 1 1 p p p p p p p n V V VV I V S V V R S v v x x ， ， ， ， 样本 样本协方差矩阵 的谱分解： 样本 ： T T T ( . )/( . ) 0.8. ( ) ) ( ) 1 1 k p n p c i i i i i i p i k Y X V Y XV V V R V 取前 个主成分使得 主成分矩阵： 或 （ 也可以定义为主成分，与 相差一个平移常数 的主成分变换： ， y x y x x x y x x T T T 在二维散点图上画出前 个主成分的两两散点图。 散点图： k PC R命令：> princomp(x) # 对数据x进行PCA分析 > prcomp(x) #不允许总体PCA，允许𝑝 > 𝑛
样本PCA例子第11讲的几个例子都是总体PCA（只有协方差矩阵或相关系数，而没有完整的数据)，R函数为princomp(covmat=)总体PCA只能得到载荷、特征根，因为没有原始数据X无法计算主成分。样本PCA是指原始数据X可用的情形，样本PCA的R函数：princomp(X)：总体PCA和样本PCA，只能处理n>p情形；prcomp：允许n≤p,不能做总体PCASvd(X)：奇异值分解，可以做PCA。注意：princomp和prcomp的载荷符号相反。样本PCA输出结果与总体PCA基本相同，因为有原始数据X，可以计算每个样本的主成分（princomp输出结果中的score，prcomp输出的x）。4
4 princomp(X): 总体PCA和样本PCA , 只能处理 𝑛 > 𝑝 情形 ； prcomp: 允许𝑛 ≤ 𝑝, 不能做总体PCA. Svd(X)：奇异值分解，可以做PCA。 注意：princomp和prcomp的载荷符号相反。 样本PCA例子 第11讲的几个例子都是总体PCA (只有协方差矩阵或相关系数，而没有 完整的数据)，R函数为 princomp(covmat= ) 总体PCA只能得到载荷、特征根，因为没有原始数据 𝑋 无法计算主成分。 样本PCA是指原始数据𝑋可用的情形, 样本PCA的R 函数: 样本PCA输出结果与总体PCA基本相同，因为有原始数据𝑋，可以计算每 个样本的主成分（princomp输出结果中的score, prcomp输出的x）
例1（样本PCA，第一讲例2的一部分数据）800个人的1000个基因位点的基因型数据，取值0,1,2(aa,Aa,AA)。前两个主成分基本能将各个种族区分开：AMREASSASto#代码：gene=read.table("http://staff.ustc.edu.cn/~ynyang/vector/data/genedatal.txt")race=gene[,1#第一列是种族（辅助信息，1维，不需要降维）x-gene[,-1]#基因型x=as.matrix(x)#n=800,p=1000mypca=prcomp(x）#对基因型数据进行PCA分析pc=mypcasx#主成分矩阵800x800plot(pc[,1:2])#前两个主成分的散点图，数据大概有四个类(cluster)points(pc[1,2],col=as.factor(race)）#标识种族信息，四个类分别对应AFR,EAS,EUR(SAS+AMR）legend(5,-4,legend=c("AFR","AMR","EAS""EUR","SAS"),col=1:5,pch=rep(1,5))
5 例1 （样本PCA，第一讲例2的一部分数据）800个人的1000个基因位点的 基因型数据，取值0,1,2 (aa,Aa,AA)。前两个主成分基本能将各个种族区 分开： #代码： gene=read.table("http://staff.ustc.edu.cn/~ynyang/vector/data/genedata1.txt") race=gene[,1] #第一列是种族 （辅助信息，1维，不需要降维） x=gene[,-1] #基因型 x=as.matrix(x) #n=800,p=1000 mypca=prcomp(x) # 对基因型数据进行PCA分析 pc=mypca$x #主成分矩阵 800x800 plot(pc[,1:2]) # 前两个主成分的散点图，数据大概有四个类(cluster) points(pc[1,2], col=as.factor(race)) #标识种族信息，四个类分别对应AFR,EAS,EUR,(SAS+AMR) legend(5,-4,legend=c("AFR","AMR","EAS","EUR","SAS"), col=1:5, pch=rep(1,5))
例2.数据集temperature.csv给出了欧洲23个国家的首都城市的气温数据以及其它相关信息。变量为12个月份的月内平均气温。其它辅助信息包括，Annual,Amplitude，地理信息（经纬度、地区）。数据：temperature.csv(http://staff.ustc.edu.cn/~ynyang/vector/data/temperature.csv)来源：FrancoisHusson，SebastienLe,JeromePages，ExploratoryMultivariateAnalysisbyExampleUsingR,Chapman&Hall2017变量Jan-DecAnnualAmplitudeLatitudeLongitudeArea含义纬度经度区域1-12月平年度平平均月内均气温均气温温度极差国家首都城市国家位置首都城市位置首都城市国家WWAmsterdam（阿姆斯特丹）希腊S德国何兰Athens（雅典）Berlin（柏林）EWN比利时Budapest（布达佩斯）匈牙利Copenhagen（哥本哈根）丹麦Brussels(布鲁塞尔）NNE爱尔兰芬兰乌克兰Dublin（都柏林）Helsinki（林尔辛基）Kiev（基辅）NEs波兰旧都葡萄牙英国Krakow（克拉科大）London（伦敦）Lisbon（里斯本）sEE俄罗斯Madrid（马德里）西班牙Minsk（明斯克）白俄罗斯Moscow（莫斯科）WEN法国挪威Paris(巴黎)捷克Oslo（奥斯陆）Prague（布拉格）NsSZ冰岛意大利波黑Reykjavik（雷克雅未克）Rome（罗马）Sarajevo（拉热窝）E保加利亚瑞典Stockholm（斯德哥尔摩）Sofia（索菲亚）Annual和Amplitude是从温度记录原始数据汇总得到的，因此不作为感兴趣的变量包含进PCA.我们将利用它们以及地理位置信息对PCA得到的结果进行解释或评价。6
6 例2. 数据集temperature.csv 给出了欧洲23个国家的首都城市的气温数 据以及其它相关信息。变量为12个月份的月内平均气温。其它辅助信息 包括，Annual, Amplitude, 地理信息(经纬度、地区）。 变量 Jan-Dec Annual Amplitude Latitude Longitude Area 含义 1-12月平 均气温 年度平 均气温 平均月内 温度极差 纬度 经度 区域 数据：temperature.csv (http://staff.ustc.edu.cn/~ynyang/vector/data/temperature.csv) 来源： François Husson，Sébastien Lê，Jérôme Pagès，Exploratory Multivariate Analysis by Example Using R，Chapman &Hall 2017 Annual和 Amplitude是从温度记录原始数据汇总得到的，因此不作为感兴趣的变量包 含进PCA. 我们将利用它们以及地理位置信息对PCA得到的结果进行解释或评价
countries and capitals of the EuropeNorwegianSeaICELANDSWEDENFaroeNORWAYReykavikIslandsFINLANDTonsDend三HelsinkiPOsto*Socbolm*TanyESTaRig*LATNorthCRLITHSea*MoscowCopRUSSAUNITEDKINGDOM*WarsawBELARUSBertin*IRELANDGERMANYPOLAND*yvSeUKRAINE美*BudapestAUSTRIAFRANCEHUNGARYLutROMANIASOANBuchareasaidITALYBiackSeaATICANPodgCorsiorBomeMON(FR.TranaARMENIA*MudrdPORTUGA*AnkarardinicSRECSPAINOT)TURKEYLisbon*Cbratarux*waleAlgiersMALTANicasia*CYPRUSMediterraneanSeaIGR
7
FebMarMayJunJulAugSepOctNovDecAnnLatJanAprAmpLonArea数据Amsterda72.92.55.78.24.49.912.514.817.117.114.511.414.652.24.5Westm9.111.711Athens9.715.420.124.527.427.223.819.214.617.818.337.623.5South0.20.14.48.213.81618.31814.4104.21.2Berlin9.118.552.313.2West3.315.617.817.81511.16.74.2Brussels3.36.78.912.84.410.314.450.5West...MarchAprilJuly August SeptemberJanuaryFebruaryMayJuneOctoberNovemberDecember相关系10.640.810.910.990.960.830.640.570.570.970.99January10.69February0.990.980.880.690.620.620.850.930.970.98数矩阵10.950.8March0.960.980.720.720.780.910.960.970.96April0.830.880.9510.940.890.860.90.970.960.920.85May0.640.690.80.9410.970.940.940.940.880.790.680.570.620.720.890.9710.980.960.930.830.740.61JuneJuly0.570.620.720.860.9410.990.930.840.740.620.980.640.690.780.90.940.960.9910.960.890.790.68August0.810.850.910.970.940.930.9610.970.920.84September0.9310.910.930.960.960.880.830.840.890.970.980.93October0.970.970.970.920.790.740.740.790.920.9810.98November10.990.980.960.850.680.610.620.680.840.930.98DecemberJanuanFebruary月份之间特别是相邻的月份之间有较MarchApril老强的相关性，因此我们应用PCA方法May福June试图发现是否能够用少数几个主成分July福（类似于季度）刻画气温特征。0.2AuuSt福0.4September二October-0.6November0.8December8
8 月份之间特别是相邻的月份之间有较 强的相关性，因此我们应用PCA方法 试图发现是否能够用少数几个主成分 （类似于季度）刻画气温特征。 Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Ann Amp Lat Lon Area Amsterda m 2.9 2.5 5.7 8.2 12.5 14.8 17.1 17.1 14.5 11.4 7 4.4 9.9 14.6 52.2 4.5 West Athens 9.1 9.7 11.7 15.4 20.1 24.5 27.4 27.2 23.8 19.2 14.6 11 17.8 18.3 37.6 23.5 South Berlin -0.2 0.1 4.4 8.2 13.8 16 18.3 18 14.4 10 4.2 1.2 9.1 18.5 52.3 13.2 West Brussels 3.3 3.3 6.7 8.9 12.8 15.6 17.8 17.8 15 11.1 6.7 4.4 10.3 14.4 50.5 4.2 West ⋯ January February March April May June July August September October November December January 1 0.99 0.96 0.83 0.64 0.57 0.57 0.64 0.81 0.91 0.97 0.99 February 0.99 1 0.98 0.88 0.69 0.62 0.62 0.69 0.85 0.93 0.97 0.98 March 0.96 0.98 1 0.95 0.8 0.72 0.72 0.78 0.91 0.96 0.97 0.96 April 0.83 0.88 0.95 1 0.94 0.89 0.86 0.9 0.97 0.96 0.92 0.85 May 0.64 0.69 0.8 0.94 1 0.97 0.94 0.94 0.94 0.88 0.79 0.68 June 0.57 0.62 0.72 0.89 0.97 1 0.98 0.96 0.93 0.83 0.74 0.61 July 0.57 0.62 0.72 0.86 0.94 0.98 1 0.99 0.93 0.84 0.74 0.62 August 0.64 0.69 0.78 0.9 0.94 0.96 0.99 1 0.96 0.89 0.79 0.68 September 0.81 0.85 0.91 0.97 0.94 0.93 0.93 0.96 1 0.97 0.92 0.84 October 0.91 0.93 0.96 0.96 0.88 0.83 0.84 0.89 0.97 1 0.98 0.93 November 0.97 0.97 0.97 0.92 0.79 0.74 0.74 0.79 0.92 0.98 1 0.98 December 0.99 0.98 0.96 0.85 0.68 0.61 0.62 0.68 0.84 0.93 0.98 1 数据 相关系 数矩阵
>temperature=read.csv("http://staff.ustc.edu.cn/~ynyang/vector/data/temperature.csv,head=T,row.name=1)>mypca=prcomp（temperaturel1:12l，scale=T）#scale=T：使用相关系数矩阵PCA>summary(mypca)贡献率Standarddeviation是Pc的标PC1PC2PC3PC4PC5准差，即、凡3.2291.1710.347 0.2060.151StandarddeviationProportionofVariance0.8690.1140.0100.0040.00221=3.23,2=1.17,总方差等于相关系数矩阵的trace=12CumulativeProp0rtion0.8690.9830.9930.9960.998前两个PC的方差贡献率：1+/23.232 +1.1720.983二>v=mypcasrotation[,1:2]#载荷矩阵v的前两列1212载荷及解释载荷PC1=-0.27Jan-0.28Feb-0.3Mar-0.31Apr-PC1PC2Pc2=-0.39Jan-0.34Feb-0.21Mar+0.07Apr+..Jan0.27-0.39Feb-0.28-0.34Mar0.3-0.21前两个载荷向量的特点反映PC的意义：-0.31Apr0.07口第1列：PC1的载荷，符号相同取值近似相同，说明PC1代May0.280.34表平均气温。0.260.4JunJul0.270.37口第2列：PC2的载荷，1-3，10-12月（秋冬）的载荷与其0.290.3Aug它月份（春夏）符号相反，所以PC2可能代表温度变化。Sep0.11-0.31口Apr，Oct的PC2载荷接近于0，这是两个气温最舒适的Oct-0.31-0.06月份，它们不出现在PC2中是合理的Nov-0.3-0.21Dec-0.28-0.35下面验证上述解释是否合理。9
9 载荷 PC1 PC2 Jan -0.27 -0.39 Feb -0.28 -0.34 Mar -0.3 -0.21 Apr -0.31 0.07 May -0.28 0.34 Jun -0.26 0.4 Jul -0.27 0.37 Aug -0.29 0.3 Sep -0.31 0.11 Oct -0.31 -0.06 Nov -0.3 -0.21 Dec -0.28 -0.35 前两个载荷向量的特点反映PC的意义： 第1列: PC1的载荷,符号相同取值近似相同, 说明PC1代 表平均气温。 第2列: PC2的载荷，1-3，10-12月（秋冬）的载荷与其 它月份（春夏）符号相反，所以PC2可能代表温度变化。 Apr，Oct的PC2载荷接近于0， 这是两个气温最舒适的 月份，它们不出现在PC2中是合理的. > temperature=read.csv("http://staff.ustc.edu.cn/~ynyang/vector/data/temperature.csv", head=T, row.name=1) > mypca=prcomp( temperature[,1:12] , scale=T ) #scale=T：使用相关系数矩阵 >summary(mypca) PC1 PC2 PC3 PC4 PC5 Standard deviation 3.229 1.171 0.347 0.206 0.151 Proportion of Variance 0.869 0.114 0.010 0.004 0.002 Cumulative Proportion 0.869 0.983 0.993 0.996 0.998 𝜆1 = 3.23, 𝜆2 = 1.17, 总方差 等于相关系数矩阵的trace = 12, 前两个PC的方差贡献率： 𝜆1 + 𝜆2 12 = 3.232 + 1.172 12 = 0.983 Standard deviation是PC的标 准差，即 𝜆𝑖 PC1= −0.27𝐽𝑎𝑛 − 0.28 𝐹𝑒𝑏 − 0.3𝑀𝑎𝑟 − 0.31𝐴𝑝𝑟 − ⋯ PC2 = −0.39𝐽𝑎𝑛 − 0.34 𝐹𝑒𝑏 − 0.21𝑀𝑎𝑟 + 0.07𝐴𝑝𝑟 + ⋯ PCA 贡献率 载荷及 解释 下面验证上述解释是否合理。 > v=mypca$rotation[,1:2] #载荷矩阵V的前两列
数据中还提供了Annual和Amplitude，下面考察关于PC1，2的解验证释是否与这两个变量一致。输出结果中的mypcaSx是所有样本的所有主成分（PC或PCscore），即是35×12矩阵Y=XV。plot(temperaturel"Annual"l,mypcasxl,1])plot(temperaturel,"Amplitude"l,mypcasxl2])nPC1与年度平均气温高度相关（相关系数0.998），PC2与温度变差高度相关（相关系数0.944），这验证了上一页对PC1和PC2的解释是合理的注意：一般数据不提供可用于验证主成分含义的信息10
10 数据中还提供了Annual和Amplitude，下面考察关于PC1，2的解 释是否与这两个变量一致。 输出结果中的mypca$x是所有样本 的所有主成分（PC或PC score），即是35 × 12矩阵𝑌 = 𝑋𝑉。 plot(temperature[, "Annual"], mypca$x[,1]) plot(temperature[, "Amplitude"], mypca$x[,2]) PC1与年度平均气温高度相关(相关系数0.998),PC2与温度变差高度 相关(相关系数0.944),这验证了上一页对PC1和PC2的解释是合理的。 验证 注意：一般数据不提供可用于验证主成分含义的信息