回归，无疑是最为常用的统计模型，它能告诉我们两个变量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')

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

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

对于两个变量b₀和b₁，如果b₀的方差为V₀, b₁的方差为V₁，

1）假如b₀和b₁是彼此完全独立的，则：

b₁X对应的方差为V₁X

b₀ + b₁X对应的方差为V₀+V₁X²

2）假如b₀和b₁不是彼此完全独立的，即二者存在一定的相关性，且其协方差为V_cov, 那么这时：

b₁X对应的方差=(V_cov+V₁X)X

          b₀ + b₁X对应的方差=(V₀+VcovX)+(Vcov+V₁X)X

其中(V₀+VcovX)为由于b₀和b₁的不独立性及X的变化而造成的方差。

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

首先用vcov函数查看m1中b₀和b₁中的VCV矩阵：

其中矩阵的对角线位置为分别为截距b₀和斜率b₁的方差，即上面提到的V₀和V₁, 非对角线位置为二者的协方差，即上面提到的Vcov。

那么对于给定的X, 模型对于Y的估计值Y=b₀+b₁X

Y的方差为：Vpredict=(V₀+VcovX)+(Vcov+V₁X)X

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

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

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

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))

结果如图：

可见，两种方法得到的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")

同样，我们也可以通过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)

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

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

matrix = model.matrix(m1)

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

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

plot(pred$se.fit, se)

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