⑴多元回归模型建立
当预测变量也即自变量不止一个时为多元线性回归（multivariable linearregression，MLR），多项式回归可以看成特殊情况下的多元线性回归。现在我们以微生物群落数据为例，探究α多样性指数与环境因子（Salinity、pH、TN、TP，在3.3.2.4VPA分析中这几个变量对微生物群落的解释量较高）之间的关系，如下所示：
代码语言：javascript复制#读取物种和环境因子信息
data=read.csv("otu_table.csv", header=TRUE, row.names=1)
otu=t(data)
envir=read.table("environment.txt", header=TRUE)
rownames(envir)=envir[,1]
env=envir[,-1]
#计算alpha多样性
library(vegan)
shannon=diversity(otu, index="shannon")
biodata=data.frame(shannon, env[,c("Salinity", "pH", "TN", "TP")])
fit=lm(shannon~Salinity pH TN TP, data=biodata)
summary(fit)回归结果如下所示：
整体回归模型检验结果显著，方差解释率60%，但是五个系数中只有3个是显著不为0的。在多元回归中，随着解释变量的增加，无论这些解释变量是否与响应变量有关，R2一般都会增加，这主要是由于随机相关的存在。因此更严谨的来说我们需要根据自由度校正R2，最常用的校正方法如下所示：
上式被称为Ezekiel公式。上面多元回归的结果中已经给出了校正后的R2（51%），我们也可以使用vegan包中的RsquareAdj()函数来校正类多元回归模型（MLR、RDA等）中的R2，如下所示：
代码语言：javascript复制library(vegan)
RsquareAdj(fit)
在上面的多元回归分析中，并没有考虑交互项，但是交互项的解释模型往往使得研究更加有趣，交互影响说明两个解释变量对响应变量的影响是非独立的，例如两种重金属浓度升高时造成的毒性大于单独存在时的毒性，如下所示：
代码语言：javascript复制fit2=lm(shannon~Salinity TN TP:pH, data=biodata)
summary(fit2)
可以看到，使用TP与pH的交互项作为一个预测变量后，系数的显著性整体增强，而在第一次回归中，TP与多样性系数不显著，交互项可以理解为TP对多样性的影响要依赖于pH的水平，也即不同pH下TP对微生物群落的影响不同。复杂的多重多元线性回归可以使用RDA分析来实现。
⑵回归诊断
我们可以使用一元回归诊断方法进行简单的诊断，结果如下：
代码语言：javascript复制par(mfrow=c(2,2))
plot(fit)
在R中car包提供了更详细的回归模型诊断函数，接下来我们对多元回归模型进行详细的评价。
①正态性
可以通过检验残差是否满足t分布来检验正态性，如下所示：
代码语言：javascript复制library(car)
par(mfrow=c(1,1))
qqPlot(fit, labels=names(shannon), id.method="identify", simulate=TRUE)其中simulate=TRUE则95%的置信区间会使用参数自助法（parametric bootstrap）生成，自助法就是假设样本是总体，然后采取有放回抽样方法来确定参数的分布。如下图所示，没有观察到超出置信区间的离群点，也即数据正态性良好：
②残差独立性
接下来检验残差是否相关，可以使用durbinWatsonTest()函数进行Durbin-Waston检验，如下所示：
代码语言：javascript复制durbinWatsonTest(fit, simulate=TRUE, reps=999)
其中参数reps设置了自助抽样的次数，结果p值刚好大于0.05，可以拒绝零假设也即残差相关，说明残差是独立的。不过这个p值很小，仍是有不独立可能的，可以想象，多样性指数越低那么误差范围越小，预测值与观察值越接近，因此残差可能存在不独立性。
③线性
因变量与自变量是否具有线性关系可以通过成分残差图来检验，方法如下：
代码语言：javascript复制crPlots(fit)如下图所示，成分残差图以每一个预测变量作为横坐标，以整体模型的残差加该预测变量和其系数的乘积（也即拟合值中该变量承担的部分）作为纵坐标，如果所有图像均为线性，说明线性关系良好；如果某一变量成分残差图为非线性，说明该变量需要添加多项式项。
④同方差性
可以使用ncvTest()函数检验方差恒定性，如下所示：
代码语言：javascript复制ncvTest(fit)
改检验零假设是误差恒定，p值大于0.05同方差性检验通过。此外，spreadLevelPlot()函数绘制残差与拟合值的关系图，并给出数据转换建议：
代码语言：javascript复制spreadLevelPlot(fit)
如下图所示，随着拟合值的变化先增大后减小（根据实际情况，这很可能是由于变量不均匀性造成的），检验结果给出的建议是对响应变量数据进行2.59次幂次变换（即power transformation）。
⑤多重共线性
在使用多个解释变量进行回归建模时，有时整个模型的显著性非常好，然而回归系数的检验却不显著，这时候很可能出现了多重共线性问题，也即解释变量之间存在较强的相关性。由于回归系数测量的实际上是固定其他解释变量时该解释变量对因变量的影响，共线性会导致回归参数的置信区间过大，是单个系数解释起来很困难。
在生态分析中，环境因子之间很可能会存在共线性问题，这对RDA、CCA、CAP等基于多元回归的模型来说非常重要，因为这些方法使用到了回归系数作为衡量解释变量影响的指标，而VPA分析若要检验每部分方差的显著性也需要消除共线性，LDA分析同样如此。在3.3.2.1RDA分析中我们使用了统计量VIF(variance inflation factor，方差膨胀因子）进行检测，VIF实际上衡量的是回归参数的置信区间能膨胀为与模型无关的解释变量的程度，一般VIF>4则认为存在多重共线性问题，检验方法如下：
代码语言：javascript复制vif(fit)
从结果可以看出，共线性问题并不严重。
⑥筛选特殊点
响应变量中模型预测效果不佳的点称之为离群点，预测变量中异常的预测变量值为高杠杆值点，对模型参数影响过大的点称之为强影响点，也即移除这一观测点模型会发生巨大的改变。对于一个模型来说，我们自然希望每个点影响是一样的，一般来说强影响点既是离群点又是高杠杆值点。在前面的诊断中，已经初步涉及了离群点和高杠杆点的检测，下面提供一些其他的检测方法：
代码语言：javascript复制#检测离群点（校正后p值小于0.05为离群点）
outlierTest(fit)
#可视化方法检测特殊点
influencePlot(fit, id.method="identify", ylim=c(-3,3))
上图横坐标为帽子统计量，越大则杠杆值越大，纵轴为残差，其绝对值越大则越可能为离群点，圆圈越大则该点影响力越大。
