统计学习：线性回归

应用场景

线性回归是研究其它一些统计学习方法的起点，其它方法可看作线性回归的扩展。

在Advertising数据集，某种商品销售额（sales）可视为电视（TV）、广播（radio）和报纸（newspaper）广告预算的函数。通过这些数据，能否解决如下问题：

1. 广告预算和销量有关吗？
2. 如果相关，这种相关性有多强？也就是给定广告预算，能否准确预测销售额？
3. 各种类型的广告对销售额都有影响么吗？
4. 如果准确评估每种广告对销售额的影响？
5. 如何准确预测未来销售额？
6. 广告和销售额是线性关系吗？
7. 广告媒体之间有协同效应（synergy effect）吗？协同效应在统计学中称为相关性（interaction effect）。例如：在电视和广播分别投放\$50000的广告和将这\$100000投放到其中一种媒体，效果一样吗？

单变量回归分析

数学上，回归问题可表示为 \begin{equation} Y\approx \beta_0 + \beta_1 X， \end{equation} $X$可以指代电视，$Y$表示销售额。$\beta_0 $和$\beta_1$分别表示截距（intercept）和斜率（slope）。模型参数用训练数据上的最小二乘法估计可得 \begin{equation} \begin{aligned} \hat\beta_1&={\sum_{i=1}^n(x_i-\bar x)(y_i-\bar y)\over\sum_{i=1}^n(x_i-\bar x)^2}\\ \hat\beta_0&=\bar y-\hat\beta_1\bar x， \end{aligned} \label{eq:lsm-coefficients} \end{equation} 销售额可以预测 \[ \hat y=\hat\beta_0+\hat\beta_1x， \] 符号$\hat{\;}$表示估计值。第$i$个样本的残差（residual）记为$e_i=y_i-\hat y_i$，$n$个样本的残差平方和（RSS，residual sum of squares）为 \begin{equation} RSS=e_i^2+e_2^2+\ldots+e_n^2， \label{eq:rss} \end{equation} RSS也称为SSE。

假设$X$和$Y$的真实关系为 \begin{equation} Y=\beta_0+\beta_1X+\epsilon， \label{eq:population-regression-line} \end{equation} 截断$\beta_0$表示当$X=0$是$Y$的取值，斜率$\beta_1$表示$X$增加一个单位$Y$的平均增加量，误差$\epsilon$代表了这个简单模型缺失的其它所有量。$X$和$Y$的真实关系可能不是线性的，还可能有其它量导致$Y$的改变，也可能有测量误差。通常假定误差项与$X$独立。

模型\eqref{eq:population-regression-line}定义了总体回归直线（population regression line），它是$X$和$Y$真实关系的最佳线性近似。最小二乘法估计的系数\eqref{eq:lsm-coefficients}确定了最小二乘直线（least squares line）。

图 1: 线性回归模型 [PNG]

上图左中，红色直线表示真实关系$f(X)=2+3X$，也就是总体回归线，深蓝色直线表示最小二乘直线，它是基于数据的最小二乘估计；上图右的浅蓝色直线是10次最小二乘法估计的结果，每次估计采用的不同数据集，但是10次估计的平均值很接近总体回归线。

每次实验的100组数据通过模型$Y=2+3X+\epsilon$随机生成，$\epsilon$由零均值的正态分布产生。在实际应用中，只能通过数据得到最小二乘直线，无法观察到总体回归直线。

在一个数据集上通过\eqref{eq:lsm-coefficients}无法得到参数的准确估计，但是通过在大量数据集上估计结果的平均，最小二乘直线就是总体回归直线的无偏（unbiased）估计。

参数估计估计的精度可以通过标准误差（SE，standard error）衡量，

\begin{equation} \begin{aligned} SE\left(\hat\beta_0\right)^2&=\sigma^2\left[{1\over n}+{\bar x^2\over\sum_{i=1}^n(x_i-\bar x)^2}\right]\\ SE\left(\hat\beta_1\right)^2&={\sigma^2\over\sum_{i=1}^n(x_i-\bar x)^2}， \end{aligned} \end{equation}

$\sigma$表示均方差（standard deviation），$\sigma^2=Var(\epsilon)$。

标准误差揭示了参数估计值和真实值之间的平均差异。公式严格成立的条件是每个样本的误差$\epsilon_i$和公共方差（common variance）无关。即使这个条件不成立，这个公式也可作为评估的很好近似。当$x_i$越分散时，$SE\left(\hat\beta_1\right)^2$越小；当$\bar x=0$时，$SE\left(\hat\beta_0\right)^2$也更小，此时$\hat\beta_0=\bar y$；随着样本数$n$的增加，$SE\left(\hat\beta_0\right)^2$也越来越小。通常$\sigma^2$未知，但可以通过数据集上的标准化残差（RSE，residual standard error）估计，

\begin{equation} RSE=\sqrt{RSS\over n-2}=\sqrt{\frac{1}{n-2}\sum_{i=1}^n\left(y_i-\hat y_i\right)^2}。 \label{eq:rse} \end{equation}

标准误差可以用于计算置信区间（confidence interval）。$95\%$的置信区间表示该区间有$95\%$的概率包括参数的真实值。类似的均方根误差（RMSE，root-mean-square error）则定义为

\[ RMSE=\sqrt{RSS\over n}。 \]

对于线性回归的参数估计，$95\%$的置信区间记为¹：

\begin{equation} \hat\beta_0\pm 2\cdot SE\left(\hat\beta_0\right)，\qquad \hat\beta_1\pm 2\cdot SE\left(\hat\beta_1\right)。 \label{eq:lr-confidence-interval} \end{equation}

标准误差也可用于系数的假设检验（hypothesis test）。最常用的两种假设检验是无效假设检验（null hypothesis）和备选假设检验（alternative hypothesis）， \[ \begin{aligned} H_0&:\mbox{There is no relationship between }X\mbox{ and }Y\\ H_a&:\mbox{There is some relationship between }X\mbox{ and }Y， \end{aligned} \] 数学上表示为 \[ H_0:\beta_1=0，\qquad H_a:\beta_1\neq0。 \] 若$\beta_1=0$，那么$X$和$Y$无关。进行无效假设检验时，需要估计$\hat\beta_1$是否离0足够远，这需要通过$\hat\beta_1$的精度$SE\left(\hat\beta_1\right)$来衡量。如果$SE\left(\hat\beta_1\right)$很小，较小的$\hat\beta_1$可以足够确信$\beta_1\neq 0$；如果$SE\left(\hat\beta_1\right)$很大，$\hat\beta_1$的绝对值要很大才能拒绝无效假设。在实际应用中，通过$t$统计量（$t$-statistic）

\begin{equation} t={\hat\beta_1-0\over SE\left(\hat\beta_1\right)} \label{eq:t-statistic} \end{equation}

测量$\hat\beta_1$远离$0$（假设$H_0$成立）的标准差（standard deviation）。如果$X$和$Y$真的无关，\eqref{eq:t-statistic}是$n-2$自由度的$t$分布（$t$-distribution）。当$n$大约大于$30$时，$t$分布与正态分布很相似。$p$值（$p$-value）表示$\beta_1=0$时观测值大于等于$|t|$的概率。较小的$p$值表示预测（predictor）与响应（response）之间有关联。若$p$值很小时，拒绝无效假设（reject the null hypothesis）表示$X$与$Y$有关系。拒绝无效假设的典型$p$值是$5\%$或$1\%$，当$n=30$时，相应的$t$统计量分别为$2$和$2.75$。²

利用R的回归模型，可以得到如下的结果：

advertising <- read.csv('Advertising.csv')
lm.advertising <- lm(Sales ~ TV, data = advertising)
summary(lm.advertising)

# 结果如下：
Call:
lm(formula = Sales ~ TV, data = advertising)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.3860 -1.9545 -0.1913  2.0671  7.2124 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 7.032594   0.457843   15.36   <2e-16 ***
TV          0.047537   0.002691   17.67   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.259 on 198 degrees of freedom
Multiple R-squared:  0.6119,    Adjusted R-squared:  0.6099 
F-statistic: 312.1 on 1 and 198 DF,  p-value: < 2.2e-16

Coefficients的最后一列$Pr(>|t|)$表示$H_0$成立的概率，若很小就拒绝$H_0$接受$H_1$。拒绝$H_0$以为着$\beta_1\neq0$，表明$X$与$Y$相关。星号表示相关性强弱。最后一列表明销售额与电视广告有关系。$RSE$衡量了模型拟合结果与数据真实表现的差异。

$RSE$衡量了模型拟合结果与数据真实表现的差异。$R^2$统计量是度量拟合的另一种方式³， \begin{equation} R^2 = {TSS-RSS\over TSS}=1-{RSS\over TSS}，\qquad TSS=\sum_{i=1}^n\left(y_i-\bar y\right)^2， \label{eq:r2-statistic} \end{equation} TSS也记为SST。$R^2$越大越好。加入新的变量会使得$R^2$增大⁴，在测试集上$R^2$可能为负数，$R^2\in(-\infty, 1]$。若$R^2$接近$0$，这表明回归模型没有很好的解释响应变化，可能是线性模型有误或内在误差$\sigma^2$很高。本例中，$R^2\approx 0.61$表示少于三分之二的销售额变化能通过基于电视广告的回归模型解释。虽然$R^2\in[0,1]$，要比较哪个$R^2$更好还要根据具体应用而定。

$R^2$是衡量$X$与$Y$线性相关的统计量。相关系数（correlation） \begin{equation} r=Cor(X,Y)={ \sum_{i=1}^n\left(x_i-\bar x\right)\left(y_i-\bar y\right) \over \sqrt{\sum_{i=1}^n\left(x_i-\bar x\right)^2}\sqrt{\sum_{i=1}^n\left(y_i-\bar y\right)^2} } \label{eq:correlation} \end{equation} 也是衡量$X$与$Y$线性相关的统计量。事实上，$R^2=r^2$。

多变量回归分析

$p$个不同预测变量（predictor）的多变量回归模型为

\begin{equation} Y=\beta_0+\beta_1X_1+\beta_2X_2+\cdots+\beta_pX_p+\epsilon。 \label{eq:population-regression-line2} \end{equation}

这些预测变量对响应的贡献如何呢？

无效假设和备选假设分别为 \[ \begin{aligned} H_0&:\beta_1=\beta_2=\cdots=\beta_p=0\\ H_a&:\mbox{at least one }\beta_j\mbox{ is non-zero}。 \end{aligned} \]

假设检验利用$F$统计量 \begin{equation} F={(TSS-RSS)/p\over RSS/(n-p-1)}。 \label{eq:f-statistic} \end{equation}

若线性模型的假设成立，那么有

\[ E\{RSS/(n-p-1)\}=\sigma^2， \]

如果$H_0$同时也为真，那么有

\[ E\{(TSS-RSS)/p\}=\sigma^2。 \]

当预测变量和响应之间没有关系时，$F\approx 1$。若$H_a$为真，那么$E\{(TSS-RSS)/p\}>\sigma^2$，因此$F > 1$。

advertising <- read.csv('Advertising.csv')
lm.advertising <- lm(Sales ~ TV + Radio + Newspaper, data = advertising)
summary(lm.advertising)

# 结果如下：
Call:
lm(formula = Sales ~ TV + Radio + Newspaper, data = advertising)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.8277 -0.8908  0.2418  1.1893  2.8292 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.938889   0.311908   9.422   <2e-16 ***
TV           0.045765   0.001395  32.809   <2e-16 ***
Radio        0.188530   0.008611  21.893   <2e-16 ***
Newspaper   -0.001037   0.005871  -0.177     0.86    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.686 on 196 degrees of freedom
Multiple R-squared:  0.8972,    Adjusted R-squared:  0.8956 
F-statistic: 570.3 on 3 and 196 DF,  p-value: < 2.2e-16

上例中$F\approx 570$，远大于$1$，这表明销售额至少与一种媒体广告有关。

通过$F$拒绝$H_0$还需要考虑$n$和$p$的值。当$n$很大时，略大于$1$的$F$值也可拒绝$H_0$；当$n$很小时，$F$的值远大于$1$才可拒绝$H_0$。如果$H_0$为真，$\epsilon_i$是正态分布，那么$F$统计量服从$F$分布。当$n$和$p$已知，可以利用$F$分布计算$p$值（$p$-value），然后根据$p$值判断是否拒绝$H_0$。上例中，最后一个参数的$p$-value$\approx 0$，这表明销售额至少与一种媒体广告有关。

有时需要测试其中的$q$个系数是否为$0$，

\[ H_0:\beta_{p-q+1}=\beta_{p-q+2}=\cdots=\beta_p=0， \]

这可以通过对需忽略的$q$个数据排在最后实现。在这种情况下，线性回归模型需要去除最后的$q$个数据，此时的残差平方和记为$RSS_0$，那么$F$统计量为

\begin{equation} F={(RSS_0-RSS)/q\over RSS/(n-p-1)}。 \label{eq:f-statistic2} \end{equation}

单变量回归表明报纸广告与销售额也有关，但多变量回归表明报纸广告与销售额关系不大。对于单变量和多变量回归，电视的系数很相近，但是报纸的系数差异很大。虽然报纸对销售额的直接贡献很小，但从下面的相关系数可以看出，报纸和广播的相关系数高达$0.35$，也就是报纸投放的广告多广播也会多，从而报纸能间接提高销售额⁵。报纸可作为广播广告的替代品。报纸是通过广播对销售额的促进而获益的。

> cor(advertising[,2:5])

                  TV      Radio  Newspaper     Sales
TV        1.00000000 0.05480866 0.05664787 0.7822244
Radio     0.05480866 1.00000000 0.35410375 0.5762226
Newspaper 0.05664787 0.35410375 1.00000000 0.2282990
Sales     0.78222442 0.57622257 0.22829903 1.0000000

利用多变量回归分析每个变量的$t$统计量和$p$值，可以分析它们对销售额的影响，这等价于$q=1$时分别剔除每个变量的$F$统计量。这揭示了加入变量到模型的偏效应（partial effect）⁶。$p$值表明电视和广播的广告与销售额相关，但是没有证据表明在有电视和广播的广告时，报纸广告与销量相关。

为什么有了每个$p$值还需要全局的$F$统计量？如果存在一个很小的$p$值，那么至少有一个变量与响应相关——这个逻辑是有瑕疵的，特别是变量个数$p$很大时。例如：当$p=100$时，$H_0$为真，事实上没有那个变量与响应相关。在这种情况下，仍然大约$5\%$变量的$p$值会小于$0.05$，也就是即使在变量与响应无关的情况下也可能发现大约5个很小的$p$值，几乎可以保证至少存在一个$p$值小于$0.05$。因此，用每个变量对应的$t$统计量和$p$值判断变量与响应是否相关很可能得到错误的结果。但是，$F$统计量不存在这样的问题。如果$H_0$为真，不论变量有多少，只有$5\%$的机会使得$F$统计量会导致存在一个小于$5\%$的$p$值。

用$F$统计量判断变量与响应是否相关的条件是$p$（变量数）较小，一定要比$n$小。当$p>n$时，不能用最小二乘法拟合回归模型，因此也不能使用$F$统计量和前文讨论的其它概念。当$p$很大时可以采用前向选择（forward selection）。

好的回归模型要求$F$较大，$R^2$较大，在这种情况下，小的$p$值表明变量与响应相关。

编程实例

R处理因子类型变量

数据集中pisa2009test.csv和pisa2009train.csv中，raceeth的数据类型是因子，线性模型不能直接使用该数据。raceeth有7种取值White、American Indian/Alaska Native、Asian、Black、Hispanic、More than one race、Native Hawaiian/Other Pacific Islander。R会选择出现频率最高的作为参考，本例中为White，自动将raceeth用6个变量替代，形如raceethAmerican Indian/Alaska Native。如果raceeth原来取值为White，那么新的6个变量取值均为0；如果原来取值为American Indian/Alaska Native，那么新变量raceethAmerican Indian/Alaska Native取值为1，其余5个取值为0。这也是一种特征变换策略。

pisaTrain <- read.csv('Data/pisa2009train.csv')
pisaTest <- read.csv('Data/pisa2009test.csv')

# 忽略NA数据
pisaTrain <- na.omit(pisaTrain)
pisaTest <- na.omit(pisaTest)

# White出现频率最高，排在最前面
pisaTrain$raceeth = relevel(pisaTrain$raceeth, "White")
pisaTest$raceeth = relevel(pisaTest$raceeth, "White")

# 带因子（factor）变量的线性回归
lmScore <- lm(readingScore ~ ., data = pisaTrain)

step逐步回归

R中的step是MASS包中stepAIC函数的简版，它通过AIC准则进行模型选择，可以得到减缓的模型。

# 接上例
lmScoreSimple <- step(lmScore)

参考资料

脚注