应用机器学习（八）：线性模型

回归分析

在统计建模中，回归分析( Regression analysis ) 是用来刻画变量之间的统计关系的一种统计技术。当一个变量( 因变量 dependent variable )受另外一些变量( 自变量 independent variables or predictors )的强烈影响时，适宜用回归的方法。例如，”回归”的思想和方法最早由英国著名的统计学家F. Galton 和他的学生 K. Pearson 提出。他们在研究父母身高与其子女身高的遗传问题时，以每对夫妇的平均身高作为自变量 x, 而取他们的一个成年孩子的身高作为 y, 共观察了1,078对夫妇，建立了回归方程 y=33.73+0.516x, 从此产生了回归分析的方法。

回归模型

回归模型的一般形式

设因变量 yy, 自变量向量 x=(x1,x2,…,xp)\mathbf{x}=(x_1, x_2, \dots, x_p), 则刻画 yy 与 x\mathbf{x}关系的回归模型的一般形式为

y=f(x)+ε(1)

\begin{align*} y=f(\mathbf{x})+\varepsilon \tag{1} \end{align*}

其中， ε\varepsilon 为随机误差，它表示除了 x\mathbf{x} 外的其它随机干扰因素。正是因为 ε\varepsilon 的存在，给定 x\mathbf{x} 的值，不能唯一确定 yy, 称这种变量关系为统计关系。随机变量 ε\varepsilon 满足基本假设 E(ε)=0E(\varepsilon)=0, 即零均值，它表示没有系统误差。

线性回归

线性回归( Linear regression ), 指的是 yy 与 x\mathbf{x} 之间是线性关系，即

y=β0+β1x1+β2x2+⋯+βpxp+ε(2)

\begin{align*} y=\beta_0+\beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p +\varepsilon \tag{2} \end{align*}
其中的参数 β0,β1,…,βp\beta_0, \beta_1, \dots, \beta_p 称为回归系数。特别地， p=1p=1, 即一个自变量时，称模型为一元线性回归。

假设有 nn 次观测数据 (yi;xi1,xi2,…,xip),i=1,2,…,n( y_i; x_{i1}, x_{i2}, \dots, x_{ip} ),\, i=1,2,\dots, n.
记y=(y1,y2,…,yn)′\mathbf{y}=(y_1, y_2, \dots, y_n)', β=(β0,β1,…,βp)′\mathbf{\beta}=(\beta_0, \beta_1, \dots, \beta_p)',
ε=(ϵ1,ϵ2,…,ϵp)′\mathbf{\varepsilon}=(\epsilon_1, \epsilon_2, \dots, \epsilon_p)',

X=⎛⎝⎜⎜⎜⎜⎜11⋮1x11x21⋮xn1x12x22⋮xn2………x1px2p⋮xnp⎞⎠⎟⎟⎟⎟⎟

\begin{equation*} \mathbf{X}=\left( \begin{array}{ccccc} 1 & x_{11} & x_{12} & \dots & x_{1p}\\ 1 & x_{21} & x_{22} & \dots & x_{2p}\\ \vdots & \vdots & \vdots & & \vdots\\1 & x_{n1} & x_{n2} & \dots & x_{np}\\\end{array} \right) \end{equation*}

则观测数据的线性回归模型可以表示为

yi=β0+β1xi1+β2xi2+⋯+βpxip+εi,i=1,2,…,n(3)

\begin{align*} y_i=\beta_0+\beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip} +\varepsilon_i, \, i=1,2,\dots,n \tag{3} \end{align*}
等价的矩阵形式为

y=Xβ+ε(4)

\begin{align*} \mathbf{y}=\mathbf{X}\beta +\varepsilon \tag{4} \end{align*}

线性回归模型的基本假设

自变量 x\mathbf{x} 是确定性变量，不是随机变量，因变量 yy 是随机变量。要求
设计矩阵 X\mathbf{X} 是列满秩的，即 rank(X)=p+1<nrank(\mathbf{X})=p+1, 这表明 X\mathbf{X}
的列向量之间线性无关;
随机误差零均值同方差，即 E(ε)=0,cov(ε,ε)=σ2EnE(\varepsilon)=0,\, cov(\varepsilon, \varepsilon)=\sigma^2 E_n.
这个假设常称为高斯-马尔科夫条件( Gauss-Markov assumption );
通常，假定随机误差向量 ε\varepsilon 服从正态分布，即 ε∼N(0,σ2En)\varepsilon\sim N(0, \sigma^2 E_n),
因此，有 y∼N(Xβ,σ2En)\mathbf{y}\sim N(\mathbf{X}\beta, \sigma^2 E_n).

回归系数的估计

普通最小二乘估计

普通最小二乘估计( ordinary least square estimation, OLSE )的思想是，考虑 nn 次观测的 yiy_i 与均值 E(yi)=β0+β1xi1+β2xi2+⋯+βpxipE(y_i)=\beta_0+\beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip}
的离差平方和

Q(β)=∑i=1n(yi−E(yi))2=∑i=1n(yi−β0−β1xi1−β2xi2−⋯−βpxip)2(5)

\begin{align*} Q(\beta)=\sum\limits_{i=1}^n (y_i - E(y_i))^2=\sum\limits_{i=1}^n (y_i - \beta_0-\beta_1 x_{i1} - \beta_2 x_{i2} - \dots - \beta_p x_{ip})^2 \tag{5} \end{align*}

所谓最小二乘估计，就是找回归系数β^=(β0^,β1^,…,βp^)′\hat{\beta}=(\hat{\beta_0}, \hat{\beta_1}, \dots, \hat{\beta_p})',
使得

Q(β^)=minβQ(β)(6)

\begin{align*} Q(\hat{\beta})=\mathop{\min}_{\beta}Q(\beta) \tag{6} \end{align*}

注意到，函数 QQ 关于 β0,β1,…,βp\beta_0, \beta_1, \dots, \beta_p 是非负可导的二次函数，故存在最小值，且在导数为0的点处。这样，解得 β\beta 最小二乘估计为

β^=(X′X)−1X′y(7)

\begin{align*} \hat{\beta}=(\mathbf{X}'\mathbf{X})^{-1} \mathbf{X}' \mathbf{y} \tag{7} \end{align*}
y\mathbf{y} 的回归值向量

y^=Xβ^(8)

\begin{align*} \hat{\mathbf{y}}=\mathbf{X}\hat{\beta} \tag{8} \end{align*}

最大似然估计

最大似然估计( maximum likelihood estimation ), 是在正态假设下，y∼N(Xβ,σ2En)\mathbf{y}\sim N(\mathbf{X}\beta, \sigma^2 E_n). 因此， y\mathbf{y} 的似然函数为

L(β,σ2)=(2π)−n2(σ2)−n2exp(−12σ2(y−Xβ)′(y−Xβ))(9)

\begin{align*} L(\beta, \sigma^2)=(2\pi)^{-\frac{n}{2}}(\sigma^2)^{-\frac{n}{2}}exp(-\dfrac{1}{2\sigma^2} (\mathbf{y}-\mathbf{X}\beta)'(\mathbf{y}-\mathbf{X}\beta)) \tag{9} \end{align*}
等价地转换为对数似然

lnL(β,σ2)=−n2ln2π−n2lnσ2−12σ2(y−Xβ)′(y−Xβ)(10)

\begin{align*} \ln L(\beta, \sigma^2)=-\dfrac{n}{2}\ln 2\pi -\dfrac{n}{2}\ln \sigma^2 -\dfrac{1}{2\sigma^2} (\mathbf{y}-\mathbf{X}\beta)'(\mathbf{y}-\mathbf{X}\beta) \tag{10} \end{align*}
注意，这里的(对数)似然函数为 β,σ2\beta, \,\sigma^2 的函数，固定 σ2\sigma^2, 求 lnL\ln L 的最大值，等价于求 (y−Xβ)′(y−Xβ)(\mathbf{y}-\mathbf{X}\beta)'(\mathbf{y}-\mathbf{X}\beta) 的最小值，而这时离差平方和 QQ, 这样，回归系数向量 β\beta 的最大似然估计与最小二乘估计是等价的，下面，固定 β^\hat{\beta}, 得到误差项方差 σ2\sigma^2 为

σ^2=1n(y−Xβ^)′(y−Xβ^)(11)

\begin{align*} \hat{\sigma}^2=\dfrac{1}{n}(\mathbf{y}-\mathbf{X}\hat{\beta})'(\mathbf{y}-\mathbf{X}\hat{\beta}) \tag{11} \end{align*}
定义残差向量

e=y−y^=y−Xβ^

\begin{align*} \mathbf{e}=\mathbf{y}-\hat{\mathbf{y}}=\mathbf{y}-\mathbf{X}\hat{\beta} \end{align*}

残差 e\mathbf{e} 可以看作是随机误差 ε\varepsilon 的估计。那么， σ2\sigma^2 的最大似然估计还可以记为
σ^2=1ne′e\hat{\sigma}^2=\dfrac{1}{n}\mathbf{e}'\mathbf{e}.

回归方程的显著性检验

确定回归系数的估计后，在使用回归方程作预测或估计前，还需要对回归方程进行显著性检验，即，检验自变量和因变量之间是否有线性关系。下面介绍两种统计检验方法，一是回归方程的 F 检验，另一个是回归系数的 t 检验。同时介绍度量回归拟合数据程度的拟合优度。

回归方程的 F-检验

考虑所有自变量 x1,x2,…,xpx_1, x_2, \dots, x_p 作为一个整体是否与 yy 有线性关系，为此提出假设

H0:β=0H1:β≠0

\begin{align*} H_0: \beta=0 \qquad H_1: \beta \ne 0 \end{align*}

定义平方和

SST=∑i=1n(yi−y¯)2SSR=∑i=1n(y^i−y¯)2SSE=∑i=1n(yi−y^i)2

\begin{align*} & SST=\sum\limits_{i=1}^n (y_i - \bar{y})^2\\ & SSR=\sum\limits_{i=1}^n (\hat{y}_i -\bar{y})^2\\ & SSE=\sum\limits_{i=1}^n (y_i - \hat{y}_i)^2\\ \end{align*}

SSTSST 称离差平方和，它表示因变量的波动，即不确定性；
SSRSSR 称回归平方和，它表示回归方程的波动；
SSESSE 称残差平方和，它表示随机波动。
可以证明平方和分解式

SST=SSR+SSE(12)

\begin{align*} SST=SSR+SSE \tag{12} \end{align*}

令检验统计量

F=SSR/pSSE/n−p−1(13)

\begin{align*} F=\dfrac{SSR/p}{SSE/n-p-1} \tag{13} \end{align*}
H0H_0 为真时， F∼F(p,n−p−1)F\sim F(p, n-p-1), 因此，给定水平 α\alpha, 拒绝域的形式 W={F>Fα(p,n−p−1)}W=\{F> F_{\alpha}(p, n-p-1)\}. 当 F>Fα(p,n−p−1)F> F_{\alpha}(p, n-p-1) 时，拒绝 H0H_0, yy 与 x1,x2,…,xpx_1, x_2, \dots, x_p 有线性关系，即，回归方程是显著的；当 F≤Fα(p,n−p−1)F\le F_{\alpha}(p, n-p-1)时，则认为没有线性关系，回归方程不显著。

回归系数的t-检验

回归方程显著，并不意味着每个自变量都与因变量有线性关系。因此，还需要对每个自变量是否与因变量有线性关系作显著性检验。为此提出假设

H0j:βj=0H1j:βj≠0,j=1,2,…,p.

\begin{align*} H_{0j}: \beta_j = 0 \qquad H_{1j}: \beta_j \ne 0,\qquad j=1,2,\dots,p. \end{align*}

在正态假设下，可以证明 β^∼N(β,σ2(X′X)−1)\hat{\beta}\sim N(\beta, \sigma^2 (\mathbf{X}'\mathbf{X})^{-1}), 记 (X′X)−1=(cij),i,j=0,1,2,…,p(\mathbf{X}'\mathbf{X})^{-1}=(c_{ij}),\, i,j=0,1,2,\dots,p. 构造检验统计量 tj=β^jcjj−−√σ^t_j = \dfrac{\hat{\beta}_j}{\sqrt{c_{jj}}\hat{\sigma}},
σ^=1n−p−1SSE−−−−−−−−−−−−√\hat{\sigma}=\sqrt{\dfrac{1}{n-p-1}SSE}. 在 H0jH_{0j} 下， tj∼t(n−p−1)t_j\sim t(n-p-1), 因此，在水平 α\alpha 下，拒绝域
W={|tj|≥tα2}W=\{|t_j|\ge t_{\frac{\alpha}{2}}\}. 当 |tj|≥tα2|t_j|\ge t_{\frac{\alpha}{2}} 时，拒绝 H0jH_{0j}, 认为自变量 xjx_j 与 yy 有线性关系。否则，可以认为没有线性关系。

拟合优度

拟合优度( Goodness of fit ), 用于考察回归方程对观测数据的拟合程度，在回归中称决定系数，定义为

R2=SSRSST=1−SSESST(14)

\begin{align*} R^2=\dfrac{SSR}{SST}=1-\dfrac{SSE}{SST} \tag{14} \end{align*}

R2R^2 在区间 [0,1][0, 1] 内取值，R2R^2 越接近1，表明回归拟合的效果越好；相反，越接近0，表明拟合的越差。与 F-检验比， R2R^2 可以更清楚地反映回归拟合的效果，但不能作为严格的显著性检验。可以证明，当观测数 n=2 时，不论变量间有没有相关性， R2=1R^2=1.

逐步回归法

在建立一个实际问题的回归模型时，首先遇到的问题是如何确定自变量。如果遗漏了某些重要的自变量，回归效果不好；如果包括过多的自变量，这些变量可能产生共线性，回归模型的稳定性较差。首先，介绍自变量选择的重要准则，赤池信息量准则。

自变量的选择

赤池信息量准则( Akaike information criterion, AIC ), 是根据最大似然估计的原理提出的一种模型选择准则。对一般情况，设模型的似然函数为 L(θ,x)L(\mathbf{\theta}, \mathbf{x}), θ\theta 的维数为 pp, x\mathbf{x} 为随机样本，则 AIC 定义为

AIC=−2lnL(θ^,x)+2p(15)

\begin{align*} AIC=-2\ln L(\hat{\theta}, \mathbf{x}) + 2 p \tag{15} \end{align*}

其中， θ^\hat{\theta} 是 θ\theta 的最大似然估计。似然函数越大的估计量越好， AICAIC 是似然函数的对数乘 -2 再加上惩罚因子 2p2p , 因而选择使 AICAIC 达到最小的模型是”最优”模型。
下面讨论把 AICAIC 用于回归模型的选择。假定回归模型的随机误差 ε∼N(0,σ2)\varepsilon\sim N(0, \sigma^2),
那么，最大对数似然

lnLmax=−n2ln2π−n2lnσ^2−12σ^2SSE(16)

\begin{align*} \ln L_{max} =-\dfrac{n}{2}\ln 2\pi -\dfrac{n}{2}\ln \hat{\sigma}^2 - \dfrac{1}{2} \hat{\sigma}^2 SSE \tag{16} \end{align*}

将 σ^2=1nSSE\hat{\sigma}^2 = \dfrac{1}{n} SSE 代入，得

lnLmax=−n2ln2π−n2lnSSEn−n2(17)

\begin{align*} \ln L_{max} =-\dfrac{n}{2}\ln 2\pi -\dfrac{n}{2}\ln \dfrac{SSE}{n} - \dfrac{n}{2} \tag{17} \end{align*}

代入式(15)，得

AIC=nlnSSE+2p(18)

\begin{align*} AIC = n\ln SSE + 2p \tag{18} \end{align*}

在回归建模过程中，对每一个回归的自变量子集，计算 AICAIC, 其中的 AICAIC 最小的模型为”最优”的回归模型。

逐步回归法

逐步回归的基本思想是”有进有出”。在自变量的选择过程中，将变量一个一个引入，每引入一个自变量后，对已选入的变量进行逐个显著性检验。当原先引入的变量由于后面变量的引入而变得不再显著时，要将其从回归模型中剔除。引入一个变量或从回归模型中剔除一个变量，为逐步回归的一步，每一步都要进行 F 检验，以确保每次引入新的变量之前回归模型中只包含显著的变量。这个过程反复进行，直到既无显著的自变量选入，也无不显著的自变量剔除为止。这样，保证了最后所得的回归模型是最优的。在逐步回归过程，要求引入自变量的检验水平 αentry\alpha_{entry}, 小于剔除自变量的检验水平 αremoval\alpha_{removal}, 否则，可能产生”死循环”。

一个实例：带交互项的线性回归

我们以 R 包 datasets 的数据集 mtcars 为例，以 mpg (每加仑汽油行驶的英里数)为因变量，
hp (马力)、 wt (汽车重量)为自变量，使用函数 lm() 建立线性回归模型。

fit <- lm(mpg ~ hp + wt, data = mtcars)
summary(fit)

从 Pr(>|t|) 栏可以看到，马力和车重都是显著的，回归拟合优度 R-squared 和 adjusted R-squared 都超过了0.8，说明线性程度较高。进一步考虑，马力和车重是否存在交互效应，为此，建立带交互项的回归模型。

fit <- lm(mpg ~ hp + wt + hp:wt, data = mtcars)
summary(fit)

从 Pr(>|t|) 栏可以看到，马力和车重的交互项是显著的，这说明响应变量 mpg 与其中一个变量的关系依赖于另一个变量的值。即，每加仑汽油行驶的英里数与汽车马力的关系依车重不同而不同。而且，通过与不带交互项的模型对比发现，带交互项的回归模型，拟合优度 R-squared 有所提高，说明带交互项的回归模型拟合数据更佳。预测模型为 mpg = 49.81 - 0.12hp - 8.22wt + 0.03hp·wt

使用 effects 包中的 effect() 函数，你可以用图形展示交互项的结果。

library(effects)
plot(effect("hp:wt", fit,xlevels=list(wt=c(2.2, 3.2, 4.2))),multiline = TRUE)

从结果图可以很清晰地看出，随着车重的增加（2.2、3.2、4.2），马力与每加仑汽油行驶英里数的关系减弱了。特别地，当 wt = 4.2 时，直线几乎是水平的，表明随着 hp 的增加， mpg 不会发生改变。

精彩内容，请关注微信公众号“统计学习与大数据”!