第六章 线性回归分析
内容梳理
TLDR记录
- 线性回归分析的目的是探索一个因变量(结果变量)与一个或多个自变量(解释变量)的关系。当自变量只有一个时,称为简单线性回归;当自变量有多个时,称为多重线性回归;探索两个连续型自变量在第三个分类变量不同类别下的关系时,成为分层线性回归。
- 线性回归的因变量是连续型变量,适用于变量间线性关系的解释及此基础上的预测问题。
- 线性回归分析的步骤一般有三个,一是建立线性回归模型(linear models),即拟合,二是评价模型,
summary()就是最基础的模型评价,做t检验(回归分析),做F检验(残差分析,方差分析),准确理解“回归诊断”的含义。
1. 简单线性回归
- 因变量Y只受一个自变量X影响,存在着近似的线性函数关系。
- 构建模型
- 假定的线性回归模型 = 线性变化(斜率,截距)+ 随机误差(假设服从正态分布)。
- 通常用最小二乘法估计拟合
lm()。
- 查看评估拟合结果:
- 回归分析
summary(lm()):斜率(回归系数),截距,p值(t检验); - 方差分析
aov():决定系数R^2,调整的决定系数Adj R^2,残差分析:F检验
- 回归分析
- 回归诊断:
- 提取样本点的残差值
residuals():拟合值(期望值)与实际值的差距即残差值- 检验残差正态性
- 直方图
hist() - 正态QQ图
qqnorm()``qqline():期望的标准正态分值与残差的散点图;是否按直线分布 - Shapiro-Wilk检验
shapiro.test():定量检验残差的正态性,零假设是给定的数据服从正态分布 - 残差诊断图
plot(lm())
- 直方图
- 检验残差正态性
- 提取样本点的残差值
# 简单线性回归
load("UCR.rdata")
library(epiDisplay)
des(UCR)
summary(UCR)
plot(ucr ~ age, data = UCR)
# 1. 构建线性回归模型
mod <- lm(ucr ~ age, data = UCR)
mod
# mod对象有12个属性,可以分别调用
attributes(mod)
mod$fitted.values
# 2. 模型数据汇总显示
summary(mod)
# 解读:
# Call:模型调用函数;拟合的是线性函数,即ucr = A * age + B
# Residuals:残差分布;本例数据近似对称分布
# 2.1 回归分析:
# 截距B即常数项的估计1.45492 表示年龄为 0 时的尿肌酐含量,0岁的值在真实世界没有实际意义,即使p=2.87e-6<0.05。
# 斜率A即age值的估计0.14869,表示年龄每增长一岁,ucr值平均增加0.14869,其p=7.60e-7<0.05,ucr值与0有显著差异。
# 2.2 方差分析;
# 对残差的进一步详细描述:
# 决定系数R^2的值为 0.7921,表示数据中有 79.2% 的变异能被该模型解释;调整后的决定系数的值为 0.7791。
# 采用 F 统计量对变量 age 的影响做了假设检验。该检验的 p 值(7.597e-07),约等于上面对回归系数的 t 检验的 p 值(7.60e-7)。
# 回归分析和方差分析给出了相同的结论,即年龄与尿肌酐含量之间有显著的线性关系。
# 3. 回归诊断
# 目的:检验残差正态性
# 散点图、直方图(略,定性)
abline(mod)
res <- residuals(mod)
hist(res)
# Q-Q图:期望的标准正态分布于残差的散点图。分布在直线上说明残差呈正态分布(定性)。
qqnorm(res)
qqline(res)
# (定量)假设残差符合正态分布,shapiro-test检验该假设,shapiro-wilk检验的零假设是给定的数据服从正态分布。
shapiro.test(res)
# p=0.9888>0.05,残差与正态分布差异不显著。残差符合正态分布。
# 4. 残差分析:模型自带的残差诊断图
par(mfrow = c(2,2)) # 画布设置为2*2
plot(mod)
par(mfrow = c(1,1))# 画布设置为1*1
# 绘制四幅诊断图,理解如下:
# 4.1 Residuals vs Fitted残差与拟合:横轴为拟合的Y值,纵轴为残差。如果残差分布比较均匀(点均匀分布在纵轴=0线上下),说明模型的误差是随机分布的(即符合Guaasian-Markov条件);如果残差随着因变量的增加而增大(或减小)或者呈现二次曲线分布,说明原始数据可能不是线性关系。有两种处理方式:1先通过对数、指数或平方根等变换,然后再进行线性回归;2选用加权线性回归的方法进行分析。
# 4.2 Q-Q Residuals残差值Q-Q图:检验残差是否符合正态分布,若散点分布在直线上,残差呈正态分布。
# 4.3 Scale-Location位置尺度图:检验残差的方差齐性,与第一张图类似,但将残差进行了标准化处理,以更好地观察残差的偏离情况。
# 4.4 Residuals vs Leverage残差杠杆效力图:判断极端值。横轴为杠杆效力(leverage),指的是某个观测值对于回归模型的影响(即对决定系数 R^2的影响)。杠杆效力越大,说明该观测值对回归模型的影响越强。纵轴为残差,用于表示模型与观测值的实际匹配程度,残差绝对值越大,说明模型预测值与实际观测值相差越远,提示该观测值有欠拟合的风险。杠杆效力与残差的阈值通常用Cook’s distance来判断(图中标注为虚线)。通常Cook’s distance的阈值为4/n(n为样本量)。当一个观测值有很强的杠杆效力同时残差绝对值也很大时,提示此观测值可能是极端值,会对模型的稳定性造成较大的影响,此时需要对此观测值进行检查,看是否存在数据录入错误或者其它造成该数据出现极端状况的因素。如果数据异常情况无法解释,可以考虑将异常值删除再重新建模。
# 图结论:模型拟合效果比较理想,满足线性模型的假设条件。
# 4.5 综合来说
# 诊断图一--残差分布的正态QQ图:如果该图中的散点图聚集在一条直线上,表明残差呈正态分布。
# 诊断图二--残差与拟合值的散点图:观察残差的分布,查看残差的随机性。
# 诊断图三--位置尺度图:检验残差的方差齐性。
# 诊断图四--残差杠杆图:鉴别离群点、高杠杆值点、强影响点。
# 5. 结论:尿肌酐含量与儿童的年龄有关联,年龄每增长一岁,尿肌酐含量平均增加 0.14869 mmol。除去年龄这个因素,影响尿肌酐含量的其他因素可能是随机误差或是其他未考虑到的影响因素。
# 6. 番外
# 第四张图中仅标注了Cook’s distance为0.5和1的情况,不足以辅助我们进行模型判断,此时我们可以使用cooks.distance()函数进行更详细的探索。
# 基于模型计算Cook's distance
cooksD_mod <- cooks.distance(model=mod) # 返回一个数值向量# 可视化Cook's distance
plot( # 对向量进行可视化,横轴为索引值,纵轴为数值cooksD_mod, # 目标向量main = "Cook's distance" # 设置图片名称
)
abline(h=4/length(cooksD_mod), lty=2, col="red") # 添加Cook's distance水平阈值线# 一共有2个超过Cook's distance的观测值
# 如果想要进一步研究这些观测值,可以将它们提取出来
# df_darwin_extract <- df_darwin[cooksD_mod >= 4/length(cooksD_mod),]
# 如果想要移除这些值,可以进行如下操作
# df_darwin_filter <- df_darwin[cooksD_mod < 4/length(cooksD_mod),]
2. 分层线性回归
- 探索两个变量在第三个分类变量的不同类别之间的关系
- 检验两条(或多条)回归直线是否平行;若平行,再检验截距是否相等
- 构建模型
- 通常用最小二乘法估计拟合
lm(因变量 ~ 自变量1 + 自变量2)。
- 通常用最小二乘法估计拟合
- 查看拟合结果
- 绘制拟合直线
- 交互作用项有无统计学意义
- 模型间的chi-square比较
anova()- 当P值显著时,选择自变量数目更多的模型;反之,选择自变量数目更少的模型。
# 分层线性回归
## 1. 假定不管正常或疾病组,年龄对UCR的影响是恒定的,那么两组ucr ~ age的斜率是一样的
mod1 <- lm(ucr ~ age + group, data = UCR)
summary(mod1)
# 年龄每增长1岁,ucr增加0.16156mmol(p=1.83e-07<0.001);在相同的年龄,患病儿童的ucr比正常儿童低0.23256mmo1(p=0.0373)。
# 绘制散点图
col <- ifelse(UCR$group == "Normal", "blue", "red")
pch <- ifelse(UCR$group == "Normal", 1, 19)
plot(ucr ~ age, data = UCR,xlab = "Age in years", ylab = "Urine creatinine (mmol)",col = col, pch = pch)
legend("topleft",legend = c("Normal children", "Diseased children"),col = c("blue", "red"), pch = c(1,19))
# 绘制(平行的)回归直线
# coef()提取模型系数
coef(mod1)
# Intercept = 1.44893 <=> 当age = 0且group = Control(对照组)时,UCR的预测值为1.44893=>control组的截距
# age = 0.16156 <=> 年龄每增加1岁,UCR平均增加0.16156个单位(Control和Patient组一样)=>两条直线的斜率
# groupPatient = -0.23256 <=> 在相同年龄条件下,Patient患者组的UCR比Control对照组平均低0.23256个单位=>patient组与Control组的截距差
b <- coef(mod1)[2];b
a0 <- coef(mod1)[1];a0
a1 <- coef(mod1)[1] + coef(mod1)[3];a1
abline(a = a0, b = b, col = "blue")
abline(a = a1, b = b, col = "red")# 2. 假定正常或疾病组,年龄对UCR的影响是不同的.
# lm(ucr ~ age*group, data = UCR) <=> lm(ucr ~ age + group + age:group, data = UCR)
mod2 <- lm(ucr ~ age*group, data = UCR)
coef(mod2)
# 根据变量定义,ucr = 1.662 + 0.139*age - 0.566*group + 0.033*age*group
# 根据上式,可以推算相应斜率和截距,group=1 <=> group = patient
a0 <- coef(mod2)[1];a0 # normal截距,age=group=0
a1 <- coef(mod2)[1] + coef(mod2)[3];a1 # patient截距,age=0,group=1
b0 <- coef(mod2)[2];b0 # normal斜率,group=0时age的系数
b1 <- coef(mod2)[2] + coef(mod2)[4];b1 # patient斜率,group=1时age的系数
col <- ifelse(UCR$group == 0, "blue", "red")
pch <- ifelse(UCR$group == 0, 1, 19)
plot(ucr ~ age, data = UCR,xlab = "Age in years", ylab = "Urine creatinine (mmol)",col = col, pch = pch)
legend("topleft",legend = c("Normal children", "Diseased children"),col = c("blue", "red"), pch = c(1,19))
abline(a = a0, b = b0, col = "blue")
abline(a = a1, b = b1, col = "red")
summary(mod2)
# age:groupPatient项的p=0.40>0.05,差异无显著意义,说明这个不是斜率不同的原因
# 说明这个二次项没必要,综合还是选择mod1模型
# anova,F检验
anova(mod1, mod2)
# p值=0.40>0.05,两个模型差异不显著,说明模型2没比模型1好,说明二次项没有必要
# 和前面的结论一致
3. 多重线性回归
- 自变量多于一个
- 构建模型
- 假定的线性回归模型 = 线性变化(偏回归系数,常数项)+ 随机误差(假设服从正态分布)。
- 通常用最小二乘法估计拟合
lm()。
- 常用数据分析
- 多重共线性:自变量之间存在相关性,使模型估计失真或难以估计准确
- 查看多个自变量之间的相关性
cor() - 方差膨胀因子VIF:自变量之间的共线性程度
1/(1-summary(lm())$r.squared)``car::vif() - 解决办法
- 剔除造成共线性的自变量
drop(),重新建立模型; - 采用主成分回归法,将一组具有共线性的自变量合并成少数不相关的变量;
- 采用逐步回归法,限制有较强相关性的自变量同时进入回归方程。
- 剔除造成共线性的自变量
- 查看多个自变量之间的相关性
- 多重共线性:自变量之间存在相关性,使模型估计失真或难以估计准确
- 逐步回归,
drop1()``step()- 检验变量对于模型有无统计学意义,需要F检验
anova()
- 检验变量对于模型有无统计学意义,需要F检验
- 回归诊断,
gvlma::gvlma()- 如果不满足假设,需要分别分析各残差。方法同上。
# 多重线性回归
library(ISwR)
data("cystfibr")
?cystfibr
str(cystfibr)
# 总共10个变量
cystfibr$sex <- factor(cystfibr$sex, labels = c("male", "female"))
cor(cystfibr[, 6:10])
# 后5个参数变量间相关系数很多都0.6以上,相关性中等以上
# 主观判断后5个变量之间相关性大,相互之间可以互相表达,所以需建立模型寻求这些变量与前面变量之前的线性相关模型。
# 先观察后5个变量中的fev1与前5个变量之间的关系
# 这个模型可以用最小二乘法拟合模型,fit1-> fev1 ~ 参数[1:5]
fit1 <- lm(fev1 ~ age + sex + height + weight + bmp, data = cystfibr)
summary(fit1)
# 仅sex的p<0.05,仅sex性别对应的t检验有统计学意义。
# F检验的p值0.01209<0.05,是显著的。决定系数0.5123,表明fev1变异的51%可由上述5个自由量的变化来解释。
# 综合上述判断,前5个变量中只有sex对fev1的影响有统计学意义,需要比较其他自由变量的线性关系,即多重共线性问题。
# 先观察身高体重BMI值(bmp)&height&weight。结合通识背景知识,这3个参数的关联性很强,cor()也能验证
cor(cystfibr[,3:5])
# 0.4 0.6 0.9三个相关系数,说明可能存在多重共线性问题
# 多重共线性是指自变量之间由于存在相关性而使模型估计失真或难以估计准确。方差膨胀因子VIF常用于表征自变量之间的共线性程度。一般认为,VIF>5(R^2>80%)则存在多重共线性,用最小二乘法估计的回归系数不准确。
# 注意是自变量之间的,所以另外建立自变量之间的线性模型lm
# 方法1: 手搓一个自变量age做因变量得到的线性模型,求取R-squared并计算vif,得到age与其他自变量之间的vif值
lm.age <- summary(lm(age ~ sex + height + weight + bmp, data = cystfibr))
str(lm.age)
1/(1 - lm.age$r.squared) # vif = 12.71,存在多重共线性
# 这是age与其他自变量的,其他变量类似。
# 方法2: 用R包一次性获得所有自变量的vif值。即直接送入原模型,自动计算所有自变量之间的vif值
library(car)
vif(fit1)
# 结果可以看到,fit1有多个变量存在多重共线性问题(>5)。
# 解决办法:(1)剔除造成共线性的自变量,重新建立模型;(2)采取主成分回归法,将一组具有共线性的自变量合并成少数不相关的变量;(3)采用逐步回归法,限制有较强相关性的自变量同时进入回归方程。这里用逐步回归法,使用到drop1()、step()等函数
# drop1函数会每次剔除一个自变量进行逐步回归,根据赤池信息量准则AIC,AIC越小,表明模型获得了足够的拟合优度。
drop1(fit1, test = "F")
# 根据输出信息,剔除age后的AIC为111.83最小,且剔除前后无统计学差异;剔除sex后的AIC为120.36,剔除前后有统计学差异。
#剔除法的核心思想是剔除那个对模型“最不重要”的变量。而“最不重要”的体现,就是当剔除它时,对模型造成的“损害”最小,即差异最小(不明显),可以通过p值来衡量的,p值越大,差异越小,越先被剔除。以此类推,不断重建模型,计算AIC值和p值。
# 也可以使用step函数,这个函数会根据AIC值进行逐步回归自动自动选择“最优”模型,参数direction可以设定为默认的“both”(向前向后法)、“backward”(向后法)或“forward”(向前法)。
fit2 <- step(fit1)
# 初始AIC是113.81,第一步测试了单独剔除每个自变量后的AIC,哪个AIC最低就剔除哪个,以此类推。直到剔除也降低不了AIC就不剔除了。本例中最终只留下自变量sex+bmp,得到模型fit2
summary(fit2)
# fit2的决定系数(Multiple R-squared)0.507与fit1的决定系数0.512很接近,但fit2调整后的决定系数Adj R^2更高,说明在减少3个变量情况下,fit2具有更好的解释效率。
# (summary()函数提供的)t检验结果表明回归系数与0的差异有无统计学意义。
# 要检验变量对于模型有无统计学意义,尤其自变量是多个水平的分类变量时,需要用F检验(anova())。
anova(fit2)
# 结论是都有显著意义。因为去掉任何一个变量对模型都产生显著影响。
# 另外,还需要对模型做回归诊断,判断模型是否满足最小二乘法的统计假设
library(gvlma)
gvlma(fit2)
# 上述输出:全局统计量Global Stat以及其他各项检验的p>0.05,差异不显著,接受假设,满足线性回归所有的统计假设。# 假如decision出现了不接受,就需要残差分析方法来判断哪些假设没有被满足。
# 残差分析方法步骤如下;
res <- residuals(fit2)
# 先绘制Q-Q图,在斜线上则满足正态分布
qqnorm(res)
qqline(res)
# 定量进行正态假设检验,shapiro-test
shapiro.test(res)
# p=0.1012>0.05,残差与正态分布差异不显著。残差符合正态分布。
# 模型自带的残差诊断图
par(mfrow = c(2,2))
plot(mod)
par(mfrow = c(1,1))
# 诊断图解读(略)# 汇总模型主要结果
library(epiDisplay)
regress.display(fit2)#保存结果
write.csv(regress.display(fit2)$table,file='mytable.csv')
习题
6-1
ISwR 包里的数据集 thuesen 包含了 24 名糖尿病患者的空腹血糖含量(mmol/L)和心室收缩速度(%/s)的测量数据,请运用线性回归分析两个变量之间的联系(注意数据中的缺失值)。
library(ISwR)
data("thuesen")
str(thuesen)
summary(thuesen)
# 缺失值用平均值替代
short.velocity.mean <- mean(thuesen$short.velocity, na.rm = TRUE)
thuesen.amend <- thuesen
thuesen.amend$short.velocity[is.na(thuesen$short.velocity)] <- short.velocity.mean
summary(thuesen.amend)
# 构建模型
fit6.1 <- lm(blood.glucose ~ short.velocity, data = thuesen.amend)
summary(fit6.1)
# p值=0.0436<0.05,差异显著,变量之间显著相关
# 残差分析
res <- residuals(fit6.1)
# 接近正态分布
hist(res)
qqnorm(res)
qqline(res)
plot(fit6.1)
# 两个变量显著相关
6-2
在第 5 章习题 5-4 中我们探索了数据集 glucose 里各个变量之间的相关性,请使用这个数据集建立线性回归模型,研究血糖和其他几项指标的关系。
血清总胆固醇(x1,mmol/L)、甘油三酯(x2,mmol/L)、空腹胰岛素(x3,μU/mL)、糖化血红蛋白(x4,%)和空腹血糖的测量值(y,mmol/L)
# 构建数据
glucose <- data.frame(x1 = c(5.68, 3.79, 6.02, 4.85, 4.6, 6.05, 4.9, 7.08, 3.85, 4.65, 4.59, 4.29, 7.97, 6.19, 6.13, 5.71, 6.4, 6.06, 5.09, 6.13, 5.78, 5.43, 6.5, 7.98, 11.54, 5.84, 3.84),x2 = c(1.9, 1.64, 3.56, 1.07, 2.32, 0.64, 8.5, 3, 2.11, 0.63, 1.97, 1.97, 1.93, 1.18, 2.06, 1.78, 2.4, 3.67, 1.03, 1.71, 3.36, 1.13, 6.21, 7.92, 10.89, 0.92, 1.2),x3 = c(4.53, 7.32, 6.95, 5.88, 4.05, 1.42, 12.6, 6.75, 16.28, 6.59, 3.61, 6.61, 7.57, 1.42, 10.35, 8.53, 4.53, 12.79, 2.53, 5.28, 2.96, 4.31, 3.47, 3.37, 1.2, 8.61, 6.45),x4 = c(8.2, 6.9, 10.8, 8.3, 7.5, 13.6, 8.5, 11.5, 7.9, 7.1, 8.7, 7.8, 9.9, 6.9, 10.5, 8, 10.3, 7.1, 8.9, 9.9, 8, 11.3, 12.3, 9.8, 10.5, 6.4, 9.6),y = c(11.2, 8.8, 12.3, 11.6, 13.4, 18.3, 11.1, 12.1, 9.6, 8.4, 9.3, 10.6, 8.4, 9.6, 10.9, 10.1, 14.8, 9.1, 10.8, 10.2, 13.6, 14.9, 16, 13.2, 20, 13.3, 10.4)
)
str(glucose)
# 查看变量间的相关性,可以看到y和x1、x3、x4之间的相关指数都在0.5以上
cor(glucose)
# 构建模型
fit6.2 <- lm(y ~ x1 + x2 + x3 + x4, data = glucose)
# 回归分析
summary(fit6.2)
# 可知:决定系数0.6008,表明y变异的60%可由上述4个自变量的变化来解释。
# 自变量间的相关性系数
cor(glucose[,1:4])
# x1和x2相关系数0.6,相关性较强# 使用car::vif()函数计算模型中自变量对应的VIF值,一般VIF>5时,认为存在多重共线性
library(car)
vif(fit6.2)
# VIF<5,不存在多重共线性# 汇总模型主要结果
library(epiDisplay)
regress.display(fit6.2)#保存结果
write.csv(regress.display(fit6.2)$table,file='mytable.csv')
6-3
对于 ISwR 包里的数据集 cystfibr,以变量 pemax 作为因变量建立“最优”回归模型以研究患者肺功能的影响因素。
library(ISwR)
data("cystfibr")
str(cystfibr)
#将性别转化为因子型变量,避免建模时作为数字处理
cystfibr$sex <- factor(cystfibr$sex, labels = c("male", "female"))
cor(cystfibr[, 6:10])
# 后5个参数间相关系数很多都0.6以上,相关性中等以上;按题意取pemax做后5个参数的代表
# 构建模型
fit6.3 <- lm(pemax ~ age + sex + height + weight + bmp, data = cystfibr)
summary(fit6.3)
# weight和bmp的p<0.05,有统计学意义
# F检验的p值0.006<0.05,是显著的。决定系数0.548,表明pemax变异的54.8%可由上述5个自由量的变化来解释。
# 1-5个变量中有weight和bmp对pemax的影响有统计学意义,需要比较其他自由变量的线性关系,即多重共线性问题。
cor(cystfibr[, 3:5])
# 相关系数都高
# vif量化确认
library(car)
vif(fit6.3)
# age,height,weight的VIF>5,需要解决多重共线性,用逐步回归法剔除共线性因素
fit6.3.2 <- step(fit6.3)
summary(fit6.3.2)
# 剔除后只剩下weight+bmp
# fit2的决定系数0.4749与fit1的决定系数0.548很接近,fit2调整后的决定系数Adj R^2(0.427)也和fit1的0.429接近,说明在减少3个变量情况下,fit2的解释效率与fit1相当。
