UA MATH571A 多元线性回归II 变量选择

多项式回归与交互项回归阶数的确定含质量型变量的回归含质量型变量的交互项二值变量与二值变量的交互项二值变量与数量型变量的交互项变量选择的准则Rp2R^2_pRp2​与Ra,p2R^2_{a,p}Ra,p2​Mallow's CpC_pCp​AICpAIC_pAICp​与SBCpSBC_pSBCp​PRESSpPRESS_pPRESSp​模型选择的算法向前逐步回归算法（Forward Stepwise Regression）向后消元算法（Backward Elimination）模型验证（Validation）MSPR数据分割（Data Splitting）

使用多元线性回归建模时，假设X1,X2,...,Xp−1X_1,X_2,...,X_{p-1}X1​,X2​,...,Xp−1​是潜在的解释变量，我们总是想要知道到底应该将哪些变量，以什么形式（阶数或交互项等）加入到回归模型中。因此这篇博文首先介绍多项式回归与交互项回归，以及怎么用序贯方差分析的方法判断是否需要考虑高阶项与交互项。然后介绍选择解释变量的方法以及构造最优模型的算法。

多项式回归与交互项回归

以一元线性回归Y=β0+β1X1+ϵY = \beta_0 + \beta_1X_1+\epsilonY=β0​+β1​X1​+ϵ为例，假设引入X1X_1X1​的二阶项，则模型变为二次一元线性回归

Y=β0+β1X1+β11X12+ϵY = \beta_0 + \beta_1X_1 + \beta_{11}X_1^2+\epsilon Y=β0​+β1​X1​+β11​X12​+ϵ

系数β11\beta_{11}β11​的含义是

∂2Y∂X12=β11\frac{\partial^2 Y}{\partial X_1^2} = \beta_{11} ∂X12​∂2Y​=β11​

代表回归曲线的凸性，如果β11<0\beta_{11}<0β11​<0，回归曲线上凸，如果β11>0\beta_{11}>0β11​>0，回归曲线下凸。假设引入三阶项，则模型变为三次一元线性回归

Y=β0+β1X1+β11X12+β111X13+ϵY = \beta_0 + \beta_1X_1 + \beta_{11}X_1^2 + \beta_{111}X_1^3+\epsilon Y=β0​+β1​X1​+β11​X12​+β111​X13​+ϵ

但通常回归模型为了描述凸性设定到二阶就足够了。在做多项式回归时通常要将样本数据中心化，即定义x1=X1−Xˉ1x_1 = X_1 - \bar{X}_1x1​=X1​−Xˉ1​，然后用x1x_1x1​来回归。这是因为X1X_1X1​与X12X_1^2X12​等高阶项通常都是高度相关的，中心化后二者的相关性能够被削弱。以二元线性回归Y=β0+β1X1+β2X2+ϵY = \beta_0 + \beta_1X_1 + \beta_2 X_2+\epsilonY=β0​+β1​X1​+β2​X2​+ϵ为例，假设引入解释变量X1X_1X1​与X2X_2X2​的交互项（Interaction Term），则模型变为二次二元线性回归

Y=β0+β1X1+β2X2+β12X1X2+ϵY = \beta_0 + \beta_1X_1 + \beta_2 X_2 + \beta_{12} X_1X_2+ \epsilon Y=β0​+β1​X1​+β2​X2​+β12​X1​X2​+ϵ

含有交互项的回归模型系数的解释与一般的线性回归模型有所不同。

∂Y∂X1=β1+β12X2∂2Y∂X1∂X2=β12\frac{\partial Y}{\partial X_1} = \beta_1 + \beta_{12}X_2 \\ \frac{\partial^2 Y}{\partial X_1 \partial X_2} = \beta_{12} ∂X1​∂Y​=β1​+β12​X2​∂X1​∂X2​∂2Y​=β12​

从上两式可以看出，有交互项时，X1X_1X1​对YYY的效应（∂Y∂X1\frac{\partial Y}{\partial X_1}∂X1​∂Y​）可以分为直接效应（β1\beta_1β1​）与通过X2X_2X2​对YYY造成影响的间接效应（β12X2\beta_{12}X_2β12​X2​），系数β12\beta_{12}β12​衡量的是单位间接效应。

阶数的确定

考虑一个线性回归模型要不要引入高次项或者交互项的方法也是序贯方差分析。假设已有的多元线性回归是

Y=β0+β1X1+β2X2+ϵY = \beta_0 + \beta_1X_1 + \beta_2 X_2+\epsilonY=β0​+β1​X1​+β2​X2​+ϵ

如果要考虑是否需要引入二阶项X12X_1^2X12​与X1X2X_1X_2X1​X2​，只需用序贯方差分析中检验部分系数的方法，在回归

Y=β0+β1X1+β2X2+β11X12+β12X1X2+ϵY = \beta_0 + \beta_1X_1 + \beta_2 X_2+ \beta_{11}X_1^2 + \beta_{12}X_1X_2 + \epsilonY=β0​+β1​X1​+β2​X2​+β11​X12​+β12​X1​X2​+ϵ中（这个回归是完整模型FM）检验

H0:β11=β12=0Ha:notbothofβ11,β12equaltozeroH_0: \beta_{11}=\beta_{12}= 0 \\ H_a: not\ both\ of\ \beta_{11}, \beta_{12}\ equal\ to\ zero H0​:β11​=β12​=0Ha​:notbothofβ11​,β12​equaltozero

定义统计量

F∗=SSR(X12,X1X2∣X1,X2)2/SSE(FM)dfFM∼F(2,N−5)F^* = \frac{SSR(X_1^2,X_1X_2|X_1,X_2) }{2}/\frac{SSE(FM)}{df_{FM}} \sim F(2,N-5) F∗=2SSR(X12​,X1​X2​∣X1​,X2​)​/dfFM​SSE(FM)​∼F(2,N−5)

假设检验水平为α\alphaα，若F∗≤F(1−α;2,N−5)F^*\le F(1-\alpha;2,N-5)F∗≤F(1−α;2,N−5)，接受原假设，无需在模型中加入X12X_1^2X12​与X1X2X_1X_2X1​X2​，若F∗>F(1−α;2,N−5)F^*>F(1-\alpha;2,N-5)F∗>F(1−α;2,N−5)，拒绝原假设。

含质量型变量的回归

通常二值变量取值为0或1，代表的含义是是或否。比如定义代表性别的变量

F=0,1F=0,1F=0,1，等于1代表对应的样本是女性，等于0代表对应的样本是男性。当类别超过两类时，有两种方式来定义解释变量。比如车辆行驶方向可以是直行，左转，右转，要记录一个区域内的车辆行驶方向状态，可以定义变量D=0,1,2D=0,1,2D=0,1,2，等于0代表直行，等于1代表左转，等于2代表右转。赋予不同的数值不会影响对回归系数的推断。也可以定义两个二值变量D1=0,1D_1=0,1D1​=0,1，D2=0,1D_2=0,1D2​=0,1。D1D_1D1​取1代表车辆左转，D2D_2D2​取1代表车辆右转，二者均为零代表车辆直行。也就是说一个MMM值的质量型变量可以被M−1M-1M−1个二值变量代替。

含质量型变量的交互项

这里用一个例子来说明怎么解释含质量型变量的交互项的系数。血浆肌酸（PCr）是临床上用来判断肾病的指标之一，现在想要考察性别（用二值变量F表示）、血压、BMI对PCr的影响。血压用二值变量BP表示，等于1代表高血压，等于0代表没有高血压。

二值变量与二值变量的交互项

首先考察回归模型，

PCr=β0+β1F+β2BP+β12F×BP+ϵPCr = \beta_0 + \beta_1 F + \beta_2 BP + \beta_{12} F \times BP + \epsilon PCr=β0​+β1​F+β2​BP+β12​F×BP+ϵ

则β0\beta_0β0​代表没有高血压的男性血浆肌酸的平均值，β0+β1\beta_0+\beta_1β0​+β1​代表没有高血压的女性血浆肌酸的平均值，因此β1\beta_1β1​代表没有高血压的男性和女性血浆肌酸平均值之差。β0+β2\beta_0+\beta_2β0​+β2​代表有高血压的男性血浆肌酸的平均值，β0+β1+β2+β12\beta_0+\beta_1+\beta_2+\beta_{12}β0​+β1​+β2​+β12​代表有高血压的女性血浆肌酸的平均值，β1+β12\beta_1+\beta_{12}β1​+β12​代表有高血压的男性和女性血浆肌酸的平均值之差。所以如果我们检验系数β12\beta_{12}β12​，那么我们想要探索的是有无高血压是否会影响男女平均血浆肌酸的差异。

二值变量与数量型变量的交互项

现在考虑另一个回归模型，

PCr=β0+β1BMI+β2BP+β12BMI×BP+ϵPCr = \beta_0 + \beta_1 BMI + \beta_2 BP + \beta_{12} BMI \times BP + \epsilon PCr=β0​+β1​BMI+β2​BP+β12​BMI×BP+ϵ

在这个模型中，BMI对PCr的效应为

∂PCr∂BMI=β1+β12BP\frac{\partial PCr}{\partial BMI} = \beta_1 + \beta_{12}BP ∂BMI∂PCr​=β1​+β12​BP

因此BMI对PCr的效应在高血压与非高血压人群间的差异是β12\beta_{12}β12​。对高血压人群

PCr=(β0+β2)+(β1+β12)BMI+ϵPCr = (\beta_0 + \beta_2)+( \beta_1+ \beta_{12} ) BMI+ \epsilon PCr=(β0​+β2​)+(β1​+β12​)BMI+ϵ

对非高血压人群

PCr=β0+β1BMI+ϵPCr = \beta_0 + \beta_1 BMI+ \epsilon PCr=β0​+β1​BMI+ϵ

因此对于不同人群，BMI对PCr效应的差异为

ΔPCr=β2+β12BMI\Delta PCr = \beta_2 + \beta_{12} BMI ΔPCr=β2​+β12​BMI

这个式子可以解释为高血压和非高血压人群PCr差异由两部分构成，直接差异β2\beta_2β2​与通过BMI带来的间接差异β12BMI\beta_{12} BMIβ12​BMI，β12\beta_{12}β12​也可以解释为BMI对高血压与非高血压人群PCr差异的效应。多指变量系数的解释与二值变量的类似。

变量选择的准则

假设我们已经有了P−1P-1P−1个潜在的解释变量，并且想用这些解释变量构造一个二阶线性回归模型，则模型最多可能包括一个常数项，P−1P-1P−1个一阶项，P−1P-1P−1个平方项，CP−12C_{P-1}^2CP−12​个交互项，共CP+12C_{P+1}^2CP+12​项，因此所有可能的模型数量为2CP+122^{C_{P+1}^2}2CP+12​，由此可见当潜在解释变量的数量较多时，比如8个（P−1=8P-1=8P−1=8），所有可能的模型数量超过三万亿个。如果只考虑一阶的模型，那么所有可能的模型数量也有2P−1=2562^{P-1}=2562P−1=256个，所以有必要在寻找最佳的回归模型之前先对变量进行选择。假设考察的模型包含这P−1P-1P−1个潜在解释变量中的p−1p-1p−1个（0≤p−1≤P−10 \le p-1 \le P-10≤p−1≤P−1）,以下指标可以用来衡量这p−1p-1p−1个解释变量构成的模型的优劣。

Rp2R^2_pRp2​与Ra,p2R^2_{a,p}Ra,p2​

仿照可决系数，定义这p−1p-1p−1个解释变量的可决系数为

Rp2=1−SSEpSSTR^2_p = 1 - \frac{SSE_p}{SST} Rp2​=1−SSTSSEp​​

其中SSEpSSE_pSSEp​是p−1p-1p−1个解释变量的线性回归模型的残差平方和，SSEpSSE_pSSEp​越小说明不能被这p−1p-1p−1个解释变量解释的信息就越少，模型质量就越高，因此Rp2R^2_pRp2​越大说明这p−1p-1p−1个解释变量的解释力越强。与一元线性回归可决系数定义的缺陷类似，此处同样需要考虑模型的自由度，因此定义这p−1p-1p−1个解释变量的调整可决系数为

Ra,p2=1−SSEp/N−pSST/N−1R^2_{a,p} = 1 - \frac{SSE_p/N-p}{SST/N-1} Ra,p2​=1−SST/N−1SSEp​/N−p​

Mallow’s CpC_pCp​

假设μi\mu_iμi​是第i个样本的真实的mean response，则第i个样本的拟合误差为Y^i−μi\hat{Y}_i-\mu_iY^i​−μi​。计算

E[(Y^i−μi)2]=E[(Y^i−E(Y^i)+E(Y^i)−μi)2]=E[Y^i−E(Y^i)]2+E[E(Y^i)−μi]2+2E[(Y^i−E(Y^i))(E(Y^i)−μi)]=σ2(Y^i)+[E(Y^i)−μi]2E[(\hat{Y}_i-\mu_i)^2] = E[(\hat{Y}_i-E(\hat{Y}_i) +E(\hat{Y}_i)- \mu_i)^2] \\=E[\hat{Y}_i-E(\hat{Y}_i) ]^2 + E[E(\hat{Y}_i)- \mu_i]^2 + 2E[(\hat{Y}_i-E(\hat{Y}_i))(E(\hat{Y}_i)- \mu_i)]\\=\sigma^2(\hat{Y}_i) + [E(\hat{Y}_i)- \mu_i]^2 E[(Y^i​−μi​)2]=E[(Y^i​−E(Y^i​)+E(Y^i​)−μi​)2]=E[Y^i​−E(Y^i​)]2+E[E(Y^i​)−μi​]2+2E[(Y^i​−E(Y^i​))(E(Y^i​)−μi​)]=σ2(Y^i​)+[E(Y^i​)−μi​]2

定义

Γp=1σ2∑i=1N{σ2(Y^i)+[E(Y^i)−μi]2}\Gamma_p = \frac{1}{\sigma^2} \sum_{i=1}^N \{ \sigma^2(\hat{Y}_i) + [E(\hat{Y}_i)- \mu_i]^2 \} Γp​=σ21​i=1∑N​{σ2(Y^i​)+[E(Y^i​)−μi​]2}

其中σ2\sigma^2σ2是模型真实的方差。Γp\Gamma_pΓp​越小代表选定的p−1p-1p−1个解释变量的解释力越强。Mallow’s CpC_pCp​是Γp\Gamma_pΓp​的估计量

Cp=SSEpMSEP−(N−2p)C_p = \frac{SSE_p}{MSE_P} - (N-2p) Cp​=MSEP​SSEp​​−(N−2p)

其中MSEPMSE_PMSEP​是P−1P-1P−1个解释变量的回归模型的均方误差，并且ECp≈pEC_p \approx pECp​≈p。显然当p=Pp=Pp=P时，

CP=SSEPMSEP−(N−2P)=SSEPSSEP/N−P−(N−2P)=PC_P = \frac{SSE_P}{MSE_P} - (N-2P) = \frac{SSE_P}{SSE_P/N-P} - (N-2P)=P CP​=MSEP​SSEP​​−(N−2P)=SSEP​/N−PSSEP​​−(N−2P)=P

下面给出Mallow’s CpC_pCp​是Γp\Gamma_pΓp​的估计量的证明思路。首先计算得到

∑i=1Nσ2(Y^i)=pσ2\sum_{i=1}^N \sigma^2(\hat{Y}_i) = p\sigma^2 i=1∑N​σ2(Y^i​)=pσ2

然后根据（之前推导MSE是方差的无偏估计时计算过MSE的期望，SSE的与之类似）

E(SSEp)=∑i=1N[E(Y^i)−μi]2+(N−p)σ2E(SSE_p) = \sum_{i=1}^N[E(\hat{Y}_i)- \mu_i]^2 + (N-p)\sigma^2 E(SSEp​)=i=1∑N​[E(Y^i​)−μi​]2+(N−p)σ2

可以得到

Γp=E(SSEp)−(N−p)σ2+pσ2σ2=E(SSEp)σ2−(N−2p)\Gamma_p = \frac{E(SSE_p) - (N-p)\sigma^2 + p\sigma^2}{\sigma^2} =\frac{E(SSE_p) }{\sigma^2} - (N-2p) Γp​=σ2E(SSEp​)−(N−p)σ2+pσ2​=σ2E(SSEp​)​−(N−2p)

用MSEPMSE_PMSEP​作为σ2\sigma^2σ2的无偏估计替换掉它就可以得到Mallow’s CpC_pCp​。

AICpAIC_pAICp​与SBCpSBC_pSBCp​

AIC（Akaike Information Criteria）与SBC（Schwartz Bayesian Criteria）是很常用的准则，这里直接给出公式：

AICp=NlnSSEp−NlnN+2pSBCp=NlnSSEp−NlnN+[lnN]pAIC_p = NlnSSE_p - NlnN + 2p \\ SBC_p = NlnSSE_p - NlnN + [lnN]p AICp​=NlnSSEp​−NlnN+2pSBCp​=NlnSSEp​−NlnN+[lnN]p

这两个表达式中NlnNNlnNNlnN只受样本量的影响，与变量选择无关。第一项lnSSEplnSSE_plnSSEp​越小说明模型解释力越强。第三项是惩罚项，变量数量越多，惩罚项就会越大。这两个指标都是越小越好，即总是希望用尽量少的解释变量去获得更强的解释力。样本数超过8时，SBCpSBC_pSBCp​的惩罚项总是比AICpAIC_pAICp​的更大，因此相比AICpAIC_pAICp​，SBCpSBC_pSBCp​倾向于选出包含更少解释变量的模型。

PRESSpPRESS_pPRESSp​

PRESS的全称是prediction sum of square，但在没有新数据的情况下没有办法评估真正的预测能力，所以这里采用的是LOO（Leave-one-out）的思想。用Y^i(i)\hat{Y}_{i(i)}Y^i(i)​表示用除了第i个样本以外的其他N−1N-1N−1个样本估计模型，并用以计算所得的第i个样本的预测值。从而

PRESSp=∑i=1N(Yi−Y^i(i))2PRESS_p = \sum_{i=1}^{N} (Y_i - \hat{Y}_{i(i)})^2 PRESSp​=i=1∑N​(Yi​−Y^i(i)​)2

与之前三种基于SSEpSSE_pSSEp​的指标不同，PRESSpPRESS_pPRESSp​关注的不是模型的拟合能力，而是模型的预测能力；但与SSEpSSE_pSSEp​类似的是，PRESSpPRESS_pPRESSp​越小说明模型的预测能力越强。

模型选择的算法

回到最初的问题中，也就是怎么在2CP+122^{C_{P+1}^2}2CP+12​个潜在的模型中找到最优的二阶模型或者是在2P−12^{P-1}2P−1个潜在的模型中找到最优的一阶模型。现在我们有四种判断一个模型优劣的指标，接下来可以利用这些指标构建一些寻找最优模型的算法。有一个比较显然的事实是，如果算法需要通过遍历所有可能的模型来寻找最优模型，那么这个算法是一个NP的算法，只有在PPP比较小的时候才能在较快完成计算，PPP比较大时需要的时间就会非常长。例如有8个潜在的解释变量时，假设计算每个模型的CpC_pCp​需要0.1秒，那么遍历完所有可能的模型需要3千亿秒以上。显然这样的算法没有实用价值。Best Subsets算法指的是不需要遍历所有可能模型就能够找出按某种指标衡量的最好的几个模型的那些算法。但是当PPP更大，超过了四十的时候，这些算法也需要很多时间。以下两种算法即使在PPP很大的时候也能在比较短的时间内完成。

向前逐步回归算法（Forward Stepwise Regression）

第一步，用P−1P-1P−1个潜在的解释变量对被解释变量YYY做一元线性回归，对每一个线性回归，计算t统计量，选出t统计量的值最大的解释变量，若它的t统计量超过了预定义检验水平下的t值，则该解释变量应该被加入到模型中；

第二步，用第一步筛选出来的解释变量与其他所有解释变量分别组成二元回归模型，计算第二个解释变量的t统计量，选出t统计量的值最大的解释变量，若它的t统计量超过了预定义检验水平下的t值，则该解释变量应该被加入到模型中；

第三步，对于已经被添加到模型中的解释变量，计算其一元回归的t统计量，剔除掉t统计量最小且未达到检验水平解释变量；

第四步，考虑被添加在模型中且被保留下来的解释变量外的其他变量是否应该被添加到模型中，并再次考虑已经被添加到模型中的解释变量是否应该被剔除。·

向后消元算法（Backward Elimination）

向后算法与向前算法的顺序正好相反，向后算法是先用所有解释变量做多元回归，剔除掉t值最小的解释变量，再对剩下的解释变量做多元回归，剔除掉t值最小的解释变量，直到剩下的解释变量t值均高于检验水平为止。

模型验证（Validation）

通常有三种验证模型质量的方法：1）收集新数据用来检验模型的预测能力；2）与理论结果、以前的实证结果以及模拟的结果对比；3）将样本分成训练集和测试集，训练集用来估计模型，测试集用来验证模型的质量。

MSPR

建好的多元回归模型可以用来预测测试集或者新数据集。定义预测均方误（mean squared of predictive error）

MSPR=∑i=1N∗(Yi−Y^i)2N∗MSPR = \frac{\sum_{i=1}^{N^*} (Y_i - \hat{Y}_i)^2}{N^*} MSPR=N∗∑i=1N∗​(Yi​−Y^i​)2​

N∗N^*N∗是训练集或者新数据集的样本数量。预测均方误越小表示模型的预测能力越强。但事实上新数据的收集十分困难，因此大部分时候是采用把完整的数据集分成训练集和测试集的形式。

数据分割（Data Splitting）

数据分割就是把数据集分成训练集（model build-in set or training set）和测试集（validation set or prediction set）的方法，训练集用来估计模型的参数，调节模型中的超参，测试集用来检验模型的优劣。这种评估模型的方式叫做交叉验证（Cross Validation）。Data splitting的目标是，训练集要足够大，包含足够多的信息使模型具有较好的泛化能力；测试集也需要足够多的不同的case使得模型的泛化能力能够得到公正的评估。针对不同的样本寻找最优的data splitting的方案到现在也是一个值得研究的话题。