十、多元线性回归

1.多元线性回归概述

多元回归研究的是因变量YYY与mmm个自变量x1,⋯,xmx_1,\cdots,x_mx1,⋯,xm之间的相互依赖关系的，这里假定YYY是随机变量，而x1,⋯,xmx_1,\cdots,x_mx1,⋯,xm是一般变量（视为非随机的），YYY不能由x1,⋯,xmx_1,\cdots,x_mx1,⋯,xm完全决定，但是可以用x1,⋯,xmx_1,\cdots,x_mx1,⋯,xm的函数进行估计。而多元线性回归中的，指的是关于“参数”的线性，我们一般可以把数据写成这样的线性模型：
{y1=β0+β1x11+⋯+βmx1m+ε1,⋯⋯⋯⋯⋯⋯yn=β0+β1xn1+⋯+βmxnm+εn.\left\{ \begin{array}l y_1=\beta_0+\beta_1x_{11}+\cdots+\beta_mx_{1m}+\varepsilon_1, \\ \cdots \cdots \cdots \cdots \cdots \cdots \\ y_n =\beta_0+\beta_1x_{n1}+\cdots+\beta_mx_{nm}+\varepsilon_n. \end{array} \right. ⎩⎨⎧y1=β0+β1x11+⋯+βmx1m+ε1,⋯⋯⋯⋯⋯⋯yn=β0+β1xn1+⋯+βmxnm+εn.
这就是经典多元线性回归模型，这里εt\varepsilon_tεt是偏差，可以将其视为一个与变量无关的白噪声，满足E(εt)=0,D(εt)=σ2,Cov(εi,εj)=0(i≠j){\rm E}(\varepsilon_t)=0,{\rm D}(\varepsilon_t)=\sigma^2,{\rm Cov}(\varepsilon_i,\varepsilon_j)=0(i\ne j)E(εt)=0,D(εt)=σ2,Cov(εi,εj)=0(i=j)，还可以进一步假定其为正态白噪声，即εt∼N(0,σ2)\varepsilon_t\sim N(0,\sigma^2)εt∼N(0,σ2)相互独立。

大型的式子都可以用向量、矩阵形式表达，多元线性回归模型也不例外，这里进行以下定义：
Y=[y1y2⋮yn]n,β=[β0β1⋮βm]m+1,ε=[ε1ε2⋮εn]n,C=(1m∣X)=[1x11x12⋯x1m1x21x22⋯x2m⋮⋮⋮⋮1xn1xn2⋯xnm]n×m,Y=\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}_n,\beta=\begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_m \end{bmatrix}_{m+1},\varepsilon=\begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix}_n,C=(\boldsymbol 1_m|X)=\begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1m} \\ 1 & x_{21} & x_{22} & \cdots & x_{2m} \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{nm} \end{bmatrix}_{n\times m}, Y=⎣⎢⎢⎢⎡y1y2⋮yn⎦⎥⎥⎥⎤n,β=⎣⎢⎢⎢⎡β0β1⋮βm⎦⎥⎥⎥⎤m+1,ε=⎣⎢⎢⎢⎡ε1ε2⋮εn⎦⎥⎥⎥⎤n,C=(1m∣X)=⎣⎢⎢⎢⎡11⋮1x11x21⋮xn1x12x22⋮xn2⋯⋯⋯x1mx2m⋮xnm⎦⎥⎥⎥⎤n×m,
这样就有Y=Cβ+εY=C\beta +\varepsilonY=Cβ+ε，其中ε∼Nn(0,σ2In)\varepsilon\sim N_n(0,\sigma^2I_n)ε∼Nn(0,σ2In)。要能够估计出m+1m+1m+1个参数β0,⋯,βm\beta_0,\cdots,\beta_mβ0,⋯,βm，就至少要有m+1m+1m+1个样本，也就是n>mn>mn>m；同时我们要求各个自变量之间是不相关的，因为如果存在相关性，则自变量可以由其他自变量表出，所需要的回归系数就可以减少，所以我们的要求是rank(C)=m+1{\rm rank} (C)=m+1rank(C)=m+1。

2.最小二乘法

在给定模型后，我们的目的是估计回归系数β\betaβ，在多元线性回归中常用的方法是最小二乘法。最小二乘法的思想就是让偏差平方和达到最小，偏差平方和即
Q(β)=∑t=1net2=∑t=1n[yt−(β0+β1xt1+⋯+βmxmt)]2=(Y−Cβ)′(Y−Cβ).Q(\beta)=\sum_{t=1}^n e_t^2=\sum_{t=1}^n [y_t-(\beta_0+\beta_1x_{t1}+\cdots+\beta_mx_{mt})]^2=(Y-C\beta)'(Y-C\beta). Q(β)=t=1∑net2=t=1∑n[yt−(β0+β1xt1+⋯+βmxmt)]2=(Y−Cβ)′(Y−Cβ).

如果令YtY_tYt为YYY的第ttt行，CtC_tCt为CCC的第ttt行，则
yt−(β0+β1xt1+⋯+βmxmt)=Yt−Ctβ,[yt−(β0+β1xt1+⋯+βmxmt)]2=(Yt−Ctβ)′(Yt−Ctβ).y_t-(\beta_0+\beta_1x_{t1}+\cdots+\beta_mx_{mt})=Y_t-C_t\beta,\\ [y_t-(\beta_0+\beta_1x_{t1}+\cdots+\beta_mx_{mt})]^2=(Y_t-C_t\beta)'(Y_t-C_t\beta). yt−(β0+β1xt1+⋯+βmxmt)=Yt−Ctβ,[yt−(β0+β1xt1+⋯+βmxmt)]2=(Yt−Ctβ)′(Yt−Ctβ).
所以
(Y−Cβ)′(Y−Cβ)=[(Y1−C1β)′(Y2−C2β)′⋮(Yn−Cnβ)′][(Y1−C1β)(Y2−C2β)⋯(Yn−Cnβ)]=∑t=1n(Yt−Ctβ)′(Yt−Ctβ)=∑t=1net2.\begin{aligned} & (Y-C\beta)'(Y-C\beta) \\ =& \begin{bmatrix} (Y_1-C_1\beta)' \\ (Y_2-C_2\beta)' \\ \vdots \\ (Y_n-C_n\beta)' \end{bmatrix}\begin{bmatrix} (Y_1-C_1\beta) & (Y_2-C_2\beta) & \cdots & (Y_n-C_n\beta) \end{bmatrix}\\ =&\sum_{t=1}^n (Y_t-C_t\beta)'(Y_t-C_t\beta) \\ =&\sum_{t=1}^n e_t^2. \end{aligned} ===(Y−Cβ)′(Y−Cβ)⎣⎢⎢⎢⎡(Y1−C1β)′(Y2−C2β)′⋮(Yn−Cnβ)′⎦⎥⎥⎥⎤[(Y1−C1β)(Y2−C2β)⋯(Yn−Cnβ)]t=1∑n(Yt−Ctβ)′(Yt−Ctβ)t=1∑net2.

要使得Q(β)Q(\beta)Q(β)最小，就有
∂Q(β)∂β=∂(Y−Cβ)′(Y−Cβ)∂β=∂(Y′Y)∂β−∂(β′C′Y)∂β−∂(Y′Cβ)∂β+∂(β′C′Cβ)∂β=0m+1−C′Y−C′Y+2C′Cβ=2(C′Cβ−C′Y)\begin{aligned} \frac{\partial Q(\beta)}{\partial \beta}=& \frac{\partial (Y-C\beta)'(Y-C\beta)}{\partial \beta} \\ =&\frac{\partial(Y'Y)}{\partial \beta}-\frac{\partial(\beta'C'Y)}{\partial \beta}-\frac{\partial(Y'C\beta)}{\partial \beta}+\frac{\partial(\beta'C'C\beta)}{\partial \beta}\\ =&\boldsymbol 0_{m+1}-C'Y-C'Y+2C'C\beta\\ =&2(C'C\beta-C'Y) \end{aligned} ∂β∂Q(β)====∂β∂(Y−Cβ)′(Y−Cβ)∂β∂(Y′Y)−∂β∂(β′C′Y)−∂β∂(Y′Cβ)+∂β∂(β′C′Cβ)0m+1−C′Y−C′Y+2C′Cβ2(C′Cβ−C′Y)

这里用到的矩阵微商公式有：
∂A∂β=0,∂β′y∂β=y,∂y′β∂β=y或∂y′Aβ∂β=A′y,∂β′Aβ∂β=(A+A′)β.\frac{\partial A}{\partial \beta}=0, \\ \frac{\partial \beta'y}{\partial \beta}=y, \\ \frac{\partial y'\beta}{\partial\beta}=y或\frac{\partial y'A\beta}{\partial \beta}=A'y,\\ \frac{\partial \beta'A\beta}{\partial \beta}=(A+A')\beta. ∂β∂A=0,∂β∂β′y=y,∂β∂y′β=y或∂β∂y′Aβ=A′y,∂β∂β′Aβ=(A+A′)β.

所以让导数为0，就得到
β^=(C′C)−1C′Y.\hat \beta=(C'C)^{-1}C'Y. β^=(C′C)−1C′Y.
在刚才的推导中，我们得到了最小二乘法估计的参数值β^=(C′C)−1C′Y\hat\beta=(C'C)^{-1}C'Yβ^=(C′C)−1C′Y，结合我们的多元线性回归模型Y=Cβ+εY=C\beta+\varepsilonY=Cβ+ε，在已知XXX时YYY的预测值就应该是
Y^=Cβ^=C(C′C)−1C′Y=dHY.\hat Y=C\hat \beta=C(C'C)^{-1}C' Y\stackrel {\rm d}=HY. Y^=Cβ^=C(C′C)−1C′Y=dHY.
这里的H=C(C′C)−1C′H=C(C'C)^{-1}C'H=C(C′C)−1C′被称为帽子矩阵，可以看到
H′H=C(C′C)−1C′⋅C(C′C)−1C′=H2=C(C′C)−1C′=H.H'H=C(C'C)^{-1}C'\cdot C(C'C)^{-1}C'=H^2=C(C'C)^{-1}C'=H. H′H=C(C′C)−1C′⋅C(C′C)−1C′=H2=C(C′C)−1C′=H.
并且，此时的偏差为e=ε^=Y−Y^=(I−H)Ye=\hat\varepsilon=Y-\hat Y=(I-H)Ye=ε^=Y−Y^=(I−H)Y，残差平方和就是
Q(β^)=ε^′ε^=Y′(I−H)′(I−H)Y=Y′(I−H)Y=Y′Y−Y′HY=Y′Y−Y′Cβ^.\begin{aligned} Q(\hat\beta)=&\hat\varepsilon'\hat\varepsilon \\ =&Y'(I-H)'(I-H)Y \\ =&Y'(I-H)Y\\ =&Y'Y-Y'HY\\ =&Y'Y-Y'C\hat \beta. \end{aligned} Q(β^)=====ε^′ε^Y′(I−H)′(I−H)YY′(I−H)YY′Y−Y′HYY′Y−Y′Cβ^.

3.最小二乘估计量的性质

最小二乘估计量β^=(C′C)−1C′Y\hat \beta=(C'C)^{-1}C'Yβ^=(C′C)−1C′Y是β\betaβ的最小方差线性无偏估计(BLUE)，对这个词做一些说明。

首先，很明显β^\hat \betaβ^是一个随机变量（因为YYY被假定为随机向量，XXX被视为普通向量，所以CCC应该被视为常数阵，YYY视为随机向量），而线性性，指的就是β^\hat \betaβ^可以被Y1,⋯,YnY_1,\cdots,Y_nY1,⋯,Yn所线性表示，这里的(C′C)−1C′(C'C)^{-1}C'(C′C)−1C′就是其线性表示系数。线性性的好处是，一些对于随机变量数字特征的对于独立随机变量的线性情形是有计算公式的，这样就可以很好地写出β^\hat \betaβ^的表现形式及其数字特征。

其次，无偏性指的是Eβ^=β{\rm E}\hat \beta=\betaEβ^=β，注意这里的β\betaβ是客观存在的，是不随抽样改变的常数，而β^\hat \betaβ^是随着我们抽取样本进行观测，根据样本的观测值算出的统计量，它是一个具有两重性的统计量。而
Y=Cβ+ε,E(Y)=Cβ,Eβ^=(C′C)−1C′E(Y)=(C′C)−1(C′C)β=β.Dβ^=E[(β^−β)(β^−β)′]=E[(C′C)−1C′μμ′C(C′C)−1]=(C′C)−1C′Eμμ′C(C′C)−1=σ2(C′C)−1.Y=C\beta+\varepsilon,\quad {\rm E}(Y)=C\beta, \\ {\rm E}\hat \beta=(C'C)^{-1}C'{\rm E}(Y)=(C'C)^{-1}(C'C)\beta=\beta. \\ \begin{aligned} {\rm D}\hat \beta=&{\rm E}[(\hat\beta-\beta)(\hat\beta-\beta)'] \\ =&{\rm E}[(C'C)^{-1}C'\mu\mu'C(C'C)^{-1}]\\ =&(C'C)^{-1}C'{\rm E}\mu\mu' C(C'C)^{-1}\\ =&\sigma^2(C'C)^{-1}. \end{aligned} Y=Cβ+ε,E(Y)=Cβ,Eβ^=(C′C)−1C′E(Y)=(C′C)−1(C′C)β=β.Dβ^====E[(β^−β)(β^−β)′]E[(C′C)−1C′μμ′C(C′C)−1](C′C)−1C′Eμμ′C(C′C)−1σ2(C′C)−1.
最后，最小方差性，指的是在所有的线性无偏估计中，β^\hat \betaβ^是方差最小的一个。假设还有其他的线性无偏估计量β∗\beta^*β∗，记作[(C′C)−1C′+D]Y=C∗Y=β^+DY[(C'C)^{-1}C'+D]Y=C^*Y=\hat \beta+DY[(C′C)−1C′+D]Y=C∗Y=β^+DY，为了满足无偏性，有
Eβ∗=E[C∗(Xβ+μ)]=E(C∗X)β=β,C∗X=I.{\rm E}\beta^*={\rm E}[C^*(X\beta+\mu)]={\rm E}(C^*X)\beta=\beta,\quad C^*X=I. Eβ∗=E[C∗(Xβ+μ)]=E(C∗X)β=β,C∗X=I.
而C∗X=(X′X)−1X′X+DX=I+DX=IC^*X=(X'X)^{-1}X'X+DX=I+DX=IC∗X=(X′X)−1X′X+DX=I+DX=I，所以DX=0DX=0DX=0。于是
Cov(β∗∣X)=E[(C∗Y−β)(C∗Y−β)′∣X]=[E(C∗μ)(C∗μ)′∣X]=E[(C∗μμ′C∗′)∣X]=[(X′X)−1X′+D]E(μμ′∣X)[X(X′X)−1+D′]=σ2[(X′X)−1X′X(X′X)−1+(X′X)−1(XD)′+DX(X′X)−1+DD′]=σ2(X′X)−1+σ2DD′.\begin{aligned} {\rm Cov}(\beta^*|X)=&{\rm E}[(C^*Y-\beta)(C^*Y-\beta)'|X]\\ =&[{\rm E}(C^*\mu)(C^*\mu)'|X]\\ =&{\rm E}[(C^*\mu\mu'{C^*}')|X]\\ =&[(X'X)^{-1}X'+D]{\rm E}(\mu\mu'|X)[X(X'X)^{-1}+D']\\ =&\sigma^2[(X'X)^{-1}X'X(X'X)^{-1}+(X'X)^{-1}(XD)'+DX(X'X)^{-1}+DD']\\ =&\sigma^2(X'X)^{-1}+\sigma^2DD'. \end{aligned} Cov(β∗∣X)======E[(C∗Y−β)(C∗Y−β)′∣X][E(C∗μ)(C∗μ)′∣X]E[(C∗μμ′C∗′)∣X][(X′X)−1X′+D]E(μμ′∣X)[X(X′X)−1+D′]σ2[(X′X)−1X′X(X′X)−1+(X′X)−1(XD)′+DX(X′X)−1+DD′]σ2(X′X)−1+σ2DD′.
由于DDD是正定的，所以σ2DD′>0\sigma^2DD'>0σ2DD′>0，所以Cov(β∗∣X)>σ2(X′X)−1=Cov(β^∣X){\rm Cov}(\beta^*|X)>\sigma^2(X'X)^{-1}={\rm Cov}(\hat \beta|X)Cov(β∗∣X)>σ2(X′X)−1=Cov(β^∣X)，这就证明了β^\hat \betaβ^的最小方差性。
这就说明β^\hat \betaβ^是最小方差的。

因为β^\hat \betaβ^是β\betaβ的最小方差线性无偏估计，所以对于β\betaβ的线性组合α′β\alpha'\betaα′β，其最小方差线性无偏估计依然是a′β^a'\hat\betaa′β^。实际应用时，我们会得到一组数xt1,⋯,xtmx_{t1},\cdots,x_{tm}xt1,⋯,xtm，要用它预报yty_tyt，实际上就是对β\betaβ分量的线性组合，这也是我们用Cβ^C\hat\betaCβ^作为YYY的预测值的理由。

接下来探讨β^\hat \betaβ^作为一个随机变量时的分布。首先由线性性，β^\hat\betaβ^是YYY的分量的线性组合，而YYY的每一个分量中随机部分都是ε\varepsilonε分量，在我们假定ε\varepsilonε是正态白噪声的情况下，YYY的各个分量应该是相互独立且服从正态分布的，即Y∼Nn(Cβ,σ2In)Y\sim N_n(C\beta,\sigma^2I_n)Y∼Nn(Cβ,σ2In)，所以
E(β^)=(C′C)−1C′E(Y)=β,D(β^)=(C′C)−1C′D(Y)C(C′C)−1=σ2(C′C)−1.{\rm E}(\hat\beta)=(C'C)^{-1}C'{\rm E}(Y)=\beta, \\ {\rm D}(\hat\beta)=(C'C)^{-1}C'{\rm D}(Y)C(C'C)^{-1}=\sigma^2(C'C)^{-1}. E(β^)=(C′C)−1C′E(Y)=β,D(β^)=(C′C)−1C′D(Y)C(C′C)−1=σ2(C′C)−1.
即β^∼Nm+1(β,σ2(C′C)−1)\hat \beta\sim N_{m+1}(\beta,\sigma^2(C'C)^{-1})β^∼Nm+1(β,σ2(C′C)−1)。

最后简要提一下白噪声方差σ2\sigma^2σ2的估计，由最大似然估计法，β\betaβ的极大似然估计依然是β^\hat \betaβ^，还可以得到σ2\sigma^2σ2的最大似然估计是
σ^2=(Y−Cβ^)′(Y−Cβ^)n=Q(β^)n.\hat \sigma^2=\frac{(Y-C\hat\beta)'(Y-C\hat\beta)}{n}=\frac{Q(\hat\beta)}{n}. σ^2=n(Y−Cβ^)′(Y−Cβ^)=nQ(β^).
但这个估计量不是无偏的，所以我们对其作修正，就得到了最小二乘估计
s2=(Y−Cβ^)′(Y−Cβ^)n−m−1,E(s2)=σ2.s^2=\frac{(Y-C\hat\beta)'(Y-C\hat\beta)}{n-m-1},\quad {\rm E}(s^2)=\sigma^2. s2=n−m−1(Y−Cβ^)′(Y−Cβ^),E(s2)=σ2.

回顾总结

多元线性回归模型指的是这样的模型：Y=Cβ+εY=C\beta+\varepsilonY=Cβ+ε，其中C=(1n∣X)C=(\boldsymbol 1_n|X)C=(1n∣X)，εt∼WN(0,σ2)\varepsilon_t\sim {\rm WN}(0,\sigma^2)εt∼WN(0,σ2)。即使εt\varepsilon_tεt不是正态白噪声，也至少应当满足零均值、同方差、序列不相关性。
使用最小二乘法求解多元线性回归模型，得到β\betaβ的最小二乘估计量是β^=(C′C)−1C′Y\hat\beta=(C'C)^{-1}C'Yβ^=(C′C)−1C′Y，基于此估计量和数据CCC估计Y^=Cβ^=C(C′C)−1C′Y=HY\hat Y=C\hat\beta=C(C'C)^{-1}C'Y=HYY^=Cβ^=C(C′C)−1C′Y=HY，将HHH称为帽子矩阵。
β^\hat\betaβ^是β\betaβ的最佳线性无偏估计(BLUE)，即满足线性性、无偏性、最小方差性。且其分布是
β^∼Nm+1(β,σ2(C′C)−1).\hat\beta\sim N_{m+1}(\beta,\sigma^2(C'C)^{-1}). β^∼Nm+1(β,σ2(C′C)−1).
在获得β^\hat\betaβ^后，残差平方和为Q(β^)=(Y−Cβ^)′(Y−Cβ^)=Y′Y−Y′Cβ^Q(\hat\beta)=(Y-C\hat\beta)'(Y-C\hat\beta)=Y'Y-Y'C\hat\betaQ(β^)=(Y−Cβ^)′(Y−Cβ^)=Y′Y−Y′Cβ^，可以用残差平方和估计σ2\sigma^2σ2，其无偏估计为s2=Q(β^)/(n−m−1)s^2=Q(\hat\beta)/(n-m-1)s2=Q(β^)/(n−m−1)。

【多元统计分析】10.多元线性回归相关推荐

python 多元线性回归_多元统计分析之多元线性回归的R语言实现
多元统计分析之多元线性回归的R语言实现多元统计分析--multivariate statistical analysis 研究客观事物中多个变量之间相互依赖的统计规律性.或从数学上说, 如果个体的观 ...
【多元统计分析】一、多元统计分析概述
〇.前情提要 b站看[厦门大学MOOC]多元统计分析,因为老师很好看. 参考: [厦门大学MOOC]多元统计分析 https://www.bilibili.com/video/BV1v7411E7PB ...
SPASS多元统计分析在无人机航线精度方面的综合评价
SPASS多元统计分析在无人机航线精度方面的综合评价无人机作为获取低空航空影像的载体,其航线飞行精度直接决定影像的质量与后期图像数据处理精度.GPS/IMU系统作为无人机目标航点及任务航点导航与定位 ...
spss典型相关分析_R语言实战多元统计分析Day10— —典型相关分析
R语言实战多元统计分析Day10-- 典型相关分析 01 前言典型相关分析是用于分析两组随机变量之间的相关性程度的一种统计方法,它能够有效的揭示两组随机变量之间的相互线性依赖关系,这种方法是由Hot ...
多元统计分析朱建平pdf_应用多元统计分析课后答案朱建平版[精心整理].doc
应用多元统计分析课后答案朱建平版[精心整理] 第二章 2.1.试叙述多元联合分布和边际分布之间的关系. 解:多元联合分布讨论多个随机变量联合到一起的概率分布状况,的联合分布密度函数是一个p维的函数,而 ...
【应用多元统计分析】-王学民Python主成分分析例题，特征值处理和可视化（2）
title: "应用多元统计分析" subtitle: "书上题目" author: | OLSRR 由于字数限制,本文省去部分数据预览. 7.6 下表中给出的 ...
【多元统计分析】均值向量和协方差阵的检验——spss上机实验
均值向量和协方差阵的检验--spss上机实验 #参考书目为<多元统计分析>(第五版)--何晓群.中国人民大学出版社 #如有错误,请指正!谢谢~ #关注公众号搜索同名文章获取数据~ 习题2. ...
【应用多元统计分析】——第三章（1）
应用多元统计分析(2022春) 3.1 设 x ∼ N p ( μ , Σ ) , rank ⁡ ( Σ ) = r , u ∼ N r ( 0 , I ) \boldsymbol{x} \sim N ...
数学建模——多元统计分析例题及程序
多元统计分析例题及程序主成分分析简述基本思想计算步骤例题程序程序理解因子分析相关性分析回归分析一元回归例题聚类分析主成分分析简述主成分分析(Principal Compo ...
多元统计分析最短距离法_(完整word版)多元统计分析模拟试题
多元统计分析模拟试题(两套:每套含填空.判断各二十道) A卷 1)判别分析常用的判别方法有距离判别法.贝叶斯判别法.费歇判别法.逐步判别法. 2)Q型聚类分析是对样品的分类,R型聚类分析是对变量_的 ...

【多元统计分析】10.多元线性回归

文章目录