R做多元线性回归全攻略
参考链接：
R中的线性回归函数比较简单，就是lm()，比较复杂的是对线性模型的诊断和调整。这里结合Statistical Learning和杜克大学的Data Analysis and Statistical Inference的章节以及《R语言实战》的OLS(Ordinary Least Square)回归模型章节来总结一下，诊断多元线性回归模型的操作分析步骤。
1、选择预测变量
因变量比较容易确定，多元回归模型中难在自变量的选择。自变量选择主要可分为向前选择（逐次加使RSS最小的自变量），向后选择（逐次扔掉p值最大的变量）。个人倾向于向后选择法，一来p值比较直观，模型返回结果直接给出了各变量的p值，却没有直接给出RSS；二来当自变量比较多时，一个个加比较麻烦。
Call:
lm(formula = Sales ~ . + Income:Advertising + Age:Price, data = Carseats)
Residuals:
Min 1Q Median 3Q Max
-2.9208 -0.7503 0.0177 0.6754 3.3413
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.5755654 1.0087470 6.519 2.22e-10 ***
CompPrice 0.0929371 0.0041183 22.567 < 2e-16 ***
Income 0.0108940 0.0026044 4.183 3.57e-05 ***
Advertising 0.0702462 0.0226091 3.107 0.002030 **
Population 0.0001592 0.0003679 0.433 0.665330
Price -0.1008064 0.0074399 -13.549 < 2e-16 ***
ShelveLocGood 4.8486762 0.1528378 31.724 < 2e-16 ***
ShelveLocMedium 1.9532620 0.1257682 15.531 < 2e-16 ***
Age -0.0579466 0.0159506 -3.633 0.000318 ***
Education -0.0208525 0.0196131 -1.063 0.288361
UrbanYes 0.1401597 0.1124019 1.247 0.213171
USYes -0.1575571 0.1489234 -1.058 0.290729
Income:Advertising 0.0007510 0.0002784 2.698 0.007290 **
Price:Age 0.0001068 0.0001333 0.801 0.423812
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.011 on 386 degrees of freedom
Multiple R-squared: 0.8761, Adjusted R-squared: 0.8719
F-statistic: 210 on 13 and 386 DF, p-value: < 2.2e-16构建一个回归模型后，先看F统计量的p值，这是对整个模型的假设检验，原假设是各系数都为0，如果连这个p值都不显著，无法证明至少有一个自变量对因变量有显著性影响，这个模型便不成立。然后看Adjusted R2，每调整一次模型，应该使它变大；Adjusted R2越大说明模型中相关的自变量对因变量可解释的变异比例越大，模型的预测性就越好。
构建了线性模型后，如果是一元线性回归，可以画模型图初步判断一下线性关系（多元回归模型不好可视化）：
par(mfrow=c(1,1))
plot(medv~lstat,Boston)
fit1=lm(medv~lstat,data=Boston)
abline(fit1,col="red")
2、模型诊断
确定了回归模型的自变量并初步得到一个线性回归模型，并不是直接可以拿来用的，还要进行验证和诊断。诊断之前，先回顾多元线性回归模型的假设前提（by Data Analysis and Statistical Inference）：
（数值型）自变量要与因变量有线性关系；
残差基本呈正态分布；
残差方差基本不变（同方差性）；
残差（样本）间相关独立。
一个好的多元线性回归模型应当尽量满足这4点假设前提。
用lm()构造一个线性模型fit后，plot(fit)即可返回4张图（可以par(mfrow=c(2,2))一次展现），这4张图可作初步检验：
左上图用来检验假设1，如果散点看不出什么规律，则表示线性关系良好，若有明显关系，则说明非线性关系明显。右上图用来检验假设2，若散点大致都集中在QQ图中的直线上，则说明残差正态性良好。左下图用来检验假设3，若点在曲线周围随机分布，则可认为假设3成立；若散点呈明显规律，比如方差随均值而增大，则越往右的散点上下间距会越大，方差差异就越明显。假设4的独立性无法通过这几张图来检验，只能通过数据本身的来源的意义去判断。
右下图是用来检验异常值。异常值与三个概念有关：
离群点：y远离散点主体区域的点
杠杆点：x远离散点主体区域的点，一般不影响回归直线的斜率
强影响点：影响回归直线的斜率，一般是高杠杆点。
对于多元线性回归，高杠杆点不一定就是极端点，有可能是各个变量的取值都正常，但仍然偏离散点主体。
对于异常值，可以谨慎地删除，看新的模型是否效果更好。
《R语言实战》里推荐了更好的诊断方法，总结如下。
1、多元线性回归假设验证：
gvlma包的gvlma()函数可对拟合模型的假设作综合验证，并对峰度、偏度进行验证。
states <- as.data.frame(state.x77[, c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income +
Frost, data = states)
library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)
Call:
gvlma(x = fit)
Value p-value Decision
Global Stat 2.7728 0.5965 Assumptions acceptable.
Skewness 1.5374 0.2150 Assumptions acceptable.
Kurtosis 0.6376 0.4246 Assumptions acceptable.
Link Function 0.1154 0.7341 Assumptions acceptable.
Heteroscedasticity 0.4824 0.4873 Assumptions acceptable.最后的Global Stat是对4个假设条件进行综合验证，通过了即表示4个假设验证都通过了。最后的Heterosceasticity是进行异方差检测。注意这里假设检验的原假设都是假设成立，所以当p>0.05时，假设才能能过验证。
如果综合验证不通过，也有其他方法对4个假设条件分别验证：
线性假设
library(car)
crPlots(fit)
返回的图是各个自变量与残差（因变量）的线性关系图，若存着明显的非线性关系，则需要对自变量作非线性转化。书中说这张图表明线性关系良好。
正态性
library(car)
qqPlot(fit,labels = row.names(states),id.method = "identify",simulate = TRUE,main = "Q-Q Plot")
qqPlot()可以生成交互式的qq图，选中异常点，就返回该点的名称。该图中除了Nevad点，其他点都在直线附近，可见正态性良好。
同方差性
library(car)
ncvTest(fit)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 1.746514 Df = 1 p = 0.1863156p值大于0.05，可认为满足方差相同的假设。
独立性
library(car)
durbinWatsonTest(fit)
lag Autocorrelation D-W Statistic p-value
1 -0.2006929 2.317691 0.258
Alternative hypothesis: rho != 0p值大于0.05，可认为误差之间相互独立。
除了以上4点基本假设，还有其他方面需要进行诊断——
2、多重共线性
理想中的线性模型各个自变量应该是线性无关的，若自变量间存在共线性，则会降低回归系数的准确性。一般用方差膨胀因子VIF(Variance Inflation Factor)来衡量共线性，《统计学习》中认为VIF超过5或10就存在共线性，《R语言实战》中认为VIF大于4则存在共线性。理想中的线性模型VIF=1，表完全不存在共线性。
library(car)
vif(fit)
Population Illiteracy Income Frost
1.245282 2.165848 1.345822 2.082547可见这4个自变量VIF都比较小，可认为不存在多重共线性的问题。
3、异常值检验
离群点
离群点有三种判断方法：一是用qqPlot()画QQ图，落在置信区间（上图中两条虚线）外的即可认为是离群点，如上图中的Nevad点；一种是判断学生标准化残差值，绝对值大于2（《R语言实战》中认为2，《统计学习》中认为3）的可认为是离群点。
plot(x=fitted(fit),y=rstudent(fit))
abline(h=3,col="red",lty=2)
abline(h=-3,col="red",lty=2)
which(abs(rstudent(fit))>3)
Nevada
28还有一种方法是利用car包里的outlierTest()函数进行假设检验：
library(car)
outlierTest(fit)
rstudent unadjusted p-value Bonferonni p
Nevada 3.542929 0.00095088 0.047544这个函数用来检验最大的标准化残差值，如果p>0.05，可以认为没有离群点；若p<0.05，则该点是离群点，但不能说明只有一个离群点，可以把这个点删除之后再作检验。第三种方法可以与第二种方法结合起来使用。
高杠杆点
高杠杆值观测点，即是与其他预测变量有关的离群点。换句话说，它们是由许多异常的预测变量值组合起来的，与响应变量值没有关系。《统计学习》中给出了一个杠杆统计量，《R语言实战》中给出了一种具体的操作方法。（两本书也稍有出入，《统计学习》中平均杠杆值为(p+1)/n，而在《R语言实战》中平均杠杆值为p/n；事实上在样本量n比较大时，几乎没有差别。）
hat.plot <- function(fit){
p <- length(coefficients(fit))
n <- length(fitted(fit))
plot(hatvalues(fit),main = "Index Plot of Hat Values")
abline(h=c(2,3)*p/n,col="red",lty=2)
identify(1:n, hatvalues(fit), names(hatvalues(fit))) #这句产生交互效果，选中某个点后，关闭后返回点的名称
}
hat.plot(fit)
超过2倍或3倍的平均杠杆值即可认为是高杠杆点，这里把Alaska和California作为高杠杆点。
强影响点
强影响点是那种若删除则模型的系数会产生明显的变化的点。一种方法是计算Cook距离，一般来说， Cook’s D值大于4/(n-k -1)，则表明它是强影响点，其中n 为样本量大小， k 是预测变量数目。
cutoff <- 4/(nrow(states-length(fit$coefficients)-2)) #coefficients加上了截距项，因此要多减1
plot(fit,which=4,cook.levels = cutoff)
abline(h=cutoff,lty=2,col="red")作者：真依然很拉风
链接：https://www.jianshu.com/p/4c6c8174f292
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。