第8章 主成分回归与偏最小二乘

8.3 对例5.5的Hald水泥问题用主成分回归方法建立模型，并与其他方法的结果进行比较。

8.4 对例5.5的Hald水泥问题用偏最小二乘方法建立模型，并与其他方法的结果进行比较。

x1 x2 x3 x4 y
7 26 6 60 78.5
1 29 15 52 74.3
11 56 8 20 104.3
11 31 8 47 87.6
7 52 6 33 95.9
11 55 9 22 109.2
3 71 17 6 102.7
1 31 22 44 72.5
2 54 18 22 93.1
21 47 4 26 115.9
1 40 23 34 83.8
11 66 9 12 113.3
10 68 8 12 109.4

rm(list=ls())

# ----- 例5.5 水泥 逐步回归法 ----
data5.5 <- read.csv('D:/rwork/应用回归/例题数据/data5.5.csv',head=TRUE)
lm5.5 <- lm(y~.,data=data5.5) #y~.表示将y对data5.5中其余的所有变量做回归
lm5.5_step <- step(lm5.5,direction='both') #初始模型包含所有变量
summary(lm5.5_step) #输出依据AIC选出的最优模型的结果
# 由输出结果可见，逐步回归筛选的最优子集为x1,x2,x4，但在显著性水为0.05时x4的回归系数不显著，
#  从上述输出结果可知，由最小的AIC值选出的模型在整体上最优，但是可能会包含不显著的变量。
# 故需要删去不显著的变量x4，得到新的回归结果如下
summary(lm(y~x1+x2,data=data5.5))
# 从上述输出结果可知，逐步回归法得到的回归方程为y^=52.577+1.468x1+0.662x2。



# ---- xt8.3 对例5.5的Hald水泥问题用主成分回归方法建立模型 ----
data5.5 <- read.csv('D:/rwork/应用回归/例题数据/data5.5.csv',head=TRUE)
datas8.3 <- data.frame(scale(data5.5)) #对原始数据进行标准化处理
pr8.3 <- princomp(~x1+x2+x3+x4,data=datas8.3,cor=T)
#对4个自变量做主成分分析，其中cor=T表明是用相关系数矩阵进行主成分分析
summary(pr8.3,loadings=TRUE) #输出主成分分析的主要结果
pr8.3$scores[,1:2] #输出前两个主成分的得分

# 由于前两个主成分的方差累计贡献率已经达到95.294%即可解释大部分变差，
#  故只需保留前两个主成分，此处只输出前两个主成分的得分。

# 现用y对前两个主成分做普通最小二乘回归
pre8.3 <- pr8.3$scores[,1:2] #将前两个主成分的得分保存在变量pre8.3中
datas8.3$z1 <- pre8.3[,1] #将第一主成分的得分添加在数据框datas中，变量名为z1
datas8.3$z2 <- pre8.3[,2] #将第一主成分的得分添加在数据框datas中，变量名为z2
pcr8.3 <- lm(y~z1+z2-1,data=datas8.3) #y对两个主成分建立回归模型
summary(pcr8.3)
# 得到y^*=0.631z1-0.008z2，
# 将z1=0.476x1*+0.564x2*-0.394x3*-0.548x4*，z2=0.509x1*-0.414x2*-0.605x3*+0.451x4*，代入上式，
# 即可得到y^*=0.296x1*+0.359x2*-0.244x3*-0.349x4*

# 两种方法的主要区别在于：普通最小二乘法认为自变量对因变量直接起作用，故要剔除对因变量作用不大
#  的自变量；而主成分回归方程则是寻找影响自变量的主要因子，关注这些因子对因变量的作用。
#  两种方法的选择应结合实际情况出发，从而选取最优的解决方案。



# ----- 例5.5 水泥 普通最小二乘法 ----
summary(lm(y~.,data=data5.5))
# 普通最小二乘得到的回归方程为y^=62.41+1.55x1+0.51x2+0.10x3-0.14x4，
#  从系数上看可以发现明显不合理的地方，x3与y是负相关，但它的系数确实正的。



# ---- xt8.4 对例5.5的Hald水泥问题用偏最小二乘方法建立模型 ----
datas8.4 <- data.frame(scale(data5.5)) #对原始数据进行标准化处理
library(pls)
pls1 <- plsr(y~.,data=datas8.4,validation='LOO',jackknife=TRUE,method='widekernelpls')
# 使用偏最小二乘法建立回归方程，其中validation='LOO'表示使用留一交叉验证计算RMSEP；
#  jackknife=TRUE表示使用jackknife方法估计回归系数方差（为后面的显著性检验做准备）。
#  在不设定主成分个数（ncomp）时，默认使用所有主成分进行分析。
summary(pls1,what='all') #输出回归结果：默认误差均方根RMSEP和变异解释度

# 从使用了所有主成分进行回归的结果可以看出，主成分个数为3个时，模型在经过留一交叉验证法后求得的
#  RMSEP总和最小，且随着成分个数的增加，RMSEP值未出现明显减少，同时3个主成分对各个因变量的累计
#  贡献率均高于98%，因此将回归的主成分个数定为m=3。

pls3 <- plsr(y~.,data=datas8.4,ncomp=3,validation='LOO',jackknife=TRUE)
coef(pls3) #得到方程的回归系数
# 得到标准化后的数据的回归方程为y*=0.51x1*+0.28x2*-0.06x3*-0.42x4*，
#  还原为原始变量为y^=85.50+1.31x1+0.27x2-0.14x3-0.38x4。

# 从系数上看，x1,x2对y起负影响，与相关分析得到的结果一致，因此偏最小二乘回归系数的解释比
#  普通最小二乘更合理，又比逐步回归保留了更多的自变量。



## 将回归方程中的变量还原为原始变量 ----
# 利用zscore(x)函数将数据标准化的核心思想是z=(x-mean(x))./std(x)。
# 得到标准化后的数据的回归方程为y*=0.51x1*+0.28x2*-0.06x3*-0.42x4*，将回归方程还原为原始方程：
# y*=(y-mean(y))/sd(y)=0.51*(x1-mean(x1))/sd(x1)+0.28*(x2-mean(x2))/sd(x2)-0.06*(x3-mean(x3))/sd(x3)-0.42*(x4-mean(x4))/sd(x4)，
# 即y=sd(y)*[0.51*(x1-mean(x1))/sd(x1)+0.28*(x2-mean(x2))/sd(x2)-0.06*(x3-mean(x3))/sd(x3)-0.42*(x4-mean(x4))/sd(x4)]+mean(y),
# 即还原为原始变量为y^=85.50+1.31x1+0.27x2-0.14x3-0.38x4。
attach(data5.5) #将该数据框添加到R的搜索路径，以便于下面直接使用数据框中包含的数组x和y
b1 <- sd(y)*0.51398740/sd(x1)
b1 #原始变量x1的系数为1.314479
b2 <-sd(y)* 0.28126396/sd(x2)
b2 #原始变量x2的系数为0.2719163
b3 <- sd(y)*(-0.05967267)/sd(x3)
b3 #原始变量x3的系数为-0.1401532
b4 <-sd(y)* (-0.42014716)/sd(x4)
b4 #原始变量x4的系数为-0.3776144
b0 <- sd(y)*(0.51398740*(-mean(x1))/sd(x1)+0.28126396*(-mean(x2))/sd(x2)-0.05967267*(-mean(x3))/sd(x3)-0.42014716*(-mean(x4))/sd(x4))+mean(y)
b0 #原始变量y与x的回归方程的截距为85.49915
detach(data5.5) #与attach()相对应，将数据框从搜索路径中移除


参考课本：应用回归分析（R语言版），何晓群编著