R语言实战8：回归

作者: 冬之心 | 来源:发表于2019-10-04 23:19 被阅读0次

title: "lesson 8 回归"
author: "wintryheart"
date: "2019/9/9"
output:
html_document:
toc: TRUE
number_sections: TRUE

knitr::opts_chunk$set(echo = TRUE)

研究的三个困难

发现一个有趣的问题
设计一个有用的、可以测量的响应变量
收集合适的数据

OLS回归

统计假设：

正态性：对于固定的自变量值，因变量值成正态分布；
独立性： $Y_i$ 值之间相互独立；
线性
同方差性：因变量的方差不随自变量的水平不同而变化。

lm()拟合模型

常用符号	用途
~	分隔符，左边为因变量，右边为自变量
+	分隔自变量
:	交互项
*	所有可能交互项
^	交互项达到某个次数
.	除因变量外的所有变量
-	从等式中移除某个变量
-1	删除截距项
I()	从算术角度解释括号中的元素
function	可用在表达式中的函数

常用函数	用途
summary()	详细的模型结果
coefficeints()	模型参数（截距和斜率）
confint()	参数的置信区间
fitted()	预测值
residuals()	残差值
anova()	方差分析表
vcov()	参数的协方差矩阵
AIC()	赤池信息统计量
plot()	拟合评价诊断图
predict()	基于拟合模型对新数据集预测

简单线性回归


fit <- lm(weight ~ height, data=women)

summary(fit)

fitted(fit)

residuals(fit)

plot(women$height, women$weight)

abline(fit)

unnamed-chunk-1-1.png

多项式回归


fit2 <- lm(weight ~ height + I(height^2), data=women)

summary(fit2)

car包scatterplot()函数绘制二元关系图

library(car)
scatterplot(weight~height,data=women, spread=F, smoother.args=list(lty=2), pch=19,
            main="Women Age 30-39", xlab="hight", ylab="weight")

unnamed-chunk-3-1.png

spread=FALSE，删除残差正负均方根在平滑曲线上的展开和非对称信息。

smoother.args=list(lty=2),设置loess拟合曲线为虚线。

多元回归

car包scatterplotMatrix()生成散点图矩阵。
scatterplotMatrix()默认在非对角线区域绘制变量间散点图，并添加平滑和线性拟合曲线。
对角线区域绘制每个变量的密度图和轴须图。


states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy","Income","Frost")])

cor(states)

library(car)

scatterplotMatrix(states, spread=FALSE, smoother.args=list(lty=2), main="Scatter Plot Matrix")

fit3 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)

summary(fit3)

unnamed-chunk-4-1.png

交互项


fit4 <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)

summary(fit4)

effects包专门用于显示线性模型、一般线性模型和其他许多模型的作用效应。
effects包中的effect()函数，可以用图形展示线性模型的交互项的结果。
plot(effect(term, mod, ,xlevels), multiline=TRUE)


library(effects)
plot(effect("hp:wt", fit4, , list(wt=c(2.2, 3.2, 4.2))))
plot(effect("hp:wt", fit4, , list(wt=c(2.2, 3.2, 4.2))), multiline=TRUE)

unnamed-chunk-6-1.png

unnamed-chunk-6-2.png

回归诊断

# 置信区间
confint(fit3)

结果表明，文盲率改变1%，谋杀率在95%的置信区间[2.38, 5.90]中变化。

因为Frost的置信区间包含0，所以可以得出结论：当其他变量不变时，温度的改变与谋杀率无关。

标准方法

对lm()函数返回的对象使用plot()函数，可以生成评价模型拟合情况的四幅图形。

正态Q-Q图，检验正态性。如果满足正态假设，图上的点应该落在呈45度角的直线上。
残差图与拟合图，检验线性。如果因变量与自变量线性相关，残差与预测值应该没有任何系统关联。
位置尺度图（Scale-Location），检验同方差性。如果满足不变方差假设，水平线周围的点应该随机分布。
残差与杠杆图（Residuals VS Leverage），检验离群点、杠杆值和强影响点。

fit <- lm(weight ~ height, data=women)
par(mfrow=c(2,2))
plot(fit)

#残差拟合图显示一个曲线关系，暗示需要对回归模型加上一个二次项。
fit2 <- lm(weight ~ height + I(height^2), data = women)
plot(fit2)

fit3 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
plot(fit3)

unnamed-chunk-8-1.png

unnamed-chunk-8-2.png

unnamed-chunk-8-3.png

改进方法

car包中的回归诊断实用函数

函数	目的
qqPlot()	分位数比较图
durbinWatsonTest()	对误差自相关做Durbin-Watson检验
crPlots()	成分与残差图
ncvTest()	对非恒定的误差方差做得分检验
spreadLevelPlot()	分散水平检验
outlierTest()	Bonferroni离群点检验
avPlots()	添加的变量图形
influencePlot()	回归影响图
scatterplot()	增强的散点图
scatterplotMatrix()	增加的散点力矩阵
vif()	方差膨胀因子

正态性检验qqPlot()


states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy","Income","Frost")])
fit3 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
par(mfrow=c(1,1))
qqPlot(fit3, labels=row.names(states), simulate=TRUE, main="Q-Q Plot")
states["Nevada",]
fitted(fit3)["Nevada"]
residuals(fit3)["Nevada"]
rstudent(fit3)["Nevada"]

unnamed-chunk-9-1.png

Default S3 method:

qqPlot(x, distribution="norm", groups, layout,
ylim=range(x, na.rm=TRUE), ylab=deparse(substitute(x)),
xlab=paste(distribution, "quantiles"), glab=deparse(substitute(groups)),
main=NULL, las=par("las"),
envelope=.95, col=carPalette()[1], col.lines=carPalette()[2],
lwd=2, pch=1, cex=par("cex"),
line=c("quartiles", "robust", "none"), id=TRUE, grid=TRUE, ...)

S3 method for class 'formula'

qqPlot(formula, data, subset, id=TRUE, ylab, glab, ...)

S3 method for class 'lm'

qqPlot(x, xlab=paste(distribution, "Quantiles"), ylab=paste("Studentized Residuals(",deparse(substitute(x)), ")", sep=""),
main=NULL, distribution=c("t", "norm"),
line=c("robust", "quartiles", "none"), las=par("las"),
simulate=TRUE, envelope=.95, reps=100,
col=carPalette()[1], col.lines=carPalette()[2], lwd=2, pch=1, cex=par("cex"),
id=TRUE, grid=TRUE, ...)

distribution:
root name of comparison distribution – e.g., "norm" for the normal distribution; t for the t-distribution.

envelope:
confidence level for point-wise confidence envelope, or FALSE for no envelope. 置信区间

id:
controls point identification; if FALSE, no points are identified; can be a list of named arguments to the showLabels function; TRUE is equivalent to list(method="y", n=2, cex=1, col=carPalette()[1], location="lr"), which identifies the 2 points with the 2 points with the most extreme verical values — studentized residuals for the "lm" method. Points labels are by default taken from the names of the variable being plotted is any, else case indices are used. Unlike most graphical functions in car, the default is id=TRUE to include point identification.

simulate:
if TRUE calculate confidence envelope by parametric bootstrap; for lm object only. The method is due to Atkinson (1985). 参数自助法生成置信区间

误差的独立性检验，Durbin-Watson检验

durbinWatsonTest(fit3)

p值不显著，说明无自相关性，误差项之间独立。

线性检验，crPlots()

通过成分残差图（component plus residual plot），或偏差图（partial residual plot），看因变量和自变量之间是否呈非线性关系。

crPlots(fit3)

unnamed-chunk-11-1.png

成分残差图说明，线性模型形式是合适的。

若图形存在非线性，说明线性建模不够充分，需要添加一些曲线成分（如多项式），或对变量进行变换（如log变换），或用其它回归变体形式而不是线性回归。

同方差性检验，判断误差方差是否恒定

ncvTest()生成一个计分检验，零假设为误差方差不变。若检验显著，说明存在异方差性（误差方差不恒定）。

spreadLevelPlot()创建一个最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系。


ncvTest(fit3)
spreadLevelPlot(fit3)

unnamed-chunk-12-1.png

ncvTest()计分检验不显著，说明满足同方差性假设。

spreadLevelPlot()建议幂次变换（suggested power transformation)，即经过p次幂变换（ $y^p$ ），非恒定的误差方差将会平稳。。由于建议幂次接近1，异方差性不明显，不需要进行变换。

线性模型假设的综合验证

gvlma包中的gvlma()函数，对线性模型假设进行综合验证，同时还能做偏斜度、峰度和异方差性的评价。


library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)

Decision列表明是否违反假设条件。

多重共线性

多重共线性用统计量VIF（variance inflation factor，方差膨胀因子）进行检测。

VIF的平方根表示，变量回归参数的置信区间能膨胀为与模型无关的预测变量的程度。

car包中的vif()函数提供VIF值。一般原则下， $\sqrt{vif}>2$ 就表明存在多重共线性问题。

library(car)

vif(fit3)

sqrt(vif(fit3))>2

异常值

离群点

离群点指那些模型预测不佳的观测点。

Q-Q图中，落在置信区间带外的点，即可被认为是离群点。
标准化残差的绝对值大于2的点可能是离群点。
car包的outlierTest()函数，根据单个绝对值最大的残差值的显著性来判断是否有离群点。

library(car)
fit5 <- lm(mpg~ wt+qsec, data = mtcars)
outlierTest(fit5)

高杠杆值点

高杠杆值点指那些与其他预测变量有关的离群点，即它们是许多异常的预测变量值组合起来的，与响应变量的值无关。

高杠杆值点可以通过帽子统计量(hat statistic)判断。帽子均值为 $p/n$ ，其中p为模型估计的参数数目（包含截距项），n是样本量。
如果观测点的帽子值大于帽子均值的2或3倍，就可以认定为高杠杆值点。

hatvalues(fit3)

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(x=1:n, y=hatvalues(fit), labels=names(hatvalues(fit)))
}
hat.plot(fit3)

unnamed-chunk-16-1.png

强影响点

强影响点指对模型参数估计值影响有些比例失衡的点。即移除该观测点，模型会发生重大改变。

有两种方法可以检测：Cook距离或称D统计量，以及变量添加图（added variable plot）。

Cook距离

一般来说，Cook's D值大于4/(n-k-1)则表明它是强影响点，其中n为样本量，k为预测变量数。


cutoff <- 4/(nrow(states)-length(fit3$coefficients)-2)
plot(fit3, which = 4, cook.levels = cutoff)
abline(h=cutoff, lty=2, col="red")

unnamed-chunk-17-1.png

其实以1为分割点，比以4/(n-k-1)为分割点更具有一般性。

变量添加图

变量添加图，即对于每个预测变量 $X_k$ ，绘制 $X_k$ 在其他k-1个预测变量上回归的残差值相对于响应变量在其他k-1预测变量上的回归残差值的关系图。

car包中的avPlots()函数，可以作变量添加图。

library(car)

avPlots(fit3, ask=FALSE, id.method="identify")

unnamed-chunk-18-1.png

car包中的influencePlot()函数，可以把离群点、杠杆值和强影响点的信息整合到一幅图形中。


influencePlot(fit3, id.method="identify", main="Influence Plot",
              sub="Circle size is proportional to Cook's distance")

unnamed-chunk-19-1.png

其中，

纵坐标超过+2或小于-2的点，可被认为离群点。
水平坐标超过0.2或0.3的点，可被认为高杠杆值。
圆圈大小与影响成比例，圆圈很大的点可能是强影响点。

改进措施

删除观测点

变量变换

car包

powerTransform()通过 $\lambda$ 的最大似然估计来正态化变量 $X^\lambda$ 。
boxTidwell()通过获得预测变量幂数的最大似然估计来改善线性关系。
spreadLevelPlot()提供幂次变换应用。


library(car)
summary(powerTransform(states$Murder))

结果表明，可用 $\lambda=0.6$ 来正态化变量Murder；不过 $\lambda=1$ 的假设也无法拒绝。因此没有强有力的证据表明需要变换变量。

boxTidwell(Murder ~ Population + Illiteracy, data=states)

结果显示，使用 $Population^{0.87}$ 和 $Illiteracy^{1.36}$ 能够改善线性关系。但是检验结果又表明变量并不需要变换。

如果变量变换得没有意义，应该避免这样做。

增删变量

尝试其他方法

如果存在离群点或强影响点，可以使用稳健回归；
如果违背了正态性假设，可以使用非参数回归；
如果存在显著的非线性，可以尝试非线性回归；
如果违背误差独立性假设，可以尝试专门研究误差结构的模型，如时间序列模型或多层回归模型。
广义线性回归模型，适用于许多OLS回归假设不成立的情况。

选择最佳模型

模型比较

anova()函数可以比较两个嵌套模型的拟合优度。


states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy","Income","Frost")])

fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)

fit2 <- lm(Murder ~ Population + Illiteracy , data=states)

anova(fit2, fit1)

结果表明，模型预测效果差异并不显著，Income和Frost两个变量没必要添加到模型中。

AIC()函数可以比较两个模型的AIC。AIC值越小越好。

AIC(fit1, fit2)

注意，anova()需要嵌套模型，而AIC()不需要。

变量选择

anova()和AIC()只能比较两个模型。
如果有更多的模型，需要用以下方法。

逐步回归

MASS包中的stepAIC()函数可以实现逐步回归模型，依据的是精确AIC准则。

library(MASS)
states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy","Income","Frost")])

fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
stepAIC(fit, direction="backward")

逐步回归法可能会找到一个好的模型，但是不能保证模型就是最佳的，因为不是每一个可能的模型都被评价。

全子集回归

leaps包中的regsubsets()函数，所有可能的模型都会被检验。

通过R方，调整R方或Mallows Cp统计量等准则选择最佳模型。

对于一个好的模型，它的Cp统计量非常接近于模型的参数数目（包括截距项）。

结果可用leaps包中的plot()绘制，也可用car包中的subsets()绘制。

library(leaps)
states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy","Income","Frost")])
leaps <- regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data=states, nbest=4)
plot(leaps, scale = "adjr2")

library(car)
subsets(leaps, statistic="cp", legend=c(3.5, 50), main="Cp Plot for All Subsets Regression")
abline(1, 1, lty=2, col="red")

subsets(leaps, statistic="bic", legend=c(3.5, 0), main="Cp Plot for All Subsets Regression")
subsets(leaps, statistic="adjr2", legend=c(3.5, 0.1), main="Cp Plot for All Subsets Regression")

unnamed-chunk-25-1.png

unnamed-chunk-25-2.png

unnamed-chunk-25-3.png

unnamed-chunk-25-4.png

在subsets()绘制的Cp图中，越好的模型离截距项和斜率均为1的直线越近。

深层次分析：相对重要性

最简单的方法是比较标准化的回归系数。

利用scale()函数将数据标准化为均值为0，标准差为1的数据集。
（注意，scale()返回的是矩阵，而lm()要求的是数据集，需要转换一下。）

states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy","Income","Frost")])

zstates <- as.data.frame(scale(states))

zfit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=zstates)

coef(zfit)

round(coef(zfit), digits = 2) #系数保留两位小数

从结果可以看出Illiteracy的标准化系数最大（r round(coef(zfit), digits = 2)[3]），可以被认为是最重要的变量。

还有其它方法可以定量分析相对重要性，如可以将相对重要性看作每个预测变量（本身或与其他预测变量组合）对R方的贡献。relaimpo包涵盖了一些相对重要性的评价方法。

library(relaimpo)

linmod <- lm(Fertility~Agriculture+Examination+Education+Catholic+Infant.Mortality,swiss)
crlm <- calc.relimp(linmod, type = c("lmg", "last", "first", "betasq", "pratt", "genizi", "car"), rela = TRUE )
crlm
plot(crlm)

unnamed-chunk-27-1.png

calc.relimp()计算线性模型的几种相对重要性指标。推荐使用lgm。该方法将R2。

lmg 将R2分解为每个自变量的贡献，参见 Lindeman, Merenda and Gold 1980, p.119ff ，Chevan and Sutherland (1991)和《线性回归模型中自变量相对重要性的衡量》(孙红卫，王玖，罗文海,2012)
pmvd 是比例边际方差分解，可以解释为带样本权重的加权自变量贡献。参见Feldman (2005)
last 是每一个变量最后被包括进来时的贡献，有时也称有用度。
first 是每一个变量最初被包括进来时的贡献，也即是y和该自变量的协方差平方。
betasq 是标准化系数的平方
pratt 是标准化系数和相关系数的乘积
genizi 是根据Genizi 1993的R2分解
car 是根据Zuber and Strimmer 2010的R2分解。

calc.relimp 中有五个指标(lmg, pmvd, pratt, genizi and car)分解R2。当rela=FALSE，它们的值相加等于R2；当rela=TRUE，它们的值相加等于100(%)