线性回归

原文地址：线性回归

4.1 一元线性模型

一元线性模型定义如下

$$y_i={\beta }_0+{\beta }_1{x}_i+{ϵ}_i,i=1,...,n$$

其中

• ${ϵ}_i$为随机变量（需要一些假设）
• $x_i$为固定值（独立/预测变量）
• $y_i$为随机变量（依赖/响应变量）
• ${\beta }_0$为截距
• ${\beta }_1$为斜率

4.1.1 最小二乘估计

选择参数${\beta }_0,{\beta }_1$最小化以下函数Q：

$$Q({\beta }_0,{\beta }_1)=\sum _{i=1}^n(y_i-{\beta }_0-{\beta }_1{x}_i)$$

函数取最小值时，参数$\widehat{{\beta }_0},\widehat{{\beta }_1}$满足：

$$\frac{\partial Q}{\partial {\beta }_0}=-2\sum _{i=1}^n(y_i-{\hat{\beta }}_0-{\hat{\beta }}_1{x}_i)=0$$

$$\frac{\partial Q}{\partial {\beta }_1}=-2\sum _{i=1}^n(y_i-{\hat{\beta }}_0-{\hat{\beta }}_1{x}_i)x_i=0$$

求解得出：

$$\widehat{{\beta }_1}=\frac{{\sum }_{i=1}^n(y_i-\overline{y})x_i}{{\sum }_{i=1}^n(x_i-\overline{x})x_i}$$

$${\beta }_0=\overline{y}-\widehat{{\beta }_1}\overline{x}$$

我们另外做出如下定义：

$${\zeta }_{xx}=\sum _{i=1}^n(x_i-\overline{x}{)}^2$$

$${\zeta }_{yy}\sum _{i=1}^n(y_i-\overline{y}{)}^2$$

$${\zeta }_{xy}\sum (x_i-\overline{x})(y_i-\overline{y})$$

则可以得出：

$$\widehat{{\beta }_1}=\frac{{\sum }_{i=1}^n(y_i-\overline{y})(x_i-\overline{x})}{{\sum }_{i=1}^n(x_i-\overline{x})(x_i-\overline{x})}=\frac{{\zeta }_{xy}}{{\zeta }_{xx}}=\frac{1}{{\zeta }_{xx}}\sum _{i=1}^n(x_i-\overline{x})y_i$$

由此得出的回归函数为：

$$\hat{y}=\widehat{{\beta }_0}+\widehat{{\beta }_1}x$$

4.1.2 期望和方差

假设 A1： $E[{ϵ}_i]=0,i=1,...,n$.

定理4.1： 假设A1成立，则$\widehat{{\beta }_0},\widehat{{\beta }_1}$是${\beta }_0,{\beta }_1$的无偏估计。

证明：

$$E[\widehat{{\beta }_1}]=\frac{1}{{\zeta }_{xx}}\sum _{i=1}^n(x_i-\overline{x})E[y_i]$$

$$=\frac{1}{{\zeta }_{xx}}\sum _{i=1}^n(x_i-\overline{x})({\beta }_0+{\beta }_1{x}_i)$$

$$=\frac{{\beta }_0}{{\zeta }_{xx}}\sum _{i=1}^n(x_i-\overline{x})+\frac{{\beta }_1}{{\zeta }_{xx}}\sum _{i=1}^n(x_i-\overline{x})x_i$$

$$=\frac{{\beta }_0}{{\zeta }_{xx}}(\sum _{i=1}^n{x}_i-\sum _{i=1}^n\overline{x})+\frac{{\beta }_1}{{\zeta }_{xx}}(x_i-\overline{x})x_i$$

$$=\frac{{\beta }_0}{{\zeta }_{xx}}(n\overline{x}-n\overline{x})+\frac{{\beta }_1}{{\zeta }_{xx}}\sum _{i=1}^n(x_i-\overline{x})x_i$$

$$=\frac{{\beta }_1}{{\zeta }_{xx}}\sum _{i=1}^n(x_i-\overline{x})x_i$$

$$=\frac{{\beta }_1}{{\zeta }_{xx}}\sum _{i=1}^n(x_i-\overline{x})x_i-\frac{\beta 1\overline{x}}{{\zeta }_{xx}}\sum _{i=1}^n(x_i-\overline{x})$$

$$=\frac{{\beta }_1}{{\zeta }_{xx}}\sum _{i=1}^n(x_i-\overline{x})(x_i-\overline{x})$$

$$={\beta }_1$$

$$E[\widehat{{\beta }_0}]=E[\overline{y}-\widehat{{\beta }_1}\overline{x}]=E[\overline{y}]-{\beta }_1\overline{x}={\beta }_0+{\beta }_1\overline{x}-{\beta }_1\overline{x}={\beta }_0$$

假设 A2： C o v ( ϵ i , ϵ j ) = σ 2 { i = j } , C o v ( ϵ i , ϵ j ) = 0 { i ≠ j } Cov(\epsilon_i,\epsilon_j)=\sigma^2\{i=j\},Cov(\epsilon_i,\epsilon_j)=0\{i\neq j\} Cov(ϵi​,ϵj​)=σ2{i=j},Cov(ϵi​,ϵj​)=0{i​=j}

定理4.2 : 假设A2成立，则有

$$Var[{\hat{\beta }}_0]=(\frac{1}{n}+\frac{{\overline{x}}^2}{{\zeta }_{xx}}){\sigma }^2$$

$$Var[\widehat{{\beta }_1}]=\frac{{\sigma }^2}{{\zeta }_{xx}}$$

$$Cov(\widehat{{\beta }_0},\widehat{{\beta }_1})=\frac{-\overline{x}}{{\zeta }_{xx}}{\sigma }^2$$

证明后补

4.1.3 误差项的方差估计

上文的A2假设中，方差${\sigma }^2$是一般是未知的，下面的理论给出${\sigma }^2$的无偏估计。

定义4.1： 我们将误差的平方和（SSE）记录为:

$$S_e^2=\sum _{i=1}^n(y_i-\widehat{{\beta }_0}-\widehat{{\beta }_1}{x}_i{)}^2$$

定理 4.3：首先我们定义

$${\hat{\sigma }}^2=\frac{Q({\hat{\beta }}_0,{\hat{\beta }}_1)}{n-2}=\frac{S_e^2}{n-2}$$

假设A1，A2成立，则有$E[{\hat{\sigma }}^2]={\sigma }^2$成立

证明后补

4.1.4 抽样分布定理

假设B：${ϵ}_i{\sim }^{iid}N(0,{\sigma }^2),i=1,...,n$

iid:independently identically distribution,表示独立同分布

定理4.4： 假设B成立，则有

1. ${\hat{\beta }}_0\sim N({\beta }_0,(\frac{1}{n}+\frac{\overline{x}}{{\zeta }_{xx}}{\sigma }^2)$
2. ${\hat{\beta }}_1\sim N({\beta }_1,\frac{{\sigma }^2}{\zeta xx})$
3. $\frac{(n-2){\hat{\sigma }}^2}{{\sigma }^2}=\frac{S_e^2}{{\sigma }^2}\sim {\chi }^2(n-2)$
4. ${\hat{\sigma }}^2$不依赖于$(\widehat{{\beta }_0},\widehat{{\beta }_1})$

4.1.5 置信区间与假设检验

给定$\sigma$的情况下，我们可以使用$\frac{\widehat{{\beta }_1}-{\beta }_1}{\sigma \sqrt{{\zeta }_{xx}}}\sim N(0,1)$确定检验和置信区间。${\beta }_1$的$100(1-\alpha )\%$置信区间为$\widehat{{\beta }_1}\pm u_{1-\alpha /2}\sigma /\sqrt{{\zeta }_{xx}}$

更一般情况是我们不知道$\sigma$只要满足$n\ge 3$，则有

$$\frac{\widehat{{\beta }_1}-{\beta }_1}{\hat{\sigma }/\sqrt{{\zeta }_{xx}}}\sim t(n-2)$$

$${\beta }_1$$的$100(1-\alpha )\%$置信区间为${\beta }_1\pm t_{1-\alpha /2}(n-2)\hat{\sigma }/\sqrt{{\zeta }_{xx}}$

${\beta }_0$的预测推到，可以使用：

$$\frac{\widehat{{\beta }_0}-{\beta }_0}{\sigma \sqrt{1/n+{\overline{x}}^2/{\zeta }_{xx}}}\sim N(0,1)$$

$$\frac{\widehat{{\beta }_0}-{\beta }_0}{\hat{\sigma }\sqrt{1/n+{\overline{x}}^2/{\zeta }_{xx}}}\sim t(n-2)$$

${\sigma }^2$的$100(1-\alpha )\%$的置信区间为:

$$[\frac{(n-2){\hat{\sigma }}^2}{{\chi }_{1-\alpha /2}^2(n-2)},\frac{(n-2){\hat{\sigma }}^2}{{\chi }_{\alpha /2}^2(n-2)}]=[\frac{S_e^2}{{\chi }_{1-\alpha /2}^2(n-2)},\frac{S_e^2}{{\chi }_{\alpha /2}^2(n-2)}]$$

4.1.6 案例分析1

一家空调商家存在空调组装问题，原因是连接杆无法满足成品重量规格。很多的连接杆在被最后加工后，由于超重问题而无法使用。为了降低成本，该商家的质量控制部门希望量化成品连接杆重量$y$与毛坯铸件重量$x$之间的关系，在进行最后加工（成本较高）之前，将超重的毛坯铸件剔除。

数据显示如下：

rod = data.frame(
  id = seq(1:25),
  rough_weight = c(2.745, 2.700, 2.690, 2.680, 2.675,                        
                   2.670, 2.665, 2.660, 2.655, 2.655,                        
                   2.650, 2.650, 2.645, 2.635, 2.630,
                   2.625, 2.625, 2.620, 2.615, 2.615,                        
                   2.615, 2.610, 2.590, 2.590, 2.565),
finished_weight = c(2.080, 2.045, 2.050, 2.005, 2.035,                        
                    2.035, 2.020, 2.005, 2.010, 2.000, 
                    2.000, 2.005, 2.015, 1.990, 1.990, 
                    1.995, 1.985, 1.970, 1.985, 1.990,                        
                    1.995, 1.990, 1.975, 1.995, 1.955)
)
knitr::kable(rod, caption = "rough weight vs. finished weight")

考虑采用线性模型拟合二者重量的关系：

$y_i={\beta }_0+{\beta }_1{x}_i+{ϵ}_i，{ϵ}_i{\sim }^{iid}N(0,{\sigma }^2)$

由观察数据得出$\overline{x}=2.643,\overline{y}=2.0048,{\zeta }_{xx}=0.0367,{\zeta }_{xy}=0.023565,\hat{\sigma }=0.0113$.

最小二乘估计为：

${\beta }_1=\frac{{\zeta }_{xy}}{{\zeta }_{xx}}=\frac{0.023565}{0.0367}=0.642,{\beta }_0=\overline{y}-\widehat{{\beta }_1}\overline{x}=0.308$

回归函数为$\hat{y}=0.308+0.642x$. 图像如下图蓝色直线：

4.1.7 拟合的评估

为了帮助评估拟合的质量，我们将广泛使用残差的概念，即观察值与拟合值之间的差：

$\widehat{{ϵ}_i}=y_i-\widehat{{\beta }_0}-\widehat{{\beta }_1}{x}_i,i=1,...,n.$

以图形方式检查残差最有用。残差与 x 值的关系可能揭示系统失配或数据不符合拟合模型的方式。理想情况下，残差应该与 x 值无关，绘图应看起来像是水平模糊。案例1的残差如下图所示。

par(mar=c(4,4,2,1))
plot(lm.rod$fitted.values,lm.rod$residuals,"p",
     xlab="Fitted values",ylab = "Residuals")

标准化残差如下图所示

par(mar=c(4,4,2,1))
plot(lm.rod$fitted.values,rstandard(lm.rod),"p",
     xlab="Fitted values",ylab = "Standardized Residuals")
abline(h=c(-2,2),lty=c(5,5))

4.1.8 预测

有关$E[y_{n+1}]$的预测

给定$x_{n+1}$，我们想要估计$y_{n+1}$的期望值，比如，$E[y_{n+1}]={\beta }_0+{\beta }_1{x}_{n+1}$。一个自然的无偏估计为${\hat{y}}_{n+1}=\widehat{{\beta }_0}+\widehat{{\beta }_1}{x}_{n+1}$. 由上文的定理4.3，我们知道：

$$Var[{\hat{y}}_{n+1}]=(\frac{1}{n}+\frac{(x_{n+1}-\overline{x}{)}^2}{{\zeta }_{xx}}){\sigma }^2$$

假设B成立，由定理4.4得出：

$${\hat{y}}_{n+1}\sim N({\beta }_0+{\beta }_1{x}_{n+1},(1/n+(x_{n+1}-\overline{x}{)}^2)/{\sigma }^2)$$

$$\frac{{\hat{y}}_{n+1}-E[y_{n+1}]}{\hat{\sigma }\sqrt{1/n+(x_{n+1}-\overline{x}{)}^2/{\zeta }_{xx}}}\sim t(n-2)$$

因此得出下列结论

定理4.5： 假设B成立，则有

$${\hat{y}}_{n+1}=\widehat{{\beta }_0}+\widehat{{\beta }_1}{x}_{n+1}\sim N({\beta }_0+{\beta }_1{x}_{n+1},[1/n+(x_{n+1}-\overline{x}{)}^2/{\zeta }_{xx}]{\sigma }^2)$$。

$E[y_{n+1}]={\beta }_0+{\beta }_1{x}_{n+1}$的100(1-\alpha)%$置信区间为

$${\hat{y}}_{n+1}\pm t_{1-\alpha /2}(n-2)\hat{\sigma }\sqrt{\frac{1}{n}+\frac{(x_{n+1}-\overline{x}{)}^2}{{\zeta }_{xx}}}$$

有关未来观察值$y_{n+1}$的预测

现在我们来给出一个未来观察值$y_{n+1}$的预测区间，而不是期望值$E[y_{n+1}]$的区间。需要注意的是$y_{n+1}$不再是一个固定的参数。一个预测区间指的是一个以指定概率包含$y_{n+1}$的区间。我们考虑 $y_{n+1}-{\hat{y}}_{n+1}$，如果A1假设成立（即$$E[{ϵ}_i]=0,i=1,..,n$$），则有

$$E[y_{n+1}-{\hat{y}}_{n+1}]=E[y_{n+1}]-E[{\hat{y}}_{n+1}]$$

$$=E[y_{n+1}]-E[{\hat{\beta }}_0+{\hat{\beta }}_1{x}_{n+1}]$$

$$=E[y_{n+1}]-[E[\widehat{{\beta }_0}]+E[{\hat{\beta }}_1{x}_{n+1}]]$$

$$=E[y_{n+1}]-[{\beta }_0+{\beta }_1{x}_{n+1}]$$

$$=E[{\beta }_0+{\beta }_1{x}_{n+1}+{ϵ}_{n+1}]-[{\beta }_0+{\beta }_1{x}_{n+1}]$$

$$={\beta }_0+{\beta }_1{x}_{n+1}+E[{ϵ}_{n+1}]-{\beta }_0-{\beta }_1{x}_{n+1}$$

$$=0$$

如果A2假设成立（即 C o v ( ϵ i , ϵ j ) = σ 2 { i = j } , = 0 { i ≠ j } Cov(\epsilon_i,\epsilon_j)=\sigma^2\{i=j\},=0\{i\ne j\} Cov(ϵi​,ϵj​)=σ2{i=j},=0{i​=j}，则有

$$Var[y_{n+1}-{\hat{y}}_{n+1}]=Var[y_{n+1}]+Var[{\hat{y}}_{n+1}]=(1+\frac{1}{n}+\frac{(x_{n+1}-\overline{x}{)}^2}{{\zeta }_{xx}}){\sigma }^2$$

如果假设B成立，则$y_{n+1}-{\hat{y}}_{n+1}$满足正太分布。

定理4.6: 假设B成立，定义$y_{n+1}={\beta }_0+{\beta }_1{x}_{n+1}+ϵ$，其中$ϵ\sim N(0,{\sigma }^2)$

则y的$100(1-\alpha )\%$的置信区间为

${\hat{y}}_{n+1}\pm t_{1-\alpha /2}(n-2)\hat{\sigma }\sqrt{1+\frac{1}{n}+\frac{(x_{n+1}-\overline{x}{)}^2}{{\zeta }_{xx}}}$

4.1.9 控制

如何控制$y_{n+1}$的取值范围呢?

再次考虑案例1。在进行最终（且昂贵）的加工过程之前，可以丢弃可能产生过重杆的铸件。公司的质量控制部门想以不低于0.95的概率生产重量不超过2.05的杆，如果选择铸件的重量。

现在我们想要$y_{n+1}\le y_0=2.05$以$1-\alpha$成立。类似于定理4.6，我们可以为$y_{n+1}$构造单边置信区间：

$(-\infty ,{\hat{y}}_{n+1}+t_{1-\alpha }(n-2)\hat{\sigma }\sqrt{1+\frac{1}{n}+\frac{(x_{n+1}-\overline{x}{)}^2}{{\zeta }_{xx}}}]$

得出:

$${\hat{\beta }}_0+{\hat{\beta }}_1{x}_{n+1}+t_{1-\alpha }(n-2)\hat{\sigma }\sqrt{1+\frac{1}{n}+\frac{(x_{n+1}-\overline{x}{)}^2}{{\zeta }_{xx}}}\le y_0$$