R语言进阶之主成分分析

2020-08-05 16:51:53 浏览数 (1)

‍今天我们将要学习R语言进阶中最重要的统计内容---主成分分析,它在我们的研究中几乎是无处不在,应用最广的就是将主成分放入回归模型进行拟合,用于矫正相关的混杂因素。

主成分分析的基本思想是将多个变量进行线性组合,在保留原数据主要特征的同时减少变量个数,从而达到降维的目的。R语言的内置函数princomp()提供了未经旋转的主成分分析。

1. 常规主成分分析

在这里,我还将以鸢尾花数据集(iris)为例介绍如何在R中进行主成分分析:

代码语言:javascript复制
# 主成分分析
# 输入原始数据并提取相关主成分
mydata <- iris # 将iris命名成mydata以便后续操作
mydata$type[which(mydata$Species== "setosa")]<- 0 # 设置setosa为0
mydata$type[which(mydata$Species== "versicolor")]<- 1 # 设置versicolor为1
mydata$type[which(mydata$Species== "virginica")]<- 2 # 设置virginica为2
attach(mydata) # 固定数据框
fit <- princomp(~ Sepal.Length Sepal.Width Petal.Length Petal.Width, cor=TRUE) # 对花萼和花瓣的相关数据进行主成分分析,cor=TRUE表示从相关系数矩阵提取主成分(实际上是对数据的一种标准化)
summary(fit) # 输出各个主成分的解释方差
代码语言:javascript复制
loadings(fit) # 输出载荷
代码语言:javascript复制
plot(fit,type="lines") # 碎石图
代码语言:javascript复制
biplot(fit) # 对前两个主成分绘图

从上面的第一幅图来看,前两个主成分的累计方差贡献率达到95.8%,并且碎石图的结果也显示前两个主成分所占方差较大,因此我们其实只要用这两个主成分就能很好描述鸢尾花的特征了。

各个主成分的载荷实际上反应的是各原始变量和主成分的关系,从图中结果我们不难看出,主成分1主要反映花萼长度、花瓣的长度和花瓣的宽度这三个原始变量,而主成分2主要反映花萼宽度这个原始变量,因此前两个主成分基本就能完全反映所有的变量特征了。

最后一幅图实际上是按照前两个主成分绘制的散点图,从图中不难看出:利用主成分1可以将“setosa”与其他两类分开,然后利用主成分2将“versicolor“和”virginica“分开。

这里我想和大家介绍一下“psych“包(一个十分强大的统计R包)的主成分函数principal( ),这个函数可以帮助我们提取和旋转主成分:

代码语言:javascript复制
# 极大方差旋转法
# 保留前两个主成分
library(psych)
fit2 <- principal(mydata[,1:4], nfactors=2,rotate="varimax") #nfactors指定主成分个数,rotate指定旋转方法
fit2 # 输出结果

关于rotate参数,我们主要有如下选项 "none", "varimax", "quatimax","promax", "oblimin", "simplimax"和"cluster",有兴趣的朋友可以自行了解,我在这儿就不赘述了。

从上图的结果来看,前两个主成分的累计方差贡献率为96%,并且第一主成分主要表征花萼长度、花瓣的长度和花瓣的宽度这三个原始变量,而第二主成分主要反映花萼宽度这个原始变量,这和之前的分析一致。

2. 主成分回归

接下来,我将以fit2的结果为例介绍本期内容最重要的部分---主成分回归:

代码语言:javascript复制
# 提取主成分
detach(mydata) # 解固定原数据
PC <- as.data.frame(fit2$scores) # 提取各观测点的主成分
PC <- cbind(PC,mydata$type) # 按列合并
colnames(PC) <- c("PC1", "PC2","type") # 重命名变量
fit3 <- lm(type ~ PC1   PC2,data=PC) # 利用前两个主成分对type进行线性回归
summary(fit3)

输出结果显示PC1和PC2都是非常显著的,而PC1的效应量为正值,PC2则是负值,这和主成分结果一致。

关于主成分分析的内容就讲到这里,希望大家能掌握主成分分析的方法以及如何用主成分进行回归分析,咱们下期再见!

0 人点赞