主成分回归分析实战教程

本文介绍主成分回归分析(Principal Components Regression)，并通过示例展示其实现过程。

给定p个预测变量和响应变量，多元线性回归使用如最小二乘法获得最小误差平方和(RSS)：

RSS = Σ(yi–y^i)2{Σ(y_i – ŷ_i)^2}Σ(yi–y^i)2

-Σ: 求和符号
-yi{y_i}yi: 第i个观测的实际响应值
-y^i{ŷ_i}y^i: 基于多重线性回归模型获得预测值

然而，当预测变量高度相关时，会产生多重共线问题，导致模型系数估计不可靠、高方差。避免该问题的一个方法是使用主成分回归分析 ———— 从p个预测变量中发现M个线性组合（主成分），然后把主成分作为预测变量使用最小二乘法拟合线性回归模型。下面通过示例带你实现主成分回归分析。

加载必要工具包

最简单方式执行主成分回归分析是使用pls包：

# 安装包
install.packages("pls")# 加载包
library(pls)

拟合模型

为了简单方便，我们使用内置的mtcars数据集，包括不同品牌汽车的数据：

head(mtcars)#                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
# Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
# Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
# Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
# Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
# Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
# Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

本例主成分回归分析使用hp作为响应变量，下面变量作为预测变量：

mpg (Miles/(US) gallon)
disp (Displacement)
drat (Rear axle ratio)
wt (Weight)
qsec (1/4 mile time)

下面代码利用上述数据拟合模型，其中两个参数解释如下：

scale = TRUE 每个预测变量都被标准为均值为0，标准差为1，这避免了模型中预测变量使用不同度量单位而产生的影响。
validation = “CV”: 使用K折交叉验证评估模型表现，默认K为10。也可以使用LOOCV参数。


# 是的示例可重现
set.seed(1)#fit PCR model
model <- pcr(hp~mpg+disp+drat+wt+qsec, data=mtcars, scale=TRUE, validation="CV")

选择主成分数

我们依据获得模型，现在需要决定保留主成分数量。下面代码是通过k次交叉验证计算的检验均方根误差(RMSE):

summary(model)# Data:    X dimension: 25 5
#   Y dimension: 25 1
# Fit method: svdpc
# Number of components considered: 5
#
# VALIDATION: RMSEP
# Cross-validated using 10 random segments.
#        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps
# CV           59.98    26.51    26.53    25.91    27.35    30.49
# adjCV        59.98    26.44    26.25    25.59    26.91    29.82
#
# TRAINING: % variance explained
#     1 comps  2 comps  3 comps  4 comps  5 comps
# X     73.45    89.53    95.68    99.04   100.00
# hp    81.11    85.56    87.53    87.65    88.25

我们解释输出中的两个表格：

VALIDATION: RMSEP

这个表告诉我们k折交叉检验的检验RMSE，解雇哦如下：

模型中仅用截距项, 检验RMSE为 69.66.
如果增加第一主成分，检验RMSE减少到 44.56.
如果增加第二主成分，检验RMSE减少到 35.64.

我们看到继续增加主成分会导致RMSE增加，似乎表示仅使用两个主成分为最优模型。

TRAINING: % variance explained

这个表结果表示响应变量中能主成分被解释方差的百分比：

通过只使用第一个主成分，响应变量中被解释方差为69.83%。
通过加入第二个主成分，响应变量中被解释方差为89.35%。

通过使用更多的主成分总是可以解释更多的方差，但我们看到添加两个以上的主成分被解释方差比例实际上增加不明显。使用validationplot()函数能够以可视化方式展示RMSE，从而更直观决定选择主成分数量。

# 2，3 两个图在第一行，1图在第二行占两列（共2行2列）
layout(matrix(c(2,3,1,1),2,2, byrow = TRUE))validationplot(model, val.type="R2")
validationplot(model, val.type="MSEP")
validationplot(model)

在上图中我们可以看到，添加两个主成分时模型的拟合度在提高，但当更多主成分加入时却更差。因此最优模型只包括前两个主成分。

使用最终模型进行预测

我们使用两个主成分模型测试新的观测数据。下面代码把原始数据分为训练集和测试集，使用PCR模型在测试集上预测：

# 定义训练集和测试集
train <- mtcars[1:25, c("hp", "mpg", "disp", "drat", "wt", "qsec")]
y_test <- mtcars[26:nrow(mtcars), c("hp")]
test <- mtcars[26:nrow(mtcars), c("mpg", "disp", "drat", "wt", "qsec")]# 在测试集上进行检验
model <- pcr(hp~mpg+disp+drat+wt+qsec, data=train, scale=TRUE, validation="CV")
pcr_pred <- predict(model, test, ncomp=2)# 计算RMSE
sqrt(mean((pcr_pred - y_test)^2))# [1] 56.86549

我们看到测试RMSE为56.86549，这是测试集中变量hp的预测值与观察值之间的平均偏差。