zhang2sha的个人博客分享 http://blog.sciencenet.cn/u/zhang2sha

博文

回归模型中,拟合线置信区间的宽窄是如何计算的?

已有 1806 次阅读 2023-4-5 17:26 |系统分类:科研笔记

回归,无疑是最为常用的统计模型,它能告诉我们两个变量X和Y之间的定量关系,并让我们根据新的X值对Y的值进行预测。

比如,我们可以拟合一个最简单的线性回归模型:

data("cars")

m1 <- lm(dist~speed, data=cars)

我们可以在R中用ggplot2很方便的做出拟合图

library(ggplot2)

ggplot(cars, aes(x=speed, y=dist)) +

  geom_point(color='#2980B9', size = 4) +

  geom_smooth(method=lm, color='#2C3E50')

图片1.png

此图中的黑色拟合线,即X, Y之间的定量关系,我们都知道是根据模型给出的回归公式,即: Y=b0+b1X, 其中b0为截距,b1为斜率。其具体值通过summary(m1)我们就能知道。那么现在问题来了,这里面的拟合线的置信区间,也就是图中的灰色阴影部分,是根据什么计算的呢?

图片2.png

要弄清楚这个问题,首先要了解一个基本原理:

对于两个变量b0b1,如果b0的方差为V0, b1的方差为V1

1)假如b0b1是彼此完全独立的,则:

b1X对应的方差为V1X

b0 + b1X对应的方差为V0+V1X2


2)假如b0b1不是彼此完全独立的,即二者存在一定的相关性,且其协方差为Vcov, 那么这时:

b1X对应的方差=(Vcov+V1X)X

          b0 + b1X对应的方差=(V0+VcovX)+(Vcov+V1X)X

其中(V0+VcovX)为由于b0b1的不独立性及X的变化而造成的方差。

     理解了这个之后,我们计算拟合线的置信区间就变得非常简单了。

首先用vcov函数查看m1b0b1中的VCV矩阵:

图片3.png

其中矩阵的对角线位置为分别为截距b0和斜率b1的方差,即上面提到的V0和V1, 非对角线位置为二者的协方差,即上面提到的Vcov

那么对于给定的X, 模型对于Y的估计值Y=b0+b1X


Y的方差为:Vpredict=(V0+VcovX)+(Vcov+V1X)X

Y的标准误 SE=sqrt(Vpredict), 也就是Vpredict的平方根。

那么这时Y的上下置信区间为:95%CI=Y +/- 1.96*SE

理解了上述原理之后,我们便可以自定义一个小函数,求给定模型和X后,YSE,

se<-function(x,model){

 

  1*vcov(model)[1,1]+x*vcov(model)[1,2]

 

  1*vcov(model)[2,1]+x*vcov(model)[2,2]

 

  vi=vcov(model)[1,1]+x*vcov(model)[1,2]*1 + (1*vcov(model)[2,1]+x*vcov(model)[2,2])*x

  se=sqrt(vi)

  return(se)

 

}

我们可以检验下用该函数求得的SE和用predict得到的SE是否相同:

pred = as.data.frame(predict(m1, newdata =cars,level = 0.05,se.fit = T))

plot(pred$se.fit, se(cars$speed, m1))

结果如图:

\"图片4.png\"

可见,两种方法得到的SE 100%完全相同。然后,我们可以通过基础作图函数把拟合线的置信区间做出:

plot(cars$speed, cars$dist)

lines(x=cars$speed, y=cars$pred,col="red")

lines(x=cars$speed, y=cars$pred+1.96*se(cars$speed, m1),col="red")

lines(x=cars$speed, y=cars$pred-1.96*se(cars$speed,m1),col="red")

 图片5.png

同样,我们也可以通过ggplot2把我们手动计算的置信区间添加到图上

ggplot(cars, aes(x=speed, y=dist)) +

  geom_ribbon(aes(ymin = cars$pred-1.96*se(cars$speed,m1),

ymax = cars$pred+1.96*se(cars$speed, m1)), fill = "grey70") +

  geom_line(aes(y = cars$pred))+

  geom_point(color='#2980B9', size = 4)

 

图片6.png

跟用ggplot内置函数做出来的置信区间,分毫不差,完全一样。

最后如果你对矩阵运算比较熟悉的话,那SE的计算就更简单了:

matrix = model.matrix(m1)

se = sqrt(diag( matrix %*% vcov(m1) %*% t(matrix)))

得到的结果依旧分毫不差:

plot(pred$se.fit, se)

 

图片7.png

正所谓万变不离其宗,对于学习、应用统计工具的我们来说,任何时候,对于基础统计概念和原理的理解,都要放在第一位。最后欢迎大家关注本人这周末将要举办的混合效应模型培训班:

图片8.png




https://blog.sciencenet.cn/blog-3442043-1383080.html

上一篇:统计扫盲:回归模型中,解释变量可以是分类变量么?
下一篇:贝叶斯统计的独门优势---模型参数估计结果的二次运算
收藏 IP: 111.197.234.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-2-26 15:40

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部