本内容为【科研私家菜】R语言机器学习与临床预测模型系列课程

你想要的R语言学习资料都在这里， 快来收藏关注【科研私家菜】

我们要建立的模型形式是Y = B0 + B1x1 + … + Bnxn + e，预测变量（特征）可以有1 ~ n个。

01 多变量线性回归

检查相关性方面的统计量，并绘出散点图矩阵。
相关系数又称Pearson’s r，可以用来测量两个变量之间线性相关性的强度和方向。相关系数是一个数值，取值在-1和1之间，-1表示完全负相关，+1表示完全正相关。要计算相关系数，
需要使用两个变量的协方差除以它们的标准差的乘积。如果取相关系数的平方，就可以得到R方。

data(water)
dim(water)
str(water)
head(water)
socal.water <- water[ ,-1]
head(socal.water)
library(corrplot)
water.cor <- cor(socal.water)
water.cor
corrplot(water.cor, method = "ellipse")

pairs(~ ., data = socal.water)

效果如下：

响应变量与那些OP开头的特征高度正相关，与OPBPC的相关系数是0.8857，与OPRC的相关系数是0.9196，与OPSLAKE的相关系数是0.9384。还可以看出，AP开头的特征彼此之间高度相关，OP开头的也是如此。我们会遇到多重共线性的问题。

02 特征选择

前向逐步选择

从一个零特征模型开始，然后每次添加一个特征，直到所有特征添加完毕。在这个过程中，被添加的选定特征建立的模型具有最小的RSS。所以理论上，第一个选定的特征应
该能最好地解释响应变量，依此类推。
添加一个特征一定会使RSS减少，使R方增加，但不一定能提高模型的拟合度和可解释性。

后向逐步回归

从一个包含所有特征的模型开始，每次删除一个起最小作用的特征。现在也有一种混合方法，这种算法先通过前向逐步回归添加特征，然后检查是否有特征不再对提高模型拟合度起作用，如果有则删除。每次建模之后，分析者都可以检查模型输出，并使用各种统计量选择能提供最佳拟合的特征。

逐步回归技术会遇到非常严重的问题。对于一个数据集，你先用前向逐步回归，然后再用后向逐步回归，可能会得到两个完全矛盾的模型。最重要的一点是，逐步回归会使回归系数发生偏离.

最优子集回归

最优子集回归是逐步回归的一个可接受的替代方案。在最优子集回归中，算法使用所有可能的特征组合来拟合模型，所以，如果有3个特征，将生成7个模型。然后和逐步回
归一样，分析者需要应用自己的判断和统计分析来选择最优模型，模型选择就是后面工作的关键。
如果数据集有多个特征，工作量就会非常大。当特征数多于观测数时（p大于n），这个方法的效果就不会好。
可以使用lm()函数。模型形式为fit = lm (y~x1 + x2 +x3 + … + xn)。

library(leaps)
fit <- lm(BSAAM ~ ., data = socal.water)
summary(fit)

与单变量回归一样，我们要检查F统计量的p值，以检验是否至少有一个非零系数。确实，p
值是高度显著的。还可以看到，OPRC和OPSLAKE这两个参数具有显著的p值。有趣的是，OPBPC并不显著，尽管它与响应变量高度相关。简言之，当我们控制其他OP开头的特征时，OPBPC无法
对预测方差提供任何有意义的解释。这就是说，模型中存在OPRC和OPSLAKE时，特征OPBPC从统计学角度来看没有任何作用。

03 多变量线性回归模型构建与模型评价

使用最优子集法。我们使用leaps包中的regsubsets()函数建立一个sub.fit对象。summary()就生成了best.summary对象，帮助我们更深入地研究模型。对于R中的所有对象，都可以使用names()函数列出输出结果。

sub.fit <- regsubsets(BSAAM ~ ., data = socal.water)
best.summary <- summary(sub.fit)
names(best.summary)
which.min(best.summary$rss)
# 有6个特征的模型具有最小的RSS

which.min()和which.max()分别给出具有某
个输出的最小值和最大值的模型.
增加特征必然会减少RSS！而且必然会增加R方。
4种用于特征选择的统计方法：赤池信息量准则、马洛斯的Cp、贝叶斯信息量准则和修正R方。前三种方法的目标是追求统计量的值最小化，修正R方的目标是追求统计量的值最大化。这些统计方法的目的是建立一个尽可能简约的模型。

在线性模型中，AIC和Cp成正比，所以我们只需关注Cp。BIC与Cp相比，更倾向于选择变量较少的模型，所以我们要对二者进行比较。

par(mfrow = c(1,2))
plot(best.summary$cp, xlab = "number of features", ylab = "cp")
plot(sub.fit, scale = "Cp")

效果如下：

在左侧的图中，可以看出有3个特征的模型具有最小的Cp值。在右侧的图中，显示了能给出
最小Cp的特征组合。这张图应该这么看：先在Y轴的最高点找到最小的Cp值，此处是1.2；然后
向右在X轴上找到与之对应的色块。通过这张图，我们可以看到这个具有最小Cp值的模型中的3个特征是APSLAKE、OPRC和OPSLAKE。通过which.min()和which.max()函数，我们可以进行Cp与BIC和修正R方的比较。

which.min(best.summary$bic)
which.max(best.summary$adjr2)

可以看出，在本例中，BIC、修正R方和Cp选择的最优模型是一致的。现在，与单变量线性
回归一样，我们需要检查模型并进行假设检验。正像之前做的那样，我们需要生成一个线性模型对象并检查统计图。
用这三个变量构建模型：

best.fit <- lm(BSAAM ~ APSLAKE + OPRC + OPSLAKE, data = socal.water)
summary(best.fit)
par(mfrow = c(2,2))
plot(best.fit)

关注R小盐，关注科研私家菜(VX_GZH: SciPrivate)，有问题请联系R小盐。让我们一起来学习 R语言机器学习与临床预测模型