残差有很多种类,如普通残差(ordinary residual),皮尔森残差(Pearson residual), 学生化残差(studentized residual)
残差都反应拟合值和观测值之间的不同,是各种回归模型诊断方法的基础统计量
实战
- 导入、查看数据集
library(car)
str(Mroz)

- 拟合模型,并绘制残差图
指导变量种类(哪些变量纳入模型)和形式(线性vs非线性)的选择
mroz.mod <- glm(lfp~k5+k618+age+wc+hc+lwg+inc,family = binomial,Mroz)
residualPlots(mroz.mod)
# Test stat Pr(>|Test stat|)
# k5 0.1157 0.73380
# k618 0.1569 0.69207
# age 1.1892 0.27549
# wc
# hc
# lwg 153.5037 < 2e-16 ***
# inc 3.5463 0.05968 .
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lwg的p值非常小,说明二者之间不是简单的线性关系;从图中可以看到,皮尔森残差和lwg之间存在非线性关系
使用非线性模型进行拟合(二次方程)
mroz.mod2 <- glm(lfp~k5+k618+age+wc+hc+lwg+I(lwg^2)+inc,family = binomial,Mroz)
residualPlots(mroz.mod2)
#此时I(lwg^2)的p值已经较大,说明使用二次方程拟合是合适的
# Test stat Pr(>|Test stat|)
# k5 0.2581 0.61145
# k618 0.1294 0.71910
# age 1.9598 0.16153
# wc
# hc
# lwg 0.0000 1.00000
# I(lwg^2) 0.2663 0.60584
# inc 3.0415 0.08116 .
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

边际图
marginalModelPlots(mroz.mod)

以上均是对模型的整体拟合情况进行评估,还需要对个别观测值进行判断。包括:
Outlier 离群值

outlierTest(mroz.mod)
# No Studentized residuals with Bonferroni p < 0.05
# Largest |rstudent|:
# rstudent unadjusted p-value Bonferroni p
# 119 2.25002 0.024448 NA
Case119的Bonferroni p 不存在,说明不存在离群值
Leverage 杠杆值

influenceIndexPlot(mroz.mod,id.n = 3) # id.n 显示3个最大的值

Influence 影响值

influencePlot(mroz.mod,col = "red",id.n = 3)
# StudRes Hat CookD
# 119 2.250020 0.02779190 0.036422614
# 220 2.234546 0.01963334 0.025498342
# 348 1.283820 0.06455634 0.010957632
# 386 1.224214 0.05792918 0.008554865

Case119很突兀,将其删除后进行模型更新
mroz.mod119 <- update(mroz.mod,subset=c(-119))
compareCoefs(mroz.mod,mroz.mod119)

可以发现上述函数输出结果中,原始模型和新模型回归系数的变化非常小,因此Case119不是影响值,最终模型不需要剔除它。
参考资料
文中代码及部分截图来自章仲恒教授的丁香园课程:回归模型诊断:残差及模型诊断
网友评论