专栏首页生物信息云R语言数据分析与挖掘(第六章):主成分分析(2)——案例讲解

R语言数据分析与挖掘(第六章):主成分分析(2)——案例讲解

这一讲通过一个案例讲解主成分分析。

在R中用于完成主成分分析的函数是princomp(),该函数有2种调用方式:

1.公式形式

基本语法为:

princomp(formula, data = NULL, subset, na.action, ...)

参数介绍:

formula:指定用于主成分的公示对象,类似于回归分析和方差分析中的公式对象;

data:指定用于主成分分析的数据对象,一般为数据框;

subset:指定可选的向量,表示选择的样本子集;

na.action:一个函数,指定缺失数据的处理方法,若为NULL,则使用函数na.omit()删除缺失数据。

2. 矩阵形式

基本语法为:

princomp(x, cor = FALSE, scores = TRUE, covmat = NULL,         
         subset = rep_len(TRUE, nrow(as.matrix(x))), fix_sign = TRUE, ...)

参数介绍:

X:指定用于主成分分析的数据对象,一般为数值矩阵的数据框:

Cor:逻辑值,指定用于主成分分析中采用的矩阵形式(相关矩阵或协方差阵),为TRUE表示用样本的相关矩阵做主成分分析,为TALSE表示用样本的协方差阵做主成分分析,默认值为FALSE;

Scores:逻辑值,指定是否计算各主成分的分量,即是否计算样本的主成分得分,默认值TRUE;

Covmati:指定协方差阵,或者为cov.wt()提供的协方差列表,当数据不用参数x提供时,可由协方差阵提供;

Subset:指定可选向量,表示选择的样本子集。

函数princomp()的返回值为-个列表,包括:

sdev表示各主成分的标准差;

loadings表示载荷矩阵;

center表示各指标的样本均值;

scale表示各指标的样本标准差;

n.obs表示观测样本的个数;

scores表示主成分得分(当参数scores= TRUE时提供)。

此外,也可以利用其他函数来提取主成分分析的结果:

  函数summary()可用于提取主成分的信息;

  函数loadings()可用于提取载荷矩阵;

  函数predict()可用于计算主成分得分;

  函数screeplot()可用于绘制主成分的碎石图;

  函数biplot()可用于绘制数据关于主成分的散点图和原坐标在主成分下的方向。

综合案例:longley数据集的变量降维及回归

下面利用第5章提到的longley数据集进行操作演练,该数据集收集了1947年至1962年16年的宏观经济数据,包含7个变量,分别为物价指数平减后的国民生产总值(GNP.deflator)、 国民生产总值(GNP)、(Amed.Frees)人口总量(Popultin)、年份(Year) 和就业人数(Emplyed)。失业人数(Unemployed)、 军队人数(Amed.Forces),人口总(Population)、年份(Year)和就业人数(Employed)。

> data(longley)
> summary(longley)
  GNP.deflator         GNP          Unemployed     Armed.Forces  
 Min.   : 83.00   Min.   :234.3   Min.   :187.0   Min.   :145.6  
 1st Qu.: 94.53   1st Qu.:317.9   1st Qu.:234.8   1st Qu.:229.8  
 Median :100.60   Median :381.4   Median :314.4   Median :271.8  
 Mean   :101.68   Mean   :387.7   Mean   :319.3   Mean   :260.7  
 3rd Qu.:111.25   3rd Qu.:454.1   3rd Qu.:384.2   3rd Qu.:306.1  
 Max.   :116.90   Max.   :554.9   Max.   :480.6   Max.   :359.4  
   Population         Year         Employed    
 Min.   :107.6   Min.   :1947   Min.   :60.17  
 1st Qu.:111.8   1st Qu.:1951   1st Qu.:62.71  
 Median :116.8   Median :1954   Median :65.50  
 Mean   :117.4   Mean   :1954   Mean   :65.32  
 3rd Qu.:122.3   3rd Qu.:1958   3rd Qu.:68.29  
 Max.   :130.1   Max.   :1962   Max.   :70.55  
> corr<-cor(longley)
> library(corrplot)
corrplot 0.84 loaded
> corrplot(corr)

该数据的相关矩阵如上图,可以很直观的看出,数据存在多重共线性问题,下面计算方差膨胀因子(VIF)来判断个变量的共线性问题的严重性:

> library(car)
> vif(lm(Employed~.,data=longley))
GNP.deflator          GNP   Unemployed Armed.Forces   Population 
   135.53244   1788.51348     33.61889      3.58893    399.15102 
        Year 
   758.98060

一般认为:

当0<VIF<10,不存在多重共线性;

当10<=VIF<100,存在较强的多重共线性;

当VIF>= 100,多重共线性非常严重。

上诉输出结果显示,变量GNP.deflator 、GNP、Population和Year存在严重的多重共线性。 为了解决多重共线性问题,下面对数据集进行主成分分析,去掉响应变量后进行分析:

> (pr1<-princomp(longley[,-7], cor = TRUE))
Call:
princomp(x = longley[, -7], cor = TRUE)

Standard deviations:
    Comp.1     Comp.2     Comp.3     Comp.4     Comp.5     Comp.6 
2.14554820 1.08413122 0.45102702 0.12218125 0.05051797 0.01940897

 6  variables and  16 observations.
> summary(pr1)
Importance of components:
                          Comp.1    Comp.2     Comp.3      Comp.4
Standard deviation     2.1455482 1.0841312 0.45102702 0.122181253
Proportion of Variance 0.7672295 0.1958901 0.03390423 0.002488043
Cumulative Proportion  0.7672295 0.9631196 0.99702383 0.999511871
                             Comp.5       Comp.6
Standard deviation     0.0505179747 1.940897e-02
Proportion of Variance 0.0004253443 6.278469e-05
Cumulative Proportion  0.9999372153 1.000000e+00

输出结果中:Standard deviation表示主成分的标准差,即特征值的开方;Proportion of Variance表示每个主成分的贡献率;Cumulative Proportion表示主成分的累计贡献率,通常情况下,参考该参数的值进行主成分的选择,一般累计贡献率超过85%的主成分,在本案例中,前2个主成分的累计贡献率已经达到96%,即前2个主成分能解释原始变量的96%的信息,故应该选择前2个主成分,剔除后面4个主成分,达到降维的目的。

下面利用函数loadings()输出载荷矩阵:

> loadings(pr1)

Loadings:
             Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
GNP.deflator  0.462         0.149  0.793  0.338  0.135
GNP           0.462         0.278 -0.122 -0.150 -0.818
Unemployed    0.321 -0.596 -0.728               -0.107
Armed.Forces  0.202  0.798 -0.562                     
Population    0.462         0.196 -0.590  0.549  0.312
Year          0.465         0.128        -0.750  0.450

               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
SS loadings     1.000  1.000  1.000  1.000  1.000  1.000
Proportion Var  0.167  0.167  0.167  0.167  0.167  0.167
Cumulative Var  0.167  0.333  0.500  0.667  0.833  1.000

可以得到|成分!原的变量的线性关系:

第1主成分:

F1=0.462*CNPdnator-0.462*GNP-0.321*Unemploycd-0.202*Amed.Forces-0.462*Ppulation-0.465*Year

第2主成分

Fn=-0.596*Unemployed+0.798* Armed.Forces

下面对选样的主成分进行命名,主成分命名的原则应注意以下几点:

由于F1中包含不同变量的信息,在进行命名时需费综合考虑不同变量的信息;需要注意线性关系中系数的符号和绝对值大小。

基于上述原则,第一主成分对应系数的符号相同,均为负号,除变最Unemployed和Amed.Forces以外,其他原始变量的系数绝对值均为0.5左有,故第1主成分反映了经济的萧条程度,经济越萧条,F1的取值越大,各原始变量的取值越小,反之,经济越不萧条,F1的取值越小,各原始变量的取值越大。

第二主成分要与原始变量Unemployed和Amed.Forces有关,且变Unemployed对应的系数为负,Amed.Forces对应的系数为正,系数的绝对值均大于0.5,故第二主成分反映了军备程度,军备程度越高,F2的取值越大,失业人数越小,且军队人数越多。反之,军备程度越低,F2的取值越小,失业人数越多,军队人数越少。

  尽管上述分析对两个主成分进行了命名,但命名结果还是不能很好地解释原版始变量的信息,因为主成分命名本身就具有较人的难度,我们能做到的就是尽量个面地涵盖原始变量的信息。

下面是利用函数predict()输出主成分的得分。

> predict(pr1)
          Comp.1     Comp.2      Comp.3       Comp.4       Comp.5
1947 -3.59294203 -0.7761110  0.31804507 -0.169624216  0.009085721
1948 -3.10924187 -0.8768853  0.66329350  0.130051755  0.063565251
1949 -2.42014988 -1.5905016 -0.50961596 -0.009106066  0.005934717
1950 -2.16257276 -1.3181781 -0.11494311 -0.063265342 -0.063873551
1951 -1.48540767  1.2763227 -0.03004947  0.100660376  0.053969833
1952 -1.04339031  1.9851409 -0.16650384  0.047914580  0.038252420
1953 -0.72546382  1.9732212  0.06933769 -0.073217347  0.022425386
1954  0.03355589  0.6124972 -1.07307182 -0.067293112 -0.002331281
1955  0.10277563  0.7162362 -0.10076710 -0.104426569 -0.102049470
1956  0.46417097  0.5658113  0.30255566  0.018137849 -0.086508832
1957  0.98638447  0.4435320  0.45984009  0.123243732 -0.024470311
1958  1.87669233 -0.8914800 -0.69963411  0.193195372  0.022381351
1959  2.00361238 -0.3992485  0.27468475  0.148640632 -0.037890240
1960  2.43855615 -0.5154723  0.37766452  0.063619308 -0.016767727
1961  3.17896368 -1.0224220 -0.20859150 -0.070338720  0.058280241
1962  3.45445683 -0.1824626  0.43775561 -0.268192230  0.059996493
           Comp.6
1947  0.002663801
1948  0.012373058
1949  0.005226048
1950 -0.014125087
1951 -0.044081763
1952 -0.013169060
1953  0.032810005
1954  0.017241531
1955 -0.019544921
1956  0.014604365
1957  0.028044545
1958  0.008369752
1959 -0.024302116
1960  0.004500429
1961 -0.001372125
1962 -0.009238462
> sort(predict(pr1)[,1],decreasing = T)
       1962        1961        1960        1959        1958        1957 
 3.45445683  3.17896368  2.43855615  2.00361238  1.87669233  0.98638447 
       1956        1955        1954        1953        1952        1951 
 0.46417097  0.10277563  0.03355589 -0.72546382 -1.04339031 -1.48540767 
       1950        1949        1948        1947 
-2.16257276 -2.42014988 -3.10924187 -3.59294203 
> sort(predict(pr1)[,2],decreasing = T)
      1952       1953       1951       1955       1954       1956 
 1.9851409  1.9732212  1.2763227  0.7162362  0.6124972  0.5658113 
      1957       1962       1959       1960       1947       1948 
 0.4435320 -0.1824626 -0.3992485 -0.5154723 -0.7761110 -0.8768853 
      1958       1961       1950       1949 
-0.8914800 -1.0224220 -1.3181781 -1.5905016

从第一主成分看,1947、1948和1949年的得分较高,说明那几年的经济较为萧条,得分较低的是1960、1961和1962年,说明那几年的经济形势较好。

从第二主成分来看,1952、1953和1951年得分较高,说明那几年的军备情况较好,1961、1950和1949年的得分较低,说明那几年的军备情况较差。

下面利用函数是绘制主成分碎石图:

screeplot(pr1,type="lines",main="")

绘图形势选择折线图,碎石图展现的是各主成分没有解释道原始数据的信息情况,折线图越靠近x轴,说明主成分解释原始变量的信息越多。可以看出选择前2个主成分是合理的。

下面利用函数biplot()绘制双坐标图,默认情况下,该函数绘制第一、二主成分样本散点图和原始坐标在第一、二主成分下的方向。

biplot(pr1)

主成分回归 主成分回归可以参考前面讲的回归分析看代码:

pre1<-predict(pr1)
longley$F1<-pre1[,1]
longley$F2<-pre1[,2]
longley$F3<-pre1[,3]
lm.pr1<-lm(Employed~F1+F2+F3,data=longley)
summary(lm.pr1)
beta<-coef(lm.pr1)
A<-loadings(pr1)[,1:3]
x.mean<-pr1$center
x.sd<-pr1$scale
coef<-A%*%beta[2:4]/x.sd
beta0<-beta[1]-x.mean%*%coef
c(beta0,coef)

本文分享自微信公众号 - MedBioInfoCloud(MedBioInfoCloud),作者:DoubleHelix

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-10-18

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • R语言数据分析与挖掘(第六章):主成分分析(1)——主成分分析概论

    在许多领域的研究与应用中,往往需要对反映事物的多个变量进行大量的观测,收集大量数据以便进行分析寻找规律。多变量大样本无疑会为研究和应用提供了丰富的信息,但也在一...

    DoubleHelix
  • 医学生信文献第10期:一篇肿瘤领域入门必读综述——新一代癌症标志物

    在人类肿瘤的多 级发展过程 形成的六 个生物功 能构成了癌症的 特征。这些 特征为分析复 杂的肿瘤性疾病 提供了一个 组织原则。这六个 特征包括:持久 的增殖信...

    DoubleHelix
  • R语言基础绘图教程——第3章:折线图和带状图

    在上一章中我们讲过plot()绘图的基本结构,主要通过type参数来设置绘制图形的类型。

    DoubleHelix
  • 如何清理photoshop cs6 被升级的烦人的adobe creative cloud组件

    安装photoshop cs6(虽然目前已经退出到cc 2015,不过因激活成熟度等,我还是偏向于使用cs6,够用!),默认安装adobe applicatio...

    williamwong
  • vue 控制某个元素的显示与隐藏之v-if属性

    解释: 默认showPrise和showRentPrise的状态是false,所以是隐藏的。  当你在changeStatus通过了某种条件,你就可以控制s...

    acoolgiser
  • ​lcc-render可调自定义渲染框架!附源码仓库

    作者:Nomat 来源:Cocos官方论坛 原文:https://forum.cocos.org/t/topic/99268

    张晓衡
  • 用 Blazor WebAssembly 实现微前端

    我聊下最近我在做的事情,然后分享下在Blazor WebAssembly 微前端的实现细节,这篇文章是我的一些心得,以及一个示例的 Demo 项目,展示了如何使...

    全球技术精选
  • 「译」 用 Blazor WebAssembly 实现微前端

    我聊下最近我在做的事情,然后分享下在Blazor WebAssembly 微前端的实现细节,这篇文章是我的一些心得,以及一个示例的 Demo 项目,展示了如何使...

    全球技术精选
  • Javascript 的addEventListener()及attachEvent()区别分析

    Mozilla中:  addEventListener的使用方式:  target.addEventListener(type, listener, useCa...

    smy
  • Java类加载基本过程

    数组类本身不通过类加载器创建,由java虚拟机直接创建,数组类的元素类型由类加载器加载。

    WindWant

扫码关注云+社区

领取腾讯云代金券