第3章
CHAPTER 3


多元线性回归






实际应用中，一个自变量同时受多个因变量的影响的情况非常普遍。因此考虑将第2章中介绍的一元线性回归拓展到多元的情形。包括多个解释变量的回归模型，就称为多元回归模型。多元线性回归分析是一元情况的简单推广，读者应该注意建立二者之间的联系。
3.1多元线性回归模型
假设因变量y与m个解释变量x1,x2,…,xm具有线性相关关系，取n组观察值，则总体线性回归模型为
yi=w0+w1xi1+w2xi2+…+wmxim+ui,i=1,2,…,n
包含m个解释变量的总体回归模型也可以表示为
E(y|xi1,xi2,…,xim)=w0+w1x1i+w2x2i+…+wmxim,i=1,2,…,n
上式表示在给定xi1,xi2,…,xim的条件下，y的条件均值或数学期望。特别地，我们称w0是截距，w1,w2,…,wm是偏回归系数。偏回归系数又称为偏斜率系数。例如，w1度量了在其他解释变量x2,x3,…,xm保持不变的情况下，x1每变化1个单位时，y的均值E(y|xi1,xi2,…,xim)的变化。换句话说，w1给出了其他解释变量保持不变时，E(y|xi1,xi2,…,xim)对x1的斜率。
不难发现，多元线性回归模型是以多个解释变量的固定值为条件的回归分析。
同一元线性回归模型一样，多元线性总体回归模型是无法得到的。所以我们只能用样本观察值进行估计。对应于前面给出的总体回归模型可知多元线性样本回归模型为
y^i=w^0+w^1xi1+w^2xi2+…+w^mxim,i=1,2,…,n
和
yi=w^0+w^1xi1+w^2xi2+…+w^mxim+ei,i=1,2,…,n
其中，y^i是总体均值E(y|xi1,xi2,…,xim)的估计; w^j是总体偏回归系数wj的估计，j=1,2,…,m； 残差项ei是对随机项ui的估计。
对多元线性总体回归方模型可以用线性方程组的形式表示为
y1=w0+w1x11+w2x12+…+wmx1m+u1
y2=w0+w1x21+w2x22+…+wmx2m+u2
︙
yn=w0+w1xn1+w2xn2+…+wmxnm+un
将上述方程组改写成矩阵的形式:
y1
y2
︙
yn=1x11x12…x1m
1x21x22…x2m
︙︙︙︙
1xn1xn2…xnmw0
w1
︙
wm+u1
u2
︙
un
或者写成如下形式:
y=Xw+u
上式就是用矩阵形式表示的多元线性总体回归模型。其中y为n阶因变量观察值向量； X表示n×m阶解释变量的观察值矩阵； u表示n阶随机扰动项向量； w表示m阶总体回归参数向量。




同理，可以得到多元线性样本回归模型的矩阵表示为
y=Xw^+e
或者
y^=Xw^
其中y^表示n阶因变量回归拟合值向量； w^表示m阶回归参数w的估计值向量； e表示n阶残差向量。
以上各向量的完整形式如下：
y=y1
y2
︙
yn,y^=y^1
y^2
︙
y^n,w^=w^0
w^1
︙
w^m,e=e1
e2
︙
en
显而易见的是，由于解释变量数量的增多，多元线性回归模型的计算要比一元的情况复杂很多。最后与一元线性回归模型一样，为了对回归模型中的参数进行估计，要求多元线性回归模型在满足线性关系之外还必须遵守以下假定。
1.  零均值假定
干扰项ui均值为零，或对每一个i，都有E(ui|xi1,xi2,…,xim)=0。
2.  同方差假定
干扰项ui的方差保持不变，即var(ui)=σ2。为了进行假设检验，我们通常认为随机扰动(噪声)符合一个均值为0、方差为σ2的正态分布，即ui~N(0,σ2)。
3.  相互独立性
随机扰动项彼此之间都是相互独立的，即cov(ui,uj)=0，其中i≠j。
4.  无多重共线性假定
解释变量之间不存在精确的线性关系，即没有一个解释变量可以被写成模型中其余解释变量的线性组合。
3.2多元回归模型估计
为了建立完整的多元回归模型，我们需要使用最小二乘法对模型中的偏回归系数进行估计，这个过程中的所用到的许多性质与一元情况下一致。
3.2.1最小二乘估计量
已知多元线性样本回归模型为
yi=w^0+w^1xi1+w^2xi2+…+w^mxim+ei,i=1,2,…,n
于是离差平方和为
∑e2i=∑(yi- y^i)2=∑(yi-w^0-w^1xi1-w^2xi2-…-w^mxim)2
现在求估计的参数w^0,w^1,…,w^m，使得离差平方和取得最小值，于是根据微积分中极值存在的条件，要解方程组
∑e2iw0=-2∑(yi-w^0-w^1xi1-…-w^mxim)=0
∑e2iw1=-2∑(yi-w^0-w^1xi1-…-w^mxim)xi1=0
︙
∑e2iwm=-2∑(yi-w^0-w^1xi1-…-w^mxim)xim=0
其解就是参数w0,w1,…,wm的最小二乘估计w^0,w^1,…,w^m。
将以上方程组改写成
nw^0+∑w^1xi1+∑w^2xi2+…+∑w^mxim=∑yi
∑w^0xi1+∑w^1x2i1+∑w^2xi1xi2+…+∑w^mxi1xim=∑xi1yi
︙
∑w^0xim+∑w^1ximxi1+∑w^2ximxi2+…+∑w^mx2im=∑ximyi
这个方程组称为正规方程组。为了把正规方程组改写成矩阵形式，记系数矩阵为A，常数项向量为b、w的估计值向量为w^，即
A=n∑xi1∑xi2…∑xim
∑xi1∑x2i1∑xi1xi2…∑xi1xim
︙︙︙︙
∑xim∑ximxi1∑ximxi2…∑x2im

=111…1
x11x21x31…xn1
︙︙︙︙
x1mx2mx3m…xnm1x11x12…x1m
1x21x22…x2m
︙︙︙︙
1xn1xn2…xnm=XTX

b=∑yi
∑xi1yi
︙
∑ximyi=11…1
x11x21…xn1
︙︙︙
x1mx2m…xnmy1
y2
︙
yn=XTy
其中w^=(w^0,w^1,…，w^m)T，y=(y1,y2,…,yn)T。所以正规方程组可以表示为
Aw^=b或(XTX)w^=XTy
当系数矩阵可逆时，正规方程组的解为
w^=A-1b=(XTX)-1XTy
进而还可以得到
y^=Xw^=X(XTX)-1XTy
令H=X(XTX)-1XT，则有y^=Hy，H是一个n阶对称矩阵，通常称为帽子矩阵。该矩阵的对角线元素记为hii，它给出了第i个观测值离其余n－1个观测值的距离有多远，我们通常称其为杠杆率。
3.2.2多元回归的实例
现在将通过一个实例来演示在Python中建立多元线性回归模型的方法。根据经验，沉淀物吸收能力是土壤的一项重要特征，因为它会影响杀虫剂和其他各种农药的有效性。在一项实验中，我们测定了若干组土壤样本的情况，数据如表31所示。其中，y表示磷酸盐吸收指标； x1和x2分别表示可提取的铁含量与可提取的铝含量。请根据这些数据建立y关于x1和x2的多元线性回归方程。


表31土壤沉淀物吸收能力采样数据




x1
61
175
111
124
130
173
169
169
160
244
257
333
199
x2
13
21
24
23
64
38
33
61
39
71
112
88
54
y
4
18
14
18
26
26
21
30
28
36
65
62
40

在进行一元线性回归分析之前往往会使用散点图来考查一下解释变量与被解释变量之间的线性关系。在进行多元线性回归分析时，也可以采用类似的图形来观察模型中解释变量与被解释变量间的关系，但所采用的统计图形要更复杂一些，它被称为是散点图阵列，如图31所示。



图31散点图阵列


从散点图中可以看出，每个解释变量都与被解释变量有一定的线性关系，而且这也是我们希望看到的。更重要的是，两个解释变量之间的线性关系并不显著，这就意味着出现多重共线性的可能性较低。在构建多元线性回归模型时，随着解释变量数目的增多，其中某两个解释变量之间产生多重共线性是很容易发生的情况。此时就需要考虑是否将其中某个变量从模型中剔除出去，甚至是重新考虑模型的构建。现在用以检验多重共线性的方法有很多，有兴趣的读者可以参阅其他相关著作，此处不再赘述。但本书后面还会向读者展示，现代回归分析是如何化解多重共线性的影响的。
下面的代码演示了在Python中构建多元线性回归模型的方法，其中用到了statsmodels这个库，它是我们在处理统计学问题时主要用到的一个模块。特别地，要开展基于普通最小二乘法的回归分析，需对OLS类进行实例化，而模型的训练过程则通过调用它的fit()方法来完成。



import numpy as np

import statsmodels.api as sm



pai = np.array(［4, 18, 14, 18, 26, 26, 21, 30, 28, 36, 65, 62, 40］)

iron = np.array(［61, 175, 111, 124, 130, 173, 169, 169, 160, 244, 257, 333, 199］)

aluminium = np.array(［13, 21, 24, 23, 64, 38, 33, 61, 39, 71, 112, 88, 54］)



X_mul = np.vstack((iron, aluminium))

X_mul = X_mul.transpose()

X_mul = sm.add_constant(X_mul)



model_multilin = sm.OLS(pai, X_mul).fit()

model_multilin.summary()




执行上述代码，系统将为刚刚构建好的模型输出一份非常全面的报告，如图32所示。根据图32所示的结果(虚线框标出的部分)，可以写出多元线性回归方程如下： 
y^=-7.3507+0.1127x1+0.349x2
同时，易见回归模型的拟合优度R2=0.948，调整判定系数R2adj=0.938，说明模型的拟合效果较好。这两个指标的意义在一元的情况下已经进行过详细的介绍，此处不再赘述。需要说明的是，在多元回归的情况下，随着自变量个数的增多，拟合优度也会提高，所以仅仅看这一个指标说服力是有限的。具体这些判定指标的意义，本章后面还会做进一步的解读。



图32多元线性回归输出结果


3.2.3总体参数估计量
由w的估计量w^的表达式可见，w^的每一个分量都是相互独立且服从正态分布的随机变量y1,y2,…,yn的线性组合，从而可知随机变量w^服从m+1维正态分布。为了求出w^的分布，首先来计算w^的期望和方差(或方差阵)。
向量w^的数学期望定义为
E(w^)=［E(w^0),E(w^1),…,E(w^m)］T
而且对任意n×(m+1)阶矩阵A，容易证明
E(Aw^)=AE(w^)
于是可得
E(w^)=E［(XTX)-1XTy］=(XTX)-1XTE(y)

=(XTX)-1XTE(Xw+u)=(XTX)-1XTXw=w
所以w^是w的无偏估计，即w^0,w^1,…,w^m依次是w0,w1,…,wm的无偏估计，为了计算w^的方差阵，先把方差阵写成向量乘积的形式：
D(w^)=D(w^0)cov(w^0,w^1)…cov(w^0,w^m)
cov(w^1,w^0)D(w^1)…cov(w^1,w^m)
︙︙︙
cov(w^m,w^0)cov(w^m,w^1)…D(w^m)

=E{［(w^0-E［w^0］),(w^1-E［w^1］),…,(w^m-E［w^m］)］T×
［(w^0-E［w^0］),(w^1-E［w^1］),…,(w^m-E［w^m］)］}

=E{［w^-E(w^)］［w^-E(w^)］T}
而且
计算过程中用到的一些矩阵计算性质如下，其中A、B是两个可以做乘积的矩阵，I是单位矩阵，则有(AB)T=BTAT，AA-1=I，AT(A-1)T=(A-1A)T=I(AT)-1=(A-1)T。
E{［w^-E(w^)］［w^-E(w^)］T}

=E{［(XTX)-1XT(y-E(y))］［(XTX)-1XT(y-E(y))］T}

=E［(XTX)-1XT(y-E(y))(y-E(y))TX(XTX)-1］

=(XTX)-1XTE［(y-E(y))(y-E(y))T］X(XTX)-1

=(XTX)-1XTE(uuT)X(XTX)-1

=(XTX)-1XTσ2IX(XTX)-1=σ2(XTX)-1
根据已经得到的计算结果，易知w^的方差阵等于σ2A-1，这个方差阵给出了w^中每个元素(即w^0,w^1,…,w^m)的方差(或标准差)，以及元素之间的协方差。当i=j时，矩阵对角线上的元素就是相应w^i的方差var(w^i)=σ2A-1ij，由此也可知道w^i的标准差为
se(w^i)=σA-1ij
当i≠j时，矩阵对角线以外的元素就表示相应w^i与w^j的协方差，即cov(w^i,w^j)=σ2A-1ij。
例如，在土壤沉淀物吸收情况的例子中可以求得矩阵A如下：
A=XTX=102305641
2305467669133162
64113316241831
相应的逆矩阵A-1如下：
A-1=0.633138-0.0038260.002477
-0.0038260.000046-0.000088
0.002477-0.0000880.000265
而且从系统的输出中也知道残差标准误差为4.379，于是有
se(w^0)=4.379×0.633138≈3.48437

se(w^1)=4.379×0.000046≈0.02969

se(w^2)=4.379×0.000265≈0.07130
在考虑到计算过程中保留精度存在差异的条件下，上述参数的标准误差与3.2.2节中系统的输出结果是基本一致的。
注意Python中的残差标准误差(Residual Standard Error)其实就是残差的标准差(Residual Standard Deviation)，如果读者对于它的计算仍感到困惑，那么可以参看第2章中的相关结论。总的来说，在多元线性回归中，总体方差(同时也是误差项的方差)σ2的无偏估计量为
σ^2=∑e2in-k=∑(yi- y^i)2n-k
所以残差(或误差项)的标准差为
s=σ^=∑(yi- y^i)2n-k
其中，k是被估计参数的数量。
3.3从线代角度理解最小二乘
本节将从矩阵的角度来解释最小二乘问题。矩阵方法与微积分方法所得的结论都是一样的，但是基于矩阵的方法更加强大，而且更容易推广。在这个过程中，最关键的内容是理解最佳逼近原理。它是从矩阵(或线性代数)角度解释最小二乘问题的核心。
3.3.1最小二乘问题的通解
非空集合SV，其中V是一个内积空间。定义S⊥={x∈V： 〈x,y〉=0,y∈S}，称S⊥为S的正交补。
关于这个定义，有以下两点值得注意： 
（1） 对于V的任何子集S，S⊥是V的一个子空间。
（2） 令S是V的一个子空间，那么S∩S⊥={0}。如果x∈S∩S⊥，那么〈x,x〉=0。因此，x=0。当然，0∈S∩S⊥。
注意： 如果S不是V的一个子空间，那么第（2）条则不成立。
定理： 令W是内积空间V的一个有限维子空间，并令y∈V。则有
（1） !u∈W，以及z∈W⊥，使得y=u+z。
（2） 如果{�瘙
經1,�瘙
經2,…,�瘙
經k}是W的一个正交标准基，那么u=∑ki=1〈y,�瘙
經i〉�瘙
經i。


图33定理的几何意义



上述定理的几何意义也是非常明确的，如图33所示。向量u是属于W空间的，z是属于W的正交补空间的，而且W和其正交补空间都是子空间。一定存在u和z使得y=u+z。
该定理有诸多非常重要的应用，其中之一就是给出了下面这个“最佳逼近原理”，在后面讨论最小二乘问题时，还会再用到它。
假设W是Rn空间中的一个子空间，y是Rn中的任意向量，y^是y在W上的正交投影，那么y^是W中最接近y的点，也就是说，‖y-y^‖<‖y-�瘙
經‖对所有属于W又异于y^的�瘙
經成立。这个结论也称为最佳逼近原理。其中，向量y^称为W中元素对y的最佳逼近。对于给定元素y，可以被某个给定子空间中的元素代替或“逼近”，我们用‖y-�瘙
經‖表示从y到�瘙
經的距离，则可以认为是用�瘙
經代替y的“误差”，而最佳逼近原理说明误差在�瘙
經=y^处取得最小值。
最佳逼近原理乍看起来有些抽象，但其几何意义却是非常直观的，如图34所示。不仅如此，它的证明也非常容易，只要使用勾股定理即可。
‖y-�瘙
經‖2=‖y- y^∈W⊥+ y^-�瘙
經∈W‖2
因为y-y^⊥y^-�瘙
經，所以应用勾股定理可得
‖y- y^+ y^-�瘙
經‖2=‖y- y^‖2+‖y^-�瘙
經‖2≥‖y- y^‖2
当且仅当y^=�瘙
經时取等号。

现在回到一直在讨论的回归问题上。假设在某些时间点xi上，观察到一些数值yi，于是有一组二维数据点(xi,yi)，其中i=1,2,…,n。现在的任务是找到一条最能代表这些数据的直线。如图35所示，假设备选直线的方程为y=bx+c。在xi时刻，实际观测值为yi，而根据回归直线，估计值应该是bxi+c。



图34最佳逼近原理的几何解释





图35线性回归



因此，实际观测值与估计值之间的误差就可以表示为yi-cxi-d。既然我们的任务是找到那条“最好的”直线，那么就可以用“总误差”最小来作为衡量所谓“最好的”标准。而且等价地，我们用平方来取代绝对值计算，注意(xi,yi)是已知的，于是现在的目标变为求使得下式(总误差)最小的参数b和c： 
∑ni=1(yi-bxi-c)2
一种可以想到的解法就是采用微积分的方法，这也是本书前面一直采用的方法。下面来看看如何把上述优化问题转变为矩阵问题。令
X=x11
x21
︙︙
xn1,w=b
c,y=y1
y2
︙
yn
此时问题就变成如下形式： 
E=∑ni=1(yi-cxi-d)2=‖y-Xw‖2=‖Xw-y‖2
可见这个形式也简单明了。不仅如此，矩阵也非常便于问题的推广，比如说现在不是要找直线，而是一个用多项式所表示的曲线(这时问题就变成了多元回归问题)，例如，二次曲线y=ax2+bx+c。这时其实只要像下面这样简单地改变一下矩阵的形式即可：
X=x21x11
x22x21
︙︙︙
x2nxn1,w=a
b
c,y=y1
y2
︙
yn
事实上，可以看出对于多元回归问题，最终要面对的问题都是要求一个‖Xw-y‖2的最佳逼近，采用线性代数的方法，就将一个相当复杂的问题统一到一个相对简单的形式上来。
假设有一个很大的方程组Xw=y，如果方程组的解不存在，但又需要对它进行求解时，其实要做的就是去寻找一个最佳的w，使得Xw尽量接近y。
考虑Xw作为y的一个近似，那么y到Xw的距离越小，以‖Xw-y‖来度量的近似程度也就越好。一般的最小二乘问题就是要找出使得‖Xw-y‖尽量小的w，术语“最小二乘”来源于这样的事实： ‖Xw-y‖是平方和的平方根。
定义： 如果m×n的矩阵X和向量y∈Rm，Xw-y的最小二乘解是w^∈Rn，使得‖Xw^-y‖≤‖Xw-y‖对所有的w∈Rn成立。
回忆一下列空间的概念。假设大小为m×n的矩阵X=［x1,x2,…,xn］，此处，xi表示矩阵中的一列，其中1≤i≤n，那么矩阵X的列空间(记为ColX)就是由X中各列的所有线性组合构成的集合，即ColX=Span{x1,x2,…,xn}。
很显然，ColX=Span{x1,x2,…,xn}是一个子空间，进而可以知道m×n的矩阵X的列空间是Rm的一个子空间。此外，我们还注意到ColX中的一个典型向量(或者称为元素)可以写成Xw的形式，其中w是某个向量，记号Xw表示X的列向量的一个线性组合，也就是说，ColX={y：y=Xw,w∈Rn)，于是该式告诉我们： 
（1) 记号Xw代表ColX空间中的向量； 
（2) ColX是线性变换w Xw的值域。
假设W是Rn空间中的一个子空间，y是Rn中的任意向量，y^是y在W上的正交投影，那么y^是W中最接近y的点，也就是说，‖y-y^‖≤‖y-�瘙
經‖对所有属于W又异于y^的�瘙
經成立。这也就是本节开始时介绍的最佳逼近原理。其中，向量y^称为W中元素对y的最佳逼近。还可以知道，对于给定元素y，可以被某个给定子空间中的元素代替或“逼近”，用‖y-�瘙
經‖表示的从y到�瘙
經的距离，可以认为是用�瘙
經代替y的“误差”，而最佳逼近原理说明误差在y^=�瘙
經处取得最小值。




图36最小二乘法的几何解释


最佳逼近原理看似抽象，其实它的几何意义也是相当直观的。如图36所示，所取的向量空间为R2，其中的一个子空间W就是横轴，y是R2中的任意向量，y^是y在W上的正交投影，那么y^显然是W中最接近y的点。更重要的是，这个结论不仅仅在欧几里得几何空间中成立，因为泛函中所讨论的“向量”可能是我们通常所说的向量，但也可能是函数、矩阵，甚至多项式等，但最佳逼近原理皆成立。
通过前面的分析，易知，最小二乘问题的最重要特征是无论怎么选取w，向量Xw必然位于列空间ColX中。而最小二乘问题的本质任务就是要寻找一个w，使得Xw是ColX中最接近y的向量。
还可以在复数域中把最小二乘问题重述如下： 给定X∈Mm×n(F)，以及y∈Fm，找到一个w^∈Fn，使得‖Xw^-y‖最小(特别地，对于回归问题而言，就是给定一个测度的集合(xi,yi)，其中 i=1,2,…,n，找到一个最佳的n-1阶多项式，来表示这个测度的集合)。
与之前的情况一样，只要利用“最佳逼近原理”，问题即可迎刃而解！如图37所示，子空间W是X的值域，即W中的每一个元素


图37最小二乘的几何解释




都是Xw。y是Fm空间中的一个元素，如果要在W中找一个使得‖Xw-y‖最小的w^，显然w^就是y在W中的投影。

特别地，如果rank(X)=n，则有唯一解w^=(X*X)-1X*y。这与本书前面讨论的最小二乘解相一致，但此时我们把结论拓展到一个更大的范围内(复数域)给出结论。
注意到，令wi=1,2,…,m都是互不相同的，假设m≥n，则有rank(X)=n。
下面来推导这个结论。我们声称w^是最小二乘问题的解，意思就是需要证明，对于任何的w∈Fn，有‖Xw^-y‖≤‖Xw-y‖成立。
证明： ‖Xw-y‖2=‖Xw^-y+Xw-Xw^‖2，根据勾股定理可得，‖Xw^-y+Xw-Xw^‖2=‖Xw^-y‖2+‖Xw-Xw^‖2≥‖Xw^-y‖2 。
对于所有的w∈Fn，
y-Xw^⊥W=R(X)〈Xw,y-Xw^〉=0

〈w,X*(y-Xw^)〉=0
因为对于所有的w∈Fn上式都成立，所以有
X*(y-Xw^)=0 X*Xw^=X*y
因为X*的大小是n×m，X的大小是m×n，所以X*X∈Mn×n(F)，如果rank(X*X)=n，那么X*X是可逆的，则有
w^=(X*X)-1X*y
上述过程中的几个步骤可以实现，由下面两个引理来保证。
引理1： 〈Xw,y〉m=〈w,X*y〉n，这里X∈Mm×n(F)，w∈Fn以及y∈Fm。
引理2： rank(X*X)=rank(X)，X∈Mm×n(F)。
综上所述，结论得证。
3.3.2最小二乘问题的计算
对于给定的X和y，运用最佳逼近定理于ColX空间，取
y^=projColXy
由于y属于X的列空间，方程Xw=y^是相容的且存在一个属于Rn的w^使得Xw^=y^。由于y^是ColX中最接近y的点，一个向量w^是Xw=y的一个最小二乘解的充分必要条件是w^满足Xw^=y^。这个属于Rn的w^是一系列由X构造y^的权值(显然如果方程Xw^=y^中存在自由变量，则这个方程会有多个解)。
如果w^满足Xw^=y^，根据正交分解定理，投影y^具有这样的一个性质： (y-y^)与ColX正交，即(y-Xw^)正交于X的每一列，如果xj是X的任意列，那么xj·(y-Xw^)=0且xTj(y-Xw^)=0(注意前者是内积运算，后者是矩阵乘法)，由于每个xTj是XT的一行，于是可得
XT(Xw^-y)=0
因此有
XTXw^-XTy=0
或者
XTXw^=XTy
这其实表明Xw=y的每个最小二乘解满足方程
XTXw=XTy
上述矩阵方程表示的线性方程组常称为Xw=y的法方程，其解通常用w^表示。最后可以得出如下定理(证明从略)。
定理： 方程Xw=y的最小二乘解集和法方程XTXw=XTy的非空解集一致。
下面通过一个例子来说明最小二乘问题的解法。
例： 求不相容方程Xw=y的最小二乘解，其中
X=40
02
11,y=2
0
11
利用前面给出的公式计算
XTX=401
02140
02
11=171
15

XTy=401
0212
0
11=19
11
那么法方程XTXw=XTy就变成了
171
15w1
w2=19
11
行变换可用于解此方程组，但由于XTX是2×2的可逆矩阵，于是很容易得到
(XTX)-1=1845-1
-117
那么可以解XTXw=XTy如下： 
w^=(XTX)-1XTy=1845-1
-11719
11=1
2
在很多计算中，XT是可逆的，但也并非都是如此，下面例子中的矩阵常常出现在统计学中的方差分析问题里。
例： 求Xw=y的最小二乘解，其中
X=1100
1100
1010
1010
1001
1001,y=-3
-1
0
2
5
1
类似地，可以算得
XTX=111111
110000
001100
0000111100
1100
1010
1010
1001
1001=6222
2200
2020
2002
XTy=111111
110000
001100
000011-3
-1
0
2
5
1=4
-4
2
6
矩阵方程XTXw=XTy的增广矩阵为 
62224
2200-4
20202
20026~10013
010-1-5
001-1-2
00000
通解是： x1=3-x4，x2=-5+x4，x3=-2+x4，x4是自由变量。
所以，Xw=y的最小二乘通解具有如下形式： 
w^=3
-5
-2
0+x4-1
1
1
1
在什么条件下，方程Xw=y的最小二乘解是唯一的？下面的定理给出了判断准则(当然，正交投影总是唯一的)，我们不具体讨论该定理的证明。
定理： 矩阵XTX是可逆的，其充分必要条件是X的列是线性无关的，在这种情况下，方程Xw=y有唯一最小二乘解w^且它有下面的表达式
w^=(XTX)-1XTy
注意，这是一个通解，也就是多元线性回归问题的解都可以用它来表示。每种具体的线性回归问题的解都可以看成是它的特例，例如第2章中给出的一元线性回归的解是
w^0=∑x2i∑yi-∑xi∑xiyin∑x2i-∑xi2

w^1=n∑xiyi-∑xi∑yin∑x2i-∑xi2
下面就基于线性回归的通解来推导证明上述关于一元线性回归最小二乘解的结论。对于一元线性回归而言，假设训练集中有n个数据：
y1
y2
︙
yn=1x1
1x2
︙︙
1xnw0
w1
代入通解公式w^=(XTX)-1XTy，则有
w0
w1=11…1
x1x2…xn1x1
1x2
︙︙
1xn-111…1
x1x2…xny1
y2
︙
yn
其中，
11…1
x1x2…xny1
y2
︙
yn=∑yi
∑xiyi
另外，
11…1
x1x2…xn1x1
1x2
︙︙
1xn=n∑xi
∑xi∑x2i
因此，
(XTX)-1=1n∑x2i-∑xi2∑x2i-∑xi
-∑xin
综上可得
w^=(XTX)-1XTy

w^0
w^1=1n∑x2i-∑xi2∑x2i-∑xi
-∑xin∑yi
∑xiyi
=1n∑x2i-∑xi2∑x2i∑yi-∑xi∑xiyi
n∑xiyi-∑xi∑yi
最终证明
w^0=∑x2i∑yi-∑xi∑xiyin∑x2i-∑xi2
w^1=n∑xiyi-∑xi∑yin∑x2i-∑xi2
下面的例子表明，当X的列向量不正交时，该如何找到Xw=y的最小二乘解。这类矩阵在线性回归中常被用到。
例： 找出Xw=y的最小二乘解，其中
X=1-6
1-2
11
17,y=-1
2
1
6
解： 由于X的列x1和x2相互正交，y在ColX的正交投影如下： 
y^=y·x1x1·x1·x1+y·x2x2·x2·x2=84·x1+4590·x2=2
2
2
2+-3
-1
1/2
7/2=-1
1
5/2
11/2
既然y^已知，我们可以解Xw^=y^。这个很容易，因为已经知道y^用X的列线性(通过线性组合来)表示时的权值。于是从上式可以立刻得到
w^=8/4
45/90=2
1/2
某些时候，最小二乘解问题的法方程可能是病态的，也就是XTX中元素在计算中出现较小的误差，可导致解w^出现较大的误差。如果X的列线性无关，最小二乘解常常可通过X的QR分解更可靠地求出。来看下面这个定理。
定理： 给定一个m×n的矩阵X，且具有线性无关的列，取X=QR是X的QR分解，那么对每一个属于Rm的y，方程Xw=y有唯一的最小二乘解，其解为w^=R-1QTy。
这个定理的证明非常简单，这里不再赘述。更进一步，基于QR分解的知识，我们知道Q的列形成ColX的正交基，因此，QQTy是y在ColX上的正交投影y^，那么Xw^=y^说明w^是Xw=y的最小二乘解。
此外，由于上述定理中的R是上三角形矩阵，w^可从方程Rw=QTy计算得到。求解该方程时，通过回代过程或行变换会比较高效。
例： 求出Xw=y的最小二乘解，其中
X=135
110
112
133,y=3
5
7
-3
解： 可以计算矩阵X的QR分解为
X=QR=1/21/21/2
1/2-1/2-1/2
1/2-1/21/2
1/21/2-1/2245
023
002
那么
QTy=1/21/21/21/2
1/2-1/2-1/21/2
1/2-1/21/21/23
5
7
-3=6
-6
4
满足Rw=QTy的最小二乘解是w^，也就是
245
023
002x1
x2
x3=6
-6
4
这个方程很容易解出得
w^=10
-6
2
3.4多元回归模型检验
借由最小二乘法所构建的线性回归模型是否给出了观察值的一种有效描述呢？或者说，所构建的模型是否具有一定的解释力？要回答这些问题，就需要对模型进行一定的检验。
3.4.1线性回归的显著性
与一元线性回归类似，要检测随机变量y和可控变量x1,x2,…,xm之间是否存在线性相关关系，即检验关系式y=w0+w1x1+…+wmxm+u是否成立，其中u~N(0,σ2)。此时主要检验m个系数w1,w2,…,wm是否全为零。如果全为零，则可认为线性回归不显著； 若系数w1,w2,…,wm不全为零，则可认为线性回归是显著的。为进行线性回归的显著性检验，在上述模型中提出原假设和备择假设：
H0： w1=w2=…=wm=0

H1： H0是错误的
设对(x1,x2,…,xm,y)已经进行了n次独立观测，得观测值(xi1,xi2,…,xim,yi)，其中i=1,2,…,n。由观测值确定的线性回归方程为
y^=w^0+w^1x1+…+w^mxm
将(x1,x2,…,xm)的观测值代入，有
y^i=w^0+w^1xi1+…+w^mxim
令
y-=1n∑ni=1yi
采用F分布检验法。首先对总离差平方和进行分解：
SStotal=∑(yi- y-)2=∑(yi- y^i)2+∑(y^i- y-)2
与前面在一元线性回归时讨论的一样，残差平方和
SSresidual=∑(yi- y^i)2
反映了试验时随机误差的影响。
回归平方和
SSregression=∑(y^i- y-)2
反映了线性回归引起的误差。
在原假设成立的条件下，可得
yi=w0+ui,i=1,2,…,n

y-=w0+ u-
观察SSregression和SSresidual的表达式，易见如果SSregression比SSresidual大得多，就不能认为所有的w1,w2,…,wm全为零，即拒绝原假设； 反之则接受原假设。从而考虑由这两项之比构造的检验统计量。
由F分布的定义可知
F=SSregression/mSSresidual/(n-m-1)~F(m,n-m-1)
给定显著水平α，由F分布表查得临界值Fα(m,n-m-1)，使得
P{F≥Fα(m,n-m-1)}=α
由采样得到观测数据，求得F统计量的数值，如果F≥Fα(m,n-m-1)，则拒绝原假设，即线性回归是显著的； 如果F<Fα(m,n-m-1)，则接受原假设，即认为线性回归方程不显著。
在土壤沉淀物吸收情况的例子中，可以算得F统计量的大小为92.03，这个值要远远大于F0.05(2,10)，所以有理由拒绝原假设，即证明回归总体是显著线性的。也可以通过与这个F统计量对应的P值来判断，在图32所示的输出结果中已经给出了该值。或者也可以调用训练所得模型的f_pvalue属性来输出该值： 



print(model_multilin.f_pvalue)




输出结果为3.63e-07，这与图32中给出的结果是一致的。另外，基于对该值计算原理的认识，也可以在Python中使用下面的代码得到相应的P值： 



from scipy.stats import f

print(f.sf(92.03, 2, 10))




可见，P值远远小于0.05，因此有足够的把握拒绝原假设，并同样得到回归总体具有显著线性的结论。
3.4.2回归系数的显著性
在多元线性回归中，若线性回归显著，回归系数不全为零，则回归方程
y^=w^0+w^1x1+w^2x2+…+w^mxm
是有意义的。但线性回归显著并不能保证每一个回归系数都足够大，或者说不能保证每一个回归系数都显著地不等于零。若某一系数等于零，如wj=0，则变量xj对y的取值就不起作用。因此，要考查每一个自变量xj对y的取值是否起作用，其中j=1,2,…，m，就需对每一个回归系数wj进行检验。为此在线性回归模型上提出原假设：
H0： wj=0,1≤j≤m
由于w^j是wj的无偏估计量，自然由w^j构造检验用的统计量。由
w^=(XTX)-1XTy
易知w^j是相互独立的正态随机变量y1,y2,…,yn的线性组合，所以w^j也服从正态分布，并且有
E(w^j)=wj,var(w^j)=σ2A-1jj
即w^j~N(wj,σ2A-1jj)，其中A-1jj是矩阵A-1的主对角线上的第j个元素，而且这里的j是从第0个算起的。于是
w^j-wjσA-1jj~N(0,1)
而
SSresidualσ2~χ2n-m-1
还可以证明w^j与SSresidual是相互独立的。因此在原假设成立的条件下，有
T=w^jA-1jjSSresidual/(n-m-1)~t(n-m-1)
给定显著水平α，查t分布表得到临界值tα/2(n-m-1)，由样本值算得T统计量的数值，若T≥tα/2(n-m-1)则拒绝原假设，即认为wj和零有显著的差异； 若|T|<tα/2(n-m-1),则接受原假设，即认为wj显著地等于零。
由于E［SSresidual/σ2］=n-m-1，所以
σ^*2=SSresidualn-m-1
是σ2的无偏估计。于是T统计量的表达式也可以简写为
T=w^jσ^*A-1jj=w^jse(w^j)~t(n-m-1)
在土壤沉淀物吸收情况的例子中，Python计算得到的各参数估计值为
w^0=-7.35066,w^1=0.11273,w^2=0.34900
于是同3个参数相对应的T统计量分别为
Tw^0=-7.350664.379×0.633138=-2.109

Tw^1=0.112734.379×0.000046=3.797

Tw^2=0.349004.379×0.000265=4.894
相应的P值可以在Python中用下列代码算得：



from scipy.stats import t

print(2*(t.sf(2.109,10)))

print(2*(t.sf(3.797,10)))

print(2*(t.sf(4.894,10)))




输出结果如下： 

0.06114493

0.00350298

0.00062871

或者，也可以直接调用训练所得模型的pvalues属性来输出该值： 



print(model_multilin.pvalues)




输出结果如下： 

［0.06110081 0.00350359 0.00062836］

以上这些与图32中输出的结果是一致的(注意，由于计算精度的问题，最后几位有效数字上的差异是正常的)，且可据此推断回归系数w1和w2是显著(不为零)的。注意,截距项w0是否为零并不是我们需要关心的。
3.5多元线性回归模型预测
对于线性回归模型
y=w0+w1x1+w2x2+…+wmxm+u
其中u~N(0,σ2)，当求得参数w的最小二乘估计w^之后，就可以建立回归方程
y^=w^0+w^1x1+w^2x2+…+w^mxm
而且在经过线性回归显著性及回归系数显著性的检验后，表明回归方程和回归系数都是显著的，那么就可以利用回归方程来进行预测。给定自变量x1,x2,…,xm的任意一组观察值x01,x02,…,x0m，由回归方程可得
y^0=w^0+w^1x01+w^2x02+…+w^mx0m
设x0=(1,x01,x02,…,x0m)，则上式可以写成
y^0=x0w^
在x=x0时，由样本回归方程计算的y^0是个别值y0和总体均值E(y0)的无偏估计，所以y^0可以作为y0和E(y0)的预测值。
与第2章中讨论的情况相同，区间预测包括两个方面： 一方面是总体个别值y0的区间预测； 另一方面是总体均值E(y0)的区间预测。设e0=y0- y^0=y0-x0w^，则有
e0~N(0,σ2［1+x0(xTx)-1xT0］)
如果w^是统计模型中某个参数w的估计值，那么T统计量的定义式就为
tβ^=w^se(w^)
所以与e0相对应的T统计量的表达式如下：
T=e0σ^1+x0(xTx)-1xT0=y0- y^0σ^1+x0(xTx)-1xT0~t(n-m-1)
在给定显著水平α的情况下，可得
y^0-tα2(n-m-1)×σ^1+x0(xTx)-1xT0≤y0

≤y^0+tα2(n-m-1)×σ^1+x0(xTx)-1xT0
总体个别值y0的区间预测就由上式给出。
在Python中如果使用statsmodels包来构建回归模型，那么调用predict()方法可以根据新的输入来进行回归预测。针对土壤沉淀物吸收的例子中，可以使用下面的命令来预测当可提取的铁含量为150，可提取的铝含量为40时，磷酸盐的吸收情况为： 



print(model_multilin.predict(［1, 150, 40］))




如果希望得到预测结果的区间估计，则需要用到get_prediction()方法： 



predictions = model_multilin.get_prediction(［1, 150, 40］)

predictions.summary_frame(alpha=0.05)




执行上述代码，程序输出如图38所示，其中被圈出的部分就是区间估计结果。




图38带区间估计的结果预测


其中点预测的结果是23.51929，在5%的显著水平下，个别值的区间预测结果是(13.33372,33.70486)。
当然也可以根据公式来尝试手动计算一下这个结果，其中的T统计量临界值可以由下面的代码求得：



t. isf(0.05/2, 10)




所得结果为2.228139。
继续计算，为了得到更精确的误差标准值，使用下面的代码进行计算。Python中自动输出的结果4.379是我们所计算结果在保留4位有效数字后得到的。



import math

math.sqrt(sum(model_multilin.resid*model_multilin.resid)/10)#4.379375




另外还可以算得x0(xTx)-1xT0的值为
1
150
40T×0.633138-0.0038260.002477
-0.0038260.000046-0.000088
0.002477-0.0000880.000265×1
150
40=0.08958766
然后在Python中使用下面的代码计算最终的预测区间，可见结果与前面给出的结果是基本一致的。



23.51929 - 2.228139 * math.sqrt(1+0.08958766) * 4.379375#13.33372

23.51929 + 2.228139 * math.sqrt(1+0.08958766) * 4.379375      #33.70486




类似地，还可以得到总体均值E(y0)的区间预测表达式为
y^0-tα2(n-m-1)×σ^x0(xTx)-1xT0≤E(y0)≤y^0+tα2(n-m-1)×σ^x0(xTx)-1xT0
前面由get_prediction()方法得到的输出中，mean_ci_lower和mean_ci_upper给出的就是总体均值的区间预测结果。而且y0期望的置信区间要比y0的置信区间更窄。
同样，下面的代码给出了包含中间过程的手动计算方法，这与刚刚得到的计算结果是一致的。



23.51929 - 2.228139 * math.sqrt(0.08958766) * 4.379375#20.59865

23.51929 + 2.228139 * math.sqrt(0.08958766) * 4.379375        #26.43993




3.6格兰杰因果关系检验
所谓因果关系，可以通过变量之间的依赖性来定义，即作为结果的变量是由作为原因的变量所决定的，原因变量的变化引起结果变量的变化。
因果关系不同于相关关系，而且从一个回归关系式我们并不能确定变量之间是否具有因果关系。虽然我们说回归方程中解释变量是被解释变量的原因，但是，这一因果关系通常是先验设定的，或者是在回归之前就已确定。
实际上，在许多情况下，变量之间的因果关系并不总像农作物产量和降雨量之间的关系那样一目了然，或者没有充分的知识使我们认清变量之间的因果关系。此外，即使某一经济理论宣称某两个变量之间存在一种因果关系，也需要给予经验上的支持。
诺贝尔经济学奖获得者、英国经济学家克莱夫·格兰杰(Clive Granger)，是著名的经济时间序列分析大师，被认为是世界上最伟大的计量经济学家之一。格兰杰从预测的角度给出了因果关系的一种描述性定义，这就是我们现在所熟知的格兰杰因果关系。
格兰杰指出，如果一个变量x无助于预测另一个变量y，则说x不是y的原因。相反，若x是y的原因，则必须满足两个条件：  第一，x应该有助于预测y，即在y关于y的过去值的回归中，添加x的过去值作为独立变量应当显著地增加回归的解释能力； 第二，y不应当有助于预测x，其原因是，如果x有助于预测y，y也有助于预测x，则很可能存在一个或几个其他变量，它们既是引起x变化的原因，也是引起y变化的原因。现在人们一般把这种从预测的角度定义的因果关系称为格兰杰因果关系。
变量x是否为变量y的格兰杰原因，是可以检验的。检验x是否为引起y变化的格兰杰原因的过程如下： 
第一步，检验原假设“H0： x不是引起y变化的格兰杰原因”。首先，估计下列两个回归模型。
无约束回归模型(u)： yt=α0+∑pi=1αiyt-i+∑qi=1βixt-i+εt
有约束回归模型(r)： yt=α0+∑pi=1αiyt-i+εt

式中，α0表示常数项； p和q分别为变量y和x的最大滞后期数，通常可以取的稍大一些； εt为白噪声。
然后，用这两个回归模型的残差平方和RSSu与RSSr构造F统计量： 
F=(RSSr-RSSu)/qRSSu/(n-p-q-1)~F(q,n-p-q-1)
其中，n为样本容量。
检验原假设“H0： x不是引起y变化的格兰杰原因”(等价于检验H0： β1=β2=…=βq=0)是否成立。如果F≥Fα(q,n-p-q-1)，则β1,β2,…,βq显著不为0，应拒绝原假设“H0： x不是引起y变化的格兰杰原因”； 反之，则不能拒绝原假设。
第二步，将y与x的位置交换，按同样的方法检验原假设“H0： y不是引起x变化的格兰杰原因”。
第三步，要得到“x是y的格兰杰原因”的结论，必须同时拒绝原假设“H0： x不是引起y变化的格兰杰原因”和接受原假设“H0： y不是引起x变化的格兰杰原因”。
最后，给出一个在Python中进行格兰杰因果关系检验的实例。从一般的认识来讲，消费与经济增长之间存在相互促进的作用。但是，相比之下，二者中哪一个对另外一个有更强的促进作用在各国经济发展过程中则呈现出不同的结论。就中国而言，自改革开放以来，我们都认为经济增长对消费的促进作用要大于消费对经济增长的促进作用。或者说，在我国经济增长可以作为消费的格兰杰原因，反之不成立。但这一结论是否能够得到计量经济研究的支持呢？下面就在Python中通过对1978—2002年的统计数据(来源为国家统计局网站)展开分析，进而来回答这个问题。
首先将GDP和消费数据存储在两个数组对象中，然后再拼接到一起。



import numpy as np

import statsmodels.api as sm

from statsmodels.tsa.stattools import grangercausalitytests



GDP = np.array(［3645.2, 4062.6, 4545.6, 4889.5, 5330.5, 5985.6, 7243.8, 9040.7, 

10274.4, 12050.6, 15036.8, 17000.9, 18718.3, 21826.2, 26937.3, 35260.0, 

48108.5, 59810.5, 70142.5, 78060.8, 83024.3, 88479.2, 98000.5, 108068.2, 119095.7］)



consumption = np.array(［1759.1, 2014, 2336.9, 2627.5, 2867.1, 3220.9, 3689.5, 

4627.4, 5293.5, 6047.6, 7532.1, 8778, 9435, 10544.5, 12312.2,  15696.2, 

21446.1, 28072.9, 33660.3, 36626.3, 38821.8, 41914.9, 46987.8, 50708.8, 55076.4］)



X_gdp_cons = np.vstack((GDP, consumption)).transpose()

X_cons_gdp = np.vstack((consumption, GDP)).transpose()




在Python中执行格兰杰因果关系检验，同样可以使用statsmodels包。具体来说需要用到的是grangercausalitytests类。实例化该类时，最重要的参数有两个。首先是一个包含有两列时间序列数据的多维数组对象，后续的格兰杰校验将检验位于第二列的时间序列数据是否是第一列时间序列数据的格兰杰原因。因此，上述代码中构建了两个输入数据对象，二者的区别在于GDP和消费数据所在的位置不同。另外一个参数用于设置最大滞后期数。如果输入的是一个整数m，那么表示要对从1~m的各个最大滞后期数都进行一次格兰杰因果关系校验。如果输入的是［m］，则表示只对m这个最大滞后期数做格兰杰因果关系校验。
首先来检验消费是否是GDP的格兰杰原因： 



gc_res_1 = grangercausalitytests(X_gdp_cons, ［2］)




上述代码的输出结果如下： 

Granger Causality

number of lags (no zero) 2

ssr based F test:F=1.6437,p=0.2210, df_denom=18, df_num=2

ssr based chi2 test:		chi2=4.2005,	p=0.1224, df=2

likelihood ratio test:	chi2=3.8581,	p=0.1453, df=2

parameter F test:			F=1.6437,		p=0.2210, df_denom=18, df_num=2

接下来检验GDP是否是消费的格兰杰原因： 



gc_res_2 = grangercausalitytests(X_cons_gdp, ［2］)




上述代码的输出结果如下： 

Granger Causality

number of lags (no zero) 2

ssr based F test:			F=13.4108,		p=0.0003, df_denom=18, df_num=2

ssr based chi2 test:		chi2=34.2720,	p=0.0000, df=2

likelihood ratio test:	chi2=20.9833,	p=0.0000, df=2

parameter F test: 		F=13.4108,		p=0.0003, df_denom=18, df_num=2

检验结果显示，当原假设为“消费不是引起GDP变化的格兰杰原因”，F测试的P值结果为0.221>0.05，我们无法拒绝原假设； 而当原假设变为“GDP不是引起消费变化的格兰杰原因”时，F测试的P值为0.0003<0.05，我们可以据此拒绝原假设。于是可以证明： GDP是消费的格兰杰原因。
之前计算结果中的P值也可以用下面的Python代码算得(注意当p=2时，观测值的数量n=25-2=23)： 



f.sf(1.6437, 2, 18)#0.2209779592

f.sf(13.411, 2, 18)    #0.0002716636