R语言进阶之主成分分析

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

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

1. 常规主成分分析

在这里，我还将以鸢尾花数据集（iris）为例介绍如何在R中进行主成分分析：

# 主成分分析
# 输入原始数据并提取相关主成分
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) # 输出各个主成分的解释方差

loadings(fit) # 输出载荷

plot(fit,type="lines") # 碎石图

biplot(fit) # 对前两个主成分绘图

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

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

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

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

# 极大方差旋转法
# 保留前两个主成分
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的结果为例介绍本期内容最重要的部分---主成分回归：

# 提取主成分
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则是负值，这和主成分结果一致。

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