26-R线性回归及回归诊断实例
1、OLS线性回归的基本原则
最优拟合曲线应该使各点到直线的距离的平方和（即残差平方和，简称RSS）最小。
2、OLS线性回归模型假设：
- 正态性：对于固定的自变量值，因变量值成正态分布；
- 独立性：个体之间相互独立;
- 线性相关：因变量和自变量之间为线性相关;
- 同方差性：因变量的方差不随自变量的水平不同而变化，即因变量的方差是不变的。
3、lm()线性回归函数
lm(formula, data)
formula中的符号注释：
~分割符号，左边为因变量，右边为自变量，例如， z~x+y，表示通过x和y来预测z
+分割预测变量
:表示预测变量的交互项，例如，z~x+y+x:y
* 表示所有可能的交互项，例如，z~x*y 展开为 z~x+y+x:y
^表示交互项的次数，例如，z ~ (x+y)^2，展开为z~x+y+x:y
.表示包含除因变量之外的所有变量，例如，如果只有三个变量x，y和z，那么代码 z~. 展开为z~x+y+x:y
-1表示删除截距项，强制回归的直线通过原点
I()从算术的角度来解释括号中的表达式，例如，z~y+I(x^2) 表示的拟合公式是 z=a+by+cx2
function 可以在表达式中应用数学函数，例如，log(z) ~x+y
对于拟合后的模型（lm函数返回的对象），可以应用下面的函数，得到模型的更多额外的信息：
summary() 展示拟合模型的详细结果
coefficients()列出模型的参数（截距项intercept和斜率）
confint()提供模型参数的置信区间
residuals() 列出拟合模型的残差值
fitted() 列出拟合模型的预测值
anova() 生成一个拟合模型的方差分析表
predict() 用拟合模型对新的数据预测响应变量
4、举个例子
4.1 数据
library(pacman)
p_load(stargazer)
x <- c(1.5,2.8,4.5,7.5,10.5,13.5,15.1,16.5,19.5,22.5,24.5,26.5)
y <- c(7.0,5.5,4.6,3.6,2.9,2.7,2.5,2.4,2.2,2.1,1.9,1.8)
4.2 模型选择
4.2.1 直线回归
lm1 <- lm(y~x)
stargazer(lm1, # 模型对象、数据框、向量、矩阵、列表
type = "text", # 可设置为text、html、latex
title = "y ~ x", # 表名
style = "default",
summary = F, # 是否做数据统计
align = T, # 是否居中
no.space = T, # 删除空行
single.row = T,# 使估计量和置信区间并排展示
keep.stat = "adj.rsq", # keep.stat="n"只展示样本量的大小，并移除其他统计量
font.size = "large",
model.names = T,
header = F # 避免输出关于包作者的一些信息
)
## y ~ x
## =======================================
## Dependent variable:
## ---------------------------
## y
## OLS
## ---------------------------------------
## x -0.170*** (0.027)
## Constant 5.603*** (0.435)
## ---------------------------------------
## Adjusted R2 0.776
## =======================================
## Note: *p<0.1; **p<0.05; ***p<0.01
修正后的R2=0.776，不是很理想。
4.2.2 多项式回归
x1=x
x2=x2
x1 <- x
x2 <- x^2
lm2 <- lm(y~x1+x2)
stargazer(lm2,
type = "text",
title = "y ~ x1 + x2",
style = "default",
summary = F,
align = T,
no.space = T,
single.row = T,
keep.stat = "adj.rsq",
font.size = "large",
model.names = T,
header = F
)
##
## y ~ x1 + x2
## =======================================
## Dependent variable:
## ---------------------------
## y
## OLS
## ---------------------------------------
## x1 -0.466*** (0.057)
## x2 0.011*** (0.002)
## Constant 6.915*** (0.332)
## ---------------------------------------
## Adjusted R2 0.941
## =======================================
## Note: *p<0.1; **p<0.05; ***p<0.01
修正后的R2=0.941，比直线回归改善很多。
画个拟合图看看效果：
library(tidyverse)
df1 <- cbind(y=y,x=x) %>% as.data.frame()
df2 <- cbind(y=fitted(lm2),x=x) %>% as.data.frame()
ggplot() +
geom_point(data=df1,aes(x,y)) +
geom_line(data=df2,aes(x,y))
多项式回归
4.2.3 对数回归
x=log(x)
lm3 <- lm(y ~ log(x))
stargazer(lm3,
type = "text",
title = "y ~ log(x)",
style = "default",
summary = F,
align = T,
no.space = T,
single.row = T,
keep.stat = "adj.rsq",
font.size = "large",
model.names = T,
header = F
)
##
## y ~ log(x)
## =======================================
## Dependent variable:
## ---------------------------
## y
## OLS
## ---------------------------------------
## log(x) -1.757*** (0.068)
## Constant 7.364*** (0.169)
## ---------------------------------------
## Adjusted R2 0.984
## =======================================
## Note: *p<0.1; **p<0.05; ***p<0.01
修正后的R2=0.984，比多项式回归又有了改善。
画个拟合图看看效果，可以看到拟合效果更好：
df3 <- cbind(y=fitted(lm3),x=x) %>% as.data.frame()
ggplot() +
geom_point(data=df1,aes(x,y)) +
geom_line(data=df3,aes(x,y))
对数回归
4.2.4 指数回归
y=log(y)
lm4 <- lm(log(y) ~ x)
stargazer(lm4,
type = "text",
title = "log(y) ~ x",
style = "default",
summary = F,
align = T,
no.space = T,
single.row = T,
keep.stat = "adj.rsq",
font.size = "large",
model.names = T,
header = F
)
##
## log(y) ~ x
## =======================================
## Dependent variable:
## ---------------------------
## log(y)
## OLS
## ---------------------------------------
## x -0.049*** (0.005)
## Constant 1.760*** (0.075)
## ---------------------------------------
## Adjusted R2 0.907
## =======================================
## Note: *p<0.1; **p<0.05; ***p<0.01
修正后的R2=0.907，比多项式回归还差。
4.2.5 幂函数回归
y=log(y)
x=log(x)
lm5 <- lm(log(y) ~ log(x))
stargazer(lm5,
type = "text",
title = "log(y) ~ log(x)",
style = "default",
summary = F,
align = T,
no.space = T,
single.row = T,
keep.stat = "adj.rsq",
font.size = "large",
model.names = T,
header = F
)
##
## log(y) ~ log(x)
## =======================================
## Dependent variable:
## ---------------------------
## log(y)
## OLS
## ---------------------------------------
## log(x) -0.472*** (0.012)
## Constant 2.191*** (0.030)
## ---------------------------------------
## Adjusted R2 0.993
## =======================================
## Note: *p<0.1; **p<0.05; ***p<0.01
修正后的R2=0.993，比对数回归又有了改善。
画个拟合图看看效果，可以看到幂函数法拟合效果最好：
df4 <- cbind(y=exp(fitted(lm5)),x=x) %>% as.data.frame()
ggplot() +
geom_point(data=df1,aes(x,y)) +
geom_line(data=df4,aes(x,y))
幂函数回归
5、回归诊断内容
- 样本是否符合正态分布假设？
- 是否存在离群值导致模型产生较大误差？
- 线性模型是否合理？
- 残差是否满足独立性、等方差、正态分布等假设条件？
- 是否存在多重共线性？
回归诊断的总结的函数为：influence.measures()
5.1 正态性检验
- 函数shapiro.test()
- 直方图或散点图目测检验
- 残差计算函数residuals()，对残差作正态性检验（残差散点图）
# 残差正态性检验，P值>0.05，W接近1说明符合正态分布
shapiro.test(lm5$residuals)
##
## Shapiro-Wilk normality test
##
## data: lm5$residuals
## W = 0.93033, p-value = 0.3837
# 回归诊断的标准方法
par(mfrow=c(2,2))
plot(lm5)
回归诊断
正态性：当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正态分布。正态QQ图（Normal Q-Q）是在正态分布对应的值下，标准化残差的概率图。图中的数据点按对角直线排列，趋于一条直线，并被对角直接穿过，直观上符合正态分布。对于近似服从正态分布的标准化残差，应该有 95% 的样本点落在 [-2,2] 区间内；若不是如此，就违反了正态假设。
独立性：无法判断因变量和自变量的值是否相互独立，只能从收集的数据中来验证。
线性相关性：若因变量与自变量线性相关，那么残差值与预测（拟合）值除了系统误差之外，就没有任何关联。在残差和拟合图（Residuals VS Fitted）中，可以看到一条曲线，这暗示着可能需要对拟合模型加上一个二次项。（所以我们最后选择幂函数回归）
同方差性：若满足不变方差假设，那么在位置尺度图（Scale-Location）中，水平线周围的点应该随机分布。
残差和杠杆图（Residuals VS Leverage）提供了特殊的单个观测点（离群点、高杠杆点、强影响点）的信息，本图中出现了红色的等高线，说明数据中存在特别影响回归结果的异常点。残差和杠杆图的可读性差，不够实用。
- 离群点：表明拟合回归模型对其预测效果不佳（产生巨大的残差）
- 高杠杆点：指自变量因子空间中的离群点（异常值），是由许多异常的自变量值组合起来的，与因变量没有关系。
- 强影响点：表明它对模型参数的估计产生的影响过大，非常不成比例。强影响点可以通过Cook距离（Cook distance）统计量来鉴别。看到上面4幅中，每幅图上都有一些点被特别的标记出来了，这些点是可能存在的异常值点，如果要对模型进行优化，我们可以从这些来入手，从图中发现，索引编号为1、3和12的3个点在多幅图中出现。这3个点可能为异常点，可以从数据中去掉这3个点，再进行显著性检验和残差分析。但本次残差分析的结果已经很好了，所以对于异常点的优化，可能并不能明显的提升模型的效果。
5.2 残差正态性检验（qqplot）
car包qqplot()函数提供了精确的正态假设检验方法，置信区间通过虚线划定，当绝大多数点都落在置信区间时，说明正态性假设符合得很好。
5.3 离群值（强影响点）
hv <- hatvalues(lm5);hv
## 1 2 3 4 5 6 7 8 9
## 0.48254231 0.26578892 0.15707660 0.09415762 0.08337304 0.09120380 0.09907002 0.10721053 0.12714202
## 10 11 12
## 0.14898865 0.16407973 0.17936679
# 如何找到最重要或最有影响的观察结果？将hat值切分为四分位数，应用95%标准过滤最异常值，将该过滤标准应用于观察结果
hv.vip <- df4$hv[df4$hv > quantile(hv,0.95)];hv.vip
# 具象化
df5 <- cbind(df1,hv=hv)
plot(df5$x,df5$y,col=ifelse(df5$hv > quantile(hv,0.95),'red','blue'))
离群点
5.4 误差的独立性
car包提供了durbinWatsonTest()函数，用于做Durbin-Watson检验，检测误差的序列相关性。
durbinWatsonTest(lm5)
lag Autocorrelation D-W Statistic p-value
1 0.1840413 1.184471 0.048
Alternative hypothesis: rho != 0
p值<0.05，不显著，误差项之间独立，不存在自相关性。
5.5 线性相关性
通过成分残差图（component + residual plots）检查因变量和自变量之间是否呈线性关系。若图形存在非线性，则说明可能对预测变量的函数形式建模不够充分，那么需要添加一些曲线成分，比如多项式，对一个或多个变量进行变换（log(x)代替x），或用其他回归变体形式而不是线性回归。
crPlots(lm5)
线性相关性
5.6 同方差性
判断方差是否恒定，nvcTest()函数生成一个记分检验，原假设为误方差不变，备择假设为误差方差随着拟合值水平的变化而变化。若检验显著，说明存在误方差不恒定。
spreadLevelPlot()函数创建了一个添加了最佳拟合曲线的散点图，展示了标准化残差绝对值与拟合值的关系。
ncvTest(lm5)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.4044358, Df = 1, p = 0.52481
记分检验不显著：p=0.52481，说明满足方差不变假设，也可以通过分布水平看到这一点，点在水平的最佳拟合曲线周围呈随机分布。
spreadLevelPlot(lm5)
水平分布
5.7 多重共线性（一个自变量不存在）
在回归分析的时候，古典模型中有一个基本的假定就是自变量之间是不相关的，但是如果我们在拟合出来的回归模型出现了自变量之间高度相关的话，可能对结果又产生影响，我们称这个问题为多重共线性。
根据经验表明，多重共线性存在的一个标志就是模型存在较大的标准差和较小的T统计量，如果一个模型的可决系数R2很大，F检验高度限制，但偏回归系数的T检验几乎都不显著，那么模型很可能是存在多重共线性。
- 利用自变量之间的简单相关系数检验，一般而言，如果两个解释变量的简单相关系数较高，则可以认为是存在着严重的多重共线性。
- 随着多重共线性的严重程度增强，方差膨胀因子会逐渐的变大，在回归中我们用VIF表示方差膨胀因子，VIF=1/(1-R2)。一般当VIF>=10的时候，我们就可以认为存在严重多重共线性。
car包中的vif()函数可以帮我们算出这个方差膨胀：
library(car)
vif(lm5)
存在多重共线性的解决方法之一：利用MASS包中的函数lm.ridge()来建立岭回归模型。
6、回归诊断（gvlma包）
gvlma()函数，用于对线性模型假设进行综合验证，同时还能验证偏斜度，峰度和异方差的评价。从输出项中可以看出，假设检验的显著性水平是5%。如果p<0.05，说明违反了假设条件。从Global Stat 的Decision 栏中，可以看到数据满足OLS回归模型所有统计假设的P值（Assumptions acceptable）。
library(gvlma)
gvmodel <- gvlma(lm5)
summary(gvmodel)
##
## Call:
## lm(formula = log(y) ~ log(x))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.054727 -0.020805 0.004548 0.024617 0.045896
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.19073 0.02951 74.23 4.81e-15 ***
## log(x) -0.47243 0.01184 -39.90 2.34e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0361 on 10 degrees of freedom
## Multiple R-squared: 0.9938, Adjusted R-squared: 0.9931
## F-statistic: 1592 on 1 and 10 DF, p-value: 2.337e-12
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = lm5)
##
## Value p-value Decision
## Global Stat 8.46263 0.076028 Assumptions acceptable.
## Skewness 0.26282 0.608186 Assumptions acceptable.
## Kurtosis 0.52317 0.469492 Assumptions acceptable.
## Link Function 7.63997 0.005709 Assumptions NOT satisfied!
## Heteroscedasticity 0.03666 0.848150 Assumptions acceptable.
