Press "Enter" to skip to content

R语言应用统计1 主成分分析

本站内容均来自兴趣收集,如不慎侵害的您的相关权益,请留言告知,我们将尽快删除.谢谢.

R语言应用统计1 主成分分析

 

这个系列就讨论应用基础,争取一条公式都不用写。当原始数据集比较庞大,并且不同变量之间存在一些相关性时,我们希望可以用更少的变量来表示原始数据集,用到的变量越少的同时,能够表示的原始数据集中的信息越多自然就更好。主成分分析就可以实现这样的目标,在主成分分析中用来表示原始数据集中的信息的变量被称为主成分。下面我们用一个例子说明R语言中进行简单的主成分分析的方法。

 

数据

 

使用 HSAUR2 包中的美国城市污染数据,代码如下

 

install.packages("HSAUR2")
data("USairpollution", package = "HSAUR2")

 

这个数据集有七个变量,分别是代表空气污染的 SO2 (空气中的二氧化硫浓度)、代表人类活动的 populmanu 以及代表天气情况的 windtempprecippredays ,这个数据集可以用多元线性回归来分析代表人类活动与天气情况的六个变量对二氧化硫浓度的影响,但现在我们尝试用PCA对六个解释变量进行分析。

 

数据的定性分析

 

首先分析一下这六个解释变量之间的相关性,使用 cor 函数,代码为

 

cor(USairpollution[,-1])

 

代码输出如下

 

temp        manu       popul        wind      precip     predays
temp     1.00000000 -0.19004216 -0.06267813 -0.34973963  0.38625342 -0.43024212
manu    -0.19004216  1.00000000  0.95526935  0.23794683 -0.03241688  0.13182930
popul   -0.06267813  0.95526935  1.00000000  0.21264375 -0.02611873  0.04208319
wind    -0.34973963  0.23794683  0.21264375  1.00000000 -0.01299438  0.16410559
precip   0.38625342 -0.03241688 -0.02611873 -0.01299438  1.00000000  0.49609671
predays -0.43024212  0.13182930  0.04208319  0.16410559  0.49609671  1.00000000

 

可以发现代表人类活动的 populmanu 的相关性系数高达0.9553,代表天气情况的四个变量之间也存在一定的相关性,这说明如果直接用多元线性回归分析这个数据集可能存在多重共线性的问题。

 

接下来我们画出这六个变量两两之间的散点图,代码如下

 

panel.hist <- function(x, ...) {
   usr <- par("usr"); on.exit(par(usr))
   par(usr = c(usr[1:2], 0, 1.5) )
   h <- hist(x, plot = FALSE)
   breaks <- h$breaks; nB <- length(breaks)
   y <- h$counts; y <- y/max(y)
   rect(breaks[-nB], 0, breaks[-1], y, col="grey", ...)
}
USairpollution$negtemp <- USairpollution$temp * (-1) 
USairpollution$temp <- NULL
pdf("Lot0.pdf") #将图片以pdf形式导出到当前工作路径
pairs(USairpollution[,-1], diag.panel = panel.hist,
       pch = ".", cex = 1.5)
dev.off() #结束导出

 

 

从这些散点图中可以看出,至少 populmanu 之间, precipnegtemp 之间存在比较明显的线性相关性。 negtemp 就是取 temp 的相反数,可以这样做但没必要,只是取相反数画出来图看上去线性性要明显一点而已。

 

主成分分析

 

使用 princomp 函数做主成分分析并用 summary 查看结果,代码如下,

 

usair_pca <- princomp(USairpollution[,-1], cor = TRUE)
summary(usair_pca, loadings = TRUE)

 

cor=TRUE 代表在计算中需要用到相关性矩阵,因为在定性分析的时候我们发现六个变量之间确实存在一定的相关性,所以需要用到,如果大家在尝试其他数据集,发现相关性系数矩阵非常接近单位矩阵,就可以用 cor=FALSE ,如果变量中存在常值变量,就只能用 cor=FALSEloadings=TRUE 代表在展示结果时同时展示各个主成分的载荷。代码输出如下

 

Importance of components:
                          Comp.1    Comp.2    Comp.3    Comp.4     Comp.5      Comp.6
Standard deviation     1.4819456 1.2247218 1.1809526 0.8719099 0.33848287 0.185599752
Proportion of Variance 0.3660271 0.2499906 0.2324415 0.1267045 0.01909511 0.005741211
Cumulative Proportion  0.3660271 0.6160177 0.8484592 0.9751637 0.99425879 1.000000000
Loadings:
        Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
temp     0.330  0.128  0.672  0.306  0.558  0.136
manu    -0.612  0.168  0.273 -0.137 -0.102  0.703
popul   -0.578  0.222  0.350               -0.695
wind    -0.354 -0.131 -0.297  0.869  0.113       
precip         -0.623  0.505  0.171 -0.568       
predays -0.238 -0.708        -0.311  0.580

 

第二行Proportion of Variance代表了每一个主成分包含的原始数据集的信息,第三行Cumulative Proportion代表前几个主成分累积包含的原始数据集的信息,比如假设选择主成分的标准是被选择的主成分包含的信息超过原始数据集的95%,则我们应该选前四个主成分,这样就实现了用四个变量代表原始的六个变量的简化目标,并且根据Loadings,我们可以用原始的六个变量的线性组合来表示这四个主成分,以第一主成分为例:

 

C o m p 1 = 0.330 t e m p − 0.612 m a n u − 0.578 p o p u l − 0.354 w i n d − 0.238 p r e d a y s Comp1 = 0.330temp-0.612manu-0.578popul \\-0.354wind-0.238 predays C o m p 1 = 0 . 3 3 0 t e m p − 0 . 6 1 2 m a n u − 0 . 5 7 8 p o p u l 0 . 3 5 4 w i n d − 0 . 2 3 8 p r e d a y s

 

接下来我们分析这个主成分分析的score,score是每一个主成分在单个数据点上的得分,得分绝对值越高说明这个主成分在这个数据点上的表现越差,我们用 MVA 包中的 bvbox 函数画出score的散点图,

 

install.packages('MVA')
library(MVA)

 

代码如下,

 

pdf("Lot1.pdf")
pairs(usair_pca$scores[,1:4], ylim = c(-6, 4), xlim = c(-6, 4), 
      panel = function(x,y, ...) {
      text(x, y, abbreviate(row.names(USairpollution)),
           cex = 0.6)
      bvbox(cbind(x,y), add = TRUE)
      }) #[,1:4]代表分析前四个主成分
dev.off()

 

 

bvbox 的作用是画二维的箱线图,在外圈虚线以外的点被认为是异常值,可以发现大部分数据点的信息都是可以被这四个主成分所代表的。

 

分析主成分与二氧化硫浓度的关系

 

现在我们用这些结果讨论主成分与二氧化硫浓度之间的关系,此时被解释变量是二氧化硫浓度,解释变量是主成分的score,代码如下

 

usair_reg <- lm(SO2 ~ usair_pca$scores, data = USairpollution)
summary(usair_reg)

 

输出如下

 

Call:
lm(formula = SO2 ~ usair_pca$scores, data = USairpollution)
Residuals:
    Min      1Q  Median      3Q     Max 
-23.004  -8.542  -0.991   5.758  48.758 
Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)              30.049      2.286  13.146 6.91e-15 ***
usair_pca$scoresComp.1   -9.942      1.542  -6.446 2.28e-07 ***
usair_pca$scoresComp.2   -2.240      1.866  -1.200  0.23845    
usair_pca$scoresComp.3    0.375      1.935   0.194  0.84752    
usair_pca$scoresComp.4   -8.549      2.622  -3.261  0.00253 ** 
usair_pca$scoresComp.5  -15.176      6.753  -2.247  0.03122 *  
usair_pca$scoresComp.6   39.271     12.316   3.189  0.00306 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 14.64 on 34 degrees of freedom
Multiple R-squared:  0.6695,Adjusted R-squared:  0.6112 
F-statistic: 11.48 on 6 and 34 DF,  p-value: 5.419e-07

 

可以发现这个多元线性回归模型是显着的,第1个主成分是显着的,在更弱一点的显着性要求下,第4、5、6个主成分也是显着的。R方说明用这六个主成分可以解释二氧化硫浓度变化信息的61.12%。

Be First to Comment

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注