线性回归模型详解
一、简介 线性回归模型是一种很常见、很常用的模型,下面我们简单介绍一下。二、介绍2.1 模型 线性回归模型的表达式为:Y=Xβ+ε(1)Y = X\beta + \varepsilon\tag1Y=Xβ+ε(1)其中,Y=(y1,y2,⋯ ,yn)TY = (y_1, y_2, \cdots, y_n)^TY=(y1,y2,⋯,yn)T, β=(β0,β1,⋯ ,βp−1)T\beta=(\
零、优缺点
优点:
- 建模速度快,不需要很复杂的计算,在数据量大的情况下依然运行速度很快。
- 可以根据系数给出每个变量的理解和解释
- 当特征数量大于样本数量时,表现更好?
缺点:
- 不能很好地拟合非线性数据
- 对异常值敏感?
不用标准化
一、简介
线性回归模型是一种很常见、很常用的模型,下面我们简单介绍一下。
二、最小二乘估计
2.1 模型
线性回归模型的表达式为: Y = X β + ε (1) Y = X\beta + \varepsilon\tag1 Y=Xβ+ε(1)
其中, Y = ( y 1 , y 2 , ⋯ , y n ) T Y = (y_1, y_2, \cdots, y_n)^T Y=(y1,y2,⋯,yn)T, β = ( β 0 , β 1 , ⋯ , β p − 1 ) T \beta=(\beta_0, \beta_1, \cdots, \beta_{p-1})^T β=(β0,β1,⋯,βp−1)T, ε = ( ε 1 , ε 2 , ⋯ , ε n ) T \varepsilon = (\varepsilon_1, \varepsilon_2, \cdots, \varepsilon_n)^T ε=(ε1,ε2,⋯,εn)T为误差项,表示除了自变量外其它因素对 Y Y Y的影响以及试验和测量误差, X = [ 1 x 11 ⋯ x 1 , p − 1 1 x 21 ⋯ x 2 , p − 1 ⋮ ⋮ ⋮ ⋮ 1 x n 1 ⋯ x n , p − 1 ] n × p X = \begin{bmatrix} 1 & x_{11} & \cdots & x_{1, p-1} \\ 1 & x_{21} & \cdots & x_{2, p-1} \\ \vdots & \vdots & \vdots &\vdots \\ 1 & x_{n1} & \cdots & x_{n, p-1}\\ \end{bmatrix}_{n\times p} X=⎣⎢⎢⎢⎡11⋮1x11x21⋮xn1⋯⋯⋮⋯x1,p−1x2,p−1⋮xn,p−1⎦⎥⎥⎥⎤n×p, n n n为样本数。
上述表达式 ( 1 ) (1) (1)也可以展开来写: y i = β 0 + β 1 x i 1 + β 2 x i 2 + ⋯ + β p − 1 x i , p − 1 + ε i , i = 1 , 2 , ⋯ , n (2) y_i = \beta_0+ \beta_1x_{i1} + \beta_2x_{i2} + \cdots + \beta_{p-1}x_{i,p-1} + \varepsilon_i, i = 1, 2, \cdots, n\tag2 yi=β0+β1xi1+β2xi2+⋯+βp−1xi,p−1+εi,i=1,2,⋯,n(2)
模型假设(Gauss-Markov假设):
- E ( ε i ) = 0 E(\varepsilon_i)=0 E(εi)=0(无系统性偏差)
- V a r ( ε i ) = σ 2 Var(\varepsilon_i)=\sigma^2 Var(εi)=σ2(等方差)
- C o v ( ε i , ε j ) = 0 , i ≠ j Cov(\varepsilon_i, \varepsilon_j)=0, i\ne j Cov(εi,εj)=0,i=j(不相关)
更严格的假设: ε i ∼ N ( 0 , σ 2 ) , i = 1 , 2 , ⋯ , n \varepsilon_i\sim N(0, \sigma^2), i = 1, 2, \cdots, n εi∼N(0,σ2),i=1,2,⋯,n, ϵ ∼ N ( 0 , σ 2 I n ) \epsilon\sim N(0,\sigma^2I_n) ϵ∼N(0,σ2In)
将Gauss-Markov假设写成矩阵形式,并结合 ( 1 ) (1) (1)即可得到最基本、最重要的线性回归模型: Y = X β + ε , E ( ε ) = 0 , C o v ( ε ) = σ 2 I n (*) Y = X\beta + \varepsilon,E(\varepsilon)=0,Cov(\varepsilon)=\sigma^2I_n\tag{*} Y=Xβ+ε,E(ε)=0,Cov(ε)=σ2In(*)
2.2 参数估计(最小二乘法)
这个方法是找 β \beta β的估计,使得偏差向量 Y − X β Y-X\beta Y−Xβ的长度平方达到最小,即等价于最小化函数 Q ( β ) : Q(\beta): Q(β): Q ( β ) = ∑ i = 1 n ( y i − β 0 − β 1 x i 1 − ⋯ − β p − 1 x i , p − 1 ) 2 : = ∣ ∣ Y − X β ∣ ∣ 2 (3) Q(\beta)=\sum_{i=1}^n(y_i-\beta_0-\beta_1x_{i1}-\dots-\beta_{p-1}x_{i,p-1})^2:=||Y-X\beta||^2\tag3 Q(β)=i=1∑n(yi−β0−β1xi1−⋯−βp−1xi,p−1)2:=∣∣Y−Xβ∣∣2(3)
上述式子 ( 3 ) (3) (3)可写为: Q = ( Y − X β ) ⊤ ( Y − X β ) = Y ⊤ Y − 2 Y ⊤ X β + β ⊤ X ⊤ X β (4) Q=(Y-X\beta)^\top(Y-X\beta) = Y^\top Y-2Y^\top X\beta+\beta^\top X^\top X\beta\tag4 Q=(Y−Xβ)⊤(Y−Xβ)=Y⊤Y−2Y⊤Xβ+β⊤X⊤Xβ(4)
对其求导,并令其为0,则可以得到方程组 X ⊤ X β = X ⊤ Y X^\top X\beta=X^\top Y X⊤Xβ=X⊤Y(称为正则方程,有唯一解的充要条件是 X ⊤ X X^\top X X⊤X的秩为 p p p,等价的, X X X的秩为 p p p),进而可得参数 β \beta β的估计为: β ^ = ( X ⊤ X ) − 1 X ⊤ Y (5) \hat\beta = (X^\top X)^{-1}X^\top Y\tag5 β^=(X⊤X)−1X⊤Y(5)
当误差项 ϵ ∼ N ( 0 , σ 2 I n ) \epsilon\sim N(0,\sigma^2I_n) ϵ∼N(0,σ2In)时,最小二乘估计得到的结果和极大似然估计得到的结果相同。
2.3 期望和方差
在Gauss-Markov假设下,我们可以得到 β ^ \hat\beta β^的期望和方差为:
- E [ β ^ ] = β E[\hat\beta]=\beta E[β^]=β
- V a r [ β ^ ] = σ 2 ( X ⊤ X ) − 1 Var[\hat\beta]=\sigma^2(X^\top X)^{-1} Var[β^]=σ2(X⊤X)−1
证明:
E [ β ^ ] = E [ ( X ⊤ X ) − 1 X ⊤ Y ] = ( X ⊤ X ) − 1 X ⊤ E [ Y ] = ( X ⊤ X ) − 1 X ⊤ ( X β ) = β E[\hat\beta]= E[(X^\top X)^{-1}X^{\top}Y]= (X^\top X)^{-1}X^{\top}E[Y]=(X^\top X)^{-1}X^{\top}(X\beta)=\beta E[β^]=E[(X⊤X)−1X⊤Y]=(X⊤X)−1X⊤E[Y]=(X⊤X)−1X⊤(Xβ)=β
V a r [ β ^ ] = V a r [ ( X ⊤ X ) − 1 X ⊤ Y ] = ( X ⊤ X ) − 1 X ⊤ V a r ( Y ) X ( X ⊤ X ) − 1 = σ 2 ( X ⊤ X ) − 1 Var[\hat\beta] = Var[(X^\top X)^{-1}X^{\top}Y]=(X^\top X)^{-1}X^{\top}Var(Y)X(X^\top X)^{-1}=\sigma^2(X^\top X)^{-1} Var[β^]=Var[(X⊤X)−1X⊤Y]=(X⊤X)−1X⊤Var(Y)X(X⊤X)−1=σ2(X⊤X)−1
这个定理的第一条表明最小二乘估计 β ^ \hat\beta β^是 β \beta β的无偏估计,也就是说估计没有系统性偏差
2.4 Gauss-Markov定理
对于线性回归模型 ∗ * ∗,在 c ⊤ β c^\top\beta c⊤β( c c c为 p × 1 p\times1 p×1常数向量)的所有线性无偏估计中,最小二乘估计 c ⊤ β ^ c^\top\hat\beta c⊤β^是唯一具有最小方差的估计。
证明:
设 a ⊤ Y a^\top Y a⊤Y是 c ⊤ β c^\top\beta c⊤β的任意一个线性无偏估计,于是有 c ⊤ β = E ( a ⊤ Y ) = a ⊤ X β c^\top\beta=E(a^\top Y)=a^\top X\beta c⊤β=E(a⊤Y)=a⊤Xβ
所以 a ⊤ X = c ⊤ a^\top X=c^\top a⊤X=c⊤恒成立。
对于最小二乘估计 c ⊤ β ^ = c ⊤ ( X ⊤ X ) − 1 X ⊤ Y c^\top\hat\beta=c^\top(X^\top X)^{-1}X^\top Y c⊤β^=c⊤(X⊤X)−1X⊤Y,由2.3里的性质2可知, V a r ( c ⊤ β ^ ) = σ 2 c ⊤ ( X ⊤ X ) − 1 c Var(c^\top\hat\beta)=\sigma^2c^\top(X^\top X)^{-1}c Var(c⊤β^)=σ2c⊤(X⊤X)−1c,而 V a r ( a ⊤ Y ) = σ 2 a ⊤ a Var(a^\top Y)=\sigma^2a^\top a Var(a⊤Y)=σ2a⊤a,所以只需证明 a ⊤ a ≥ c ⊤ ( X ⊤ X ) − 1 c a^\top a\ge c^\top(X^\top X)^{-1}c a⊤a≥c⊤(X⊤X)−1c即可,结合 a ⊤ X = c ⊤ a^\top X=c^\top a⊤X=c⊤,只需证明 a ⊤ a ≥ a ⊤ X ( X ⊤ X ) − 1 X ⊤ a ⟺ a ⊤ ( I − X ( X ⊤ X ) − 1 X ⊤ ) a ≥ 0 a^\top a\ge a^\top X(X^\top X)^{-1}X^\top a\iff a^\top(I-X(X^\top X)^{-1}X^\top)a\ge0 a⊤a≥a⊤X(X⊤X)−1X⊤a⟺a⊤(I−X(X⊤X)−1X⊤)a≥0即可
而 I − X ( X ⊤ X ) − 1 X ⊤ I-X(X^\top X)^{-1}X^\top I−X(X⊤X)−1X⊤为对称幂等阵,所以上式恒成立。唯一性也易证。
这个定理说明最小二乘估计 c ⊤ β ^ c^\top\hat\beta c⊤β^在 c ⊤ β c^\top\beta c⊤β的一切线性无偏估计中是最优的,所以也称 c ⊤ β ^ c^\top\hat\beta c⊤β^是 c ⊤ β c^\top\beta c⊤β的最佳线性无偏估计。
2.5 误差项的方差的估计
σ 2 \sigma^2 σ2是模型误差项的方差,反映了模型误差以及观测误差的大小,因此估计它是一个重要的问题。
误差项是 ε = ( ε 1 , ε 2 , ⋯ , ε n ) T \varepsilon =(\varepsilon_1, \varepsilon_2, \cdots, \varepsilon_n)^T ε=(ε1,ε2,⋯,εn)T
不妨设 Y ^ = X β ^ \hat Y = X\hat\beta Y^=Xβ^,那么残差为 ϵ ^ = Y − Y ^ = Y − X β ^ \hat\epsilon = Y-\hat Y=Y-X\hat\beta ϵ^=Y−Y^=Y−Xβ^,可将其看作是残差向量 Y − X β Y-X\beta Y−Xβ的一个估计,那么很自然用残差平方和 S e 2 = Q ( β ^ ) = ∣ ∣ Y − Y ^ ∣ ∣ 2 = ∣ ∣ ϵ ^ ∣ ∣ 2 = ∑ i = 1 n ( y i − x i ⊤ β ^ ) 2 S_e^2=Q(\hat\beta)=||Y-\hat Y||^2=||\hat\epsilon||^2=\sum\limits_{i=1}^n(y_i-x_i^\top\hat\beta)^2 Se2=Q(β^)=∣∣Y−Y^∣∣2=∣∣ϵ^∣∣2=i=1∑n(yi−xi⊤β^)2来衡量 σ 2 \sigma^2 σ2的大小,越小,说明数据和模型拟合的越好。
注意到 β ^ = ( X ⊤ X ) − 1 X ⊤ Y \hat\beta = (X^\top X)^{-1}X^\top Y β^=(X⊤X)−1X⊤Y,所以 Y ^ = X β ^ = X ( X ⊤ X ) − 1 X ⊤ Y = : P Y \hat Y = X\hat\beta=X(X^\top X)^{-1}X^\top Y=:PY Y^=Xβ^=X(X⊤X)−1X⊤Y=:PY,即 P = X ( X ⊤ X ) − 1 X ⊤ P = X(X^\top X)^{-1}X^\top P=X(X⊤X)−1X⊤,而且 P P P有以下性质: P = P ⊤ = P 2 P = P^\top=P^2 P=P⊤=P2 I n − P = ( I n − P ) ⊤ = ( I n − P ) 2 I_n-P = (I_n-P)^\top=(I_n-P)^2 In−P=(In−P)⊤=(In−P)2
所以此时残差向量可写作: ϵ ^ = ( I n − P ) Y \hat\epsilon=(I_n-P)Y ϵ^=(In−P)Y,残差平方和为: S e 2 : = ∣ ∣ ϵ ^ ∣ ∣ 2 = Y ⊤ ( I n − P ) ⊤ ( I n − P ) Y = Y ⊤ ( I n − P ) Y (6) S_e^2 := ||\hat \epsilon||^2 = Y^\top(I_n-P)^\top(I_n-P)Y=Y^\top(I_n-P)Y\tag6 Se2:=∣∣ϵ^∣∣2=Y⊤(In−P)⊤(In−P)Y=Y⊤(In−P)Y(6)
结论: σ ^ 2 = S e 2 n − p \hat\sigma^2 = \frac{S_e^2}{n-p} σ^2=n−pSe2是 σ 2 \sigma^2 σ2的无偏估计。
证明: E [ S e 2 ] = E [ Y ⊤ ( I n − P ) Y ] = E [ t r ( Y ⊤ ( I n − P ) Y ) ] E[S_e^2]= E[Y^\top(I_n-P)Y]=E[tr(Y^\top(I_n-P)Y)] E[Se2]=E[Y⊤(In−P)Y]=E[tr(Y⊤(In−P)Y)]
= E [ t r ( ( I n − P ) Y Y ⊤ ) ] = t r ( ( I n − P ) E [ Y Y ⊤ ] ) = E[tr((I_n-P)YY^\top)]=tr((I_n-P)E[YY^\top]) =E[tr((In−P)YY⊤)]=tr((In−P)E[YY⊤])
= t r ( ( I n − P ) ( V a r [ Y ] + E [ Y ] E [ Y ⊤ ] ) ) =tr((I_n-P)(Var[Y]+E[Y]E[Y^\top])) =tr((In−P)(Var[Y]+E[Y]E[Y⊤]))
= t r ( ( I n − P ) ( σ 2 I n ) ) + t r ( ( I n − P ) X β β ⊤ X ⊤ ) =tr((I_n-P)(\sigma^2 I_n))+tr((I_n-P)X\beta\beta^\top X^\top) =tr((In−P)(σ2In))+tr((In−P)Xββ⊤X⊤)
= σ 2 ( n − t r ( P ) ) =\sigma^2(n-tr(P)) =σ2(n−tr(P))
其中,我们用到了 t r ( A B ) = t r ( B A ) tr(AB)=tr(BA) tr(AB)=tr(BA), ( I n − P ) X = X − X ( X ⊤ X ) − 1 X ⊤ X = 0 (I_n-P)X=X-X(X^\top X)^{-1}X^\top X=\bm0 (In−P)X=X−X(X⊤X)−1X⊤X=0
结合 t r ( P ) = t r ( X ( X ⊤ X ) − 1 X ⊤ ) = t r ( X ⊤ X ( X ⊤ X ) − 1 ) = t r ( I p ) = p tr(P)= tr(X(X^\top X)^{-1}X^\top)=tr(X^\top X(X^\top X)^{-1})=tr(I_p)=p tr(P)=tr(X(X⊤X)−1X⊤)=tr(X⊤X(X⊤X)−1)=tr(Ip)=p,所以 E [ S e 2 ] = ( n − p ) σ 2 E[S_e^2]=(n-p)\sigma^2 E[Se2]=(n−p)σ2,所以 σ ^ 2 = S e 2 n − p \hat\sigma^2 = \frac{S_e^2}{n-p} σ^2=n−pSe2是 σ 2 \sigma^2 σ2的无偏估计
2.6 抽样分布定理
对于线性回归模型 ∗ * ∗,若进一步假设误差 ε i ∼ N ( 0 , σ 2 ) , i = 1 , 2 , ⋯ , n \varepsilon_i\sim N(0, \sigma^2), i = 1, 2, \cdots, n εi∼N(0,σ2),i=1,2,⋯,n,即 ϵ ∼ N ( 0 , σ 2 I n ) \epsilon\sim N(0,\sigma^2I_n) ϵ∼N(0,σ2In),则有如下定理:
- β ^ ∼ N ( β , σ 2 ( X ⊤ X ) − 1 ) \hat\beta \sim N(\beta, \sigma^2(X^\top X)^{-1}) β^∼N(β,σ2(X⊤X)−1)
- ( n − p ) σ ^ 2 σ 2 = S e 2 σ 2 ∼ χ 2 ( n − p ) \frac{(n-p)\hat\sigma^2}{\sigma^2}=\frac{S_e^2}{\sigma^2}\sim \chi^2(n-p) σ2(n−p)σ^2=σ2Se2∼χ2(n−p)
- ϵ ^ \hat\epsilon ϵ^和 Y ^ \hat Y Y^独立
- S e 2 S_e^2 Se2或 σ ^ 2 \hat\sigma^2 σ^2与 β ^ \hat\beta β^独立
证明:
结论1比较容易证明。
结论2的思路:
S e 2 = Y ⊤ ( I n − P ) Y = ( X β + ϵ ) ⊤ ( I n − P ) ( X β + ϵ ) : = ( X β + ϵ ) ⊤ N ( X β + ϵ ) = ϵ ⊤ N ϵ S_e^2=Y^\top(I_n-P)Y=(X\beta+\epsilon)^\top(I_n-P)(X\beta+\epsilon):=(X\beta+\epsilon)^\top N(X\beta+\epsilon)=\epsilon^\top N\epsilon Se2=Y⊤(In−P)Y=(Xβ+ϵ)⊤(In−P)(Xβ+ϵ):=(Xβ+ϵ)⊤N(Xβ+ϵ)=ϵ⊤Nϵ,最后一个等号成立是因为 N X = 0 NX=0 NX=0,再结合 ϵ ∼ N ( 0 , σ 2 I n ) \epsilon\sim N(0,\sigma^2I_n) ϵ∼N(0,σ2In)即可得到。
2.7 回归诊断
我们做的假设是Gauss-Markov假设,但拿到一组数据后,如何判断是否满足Gauss-Markov假设,就需要用到残差分析(残差是误差项的一个估计)详情见回归诊断
2.8 Box-Cox变换
Box-Cox变换是对回归因变量 Y Y Y的如下变换: Y ( λ ) = { 1 λ ( Y λ − 1 ) λ ≠ 0 log Y λ = 0 Y^{(\lambda)} = \begin{cases} \frac 1\lambda (Y^{\lambda}-1) \ \lambda \neq 0\\ \log Y\ \lambda=0 \end{cases} Y(λ)={λ1(Yλ−1) λ=0logY λ=0
其中 λ \lambda λ是一个待定变换参数,所以它是一个变换族,包括了对数变换 ( λ = 0 ) (\lambda=0) (λ=0),平方根变换 λ = 1 / 2 \lambda=1/2 λ=1/2和倒数变换 λ = − 1 \lambda=-1 λ=−1等常用变换。
我们要确定变换参数 λ \lambda λ,使得 Y ( λ ) Y^{(\lambda)} Y(λ)满足 Y ( λ ) = X β + ε , ε ∼ N ( 0 , σ 2 I ) Y^{(\lambda)}=X\beta+\varepsilon, \varepsilon\sim N(0, \sigma^2I) Y(λ)=Xβ+ε,ε∼N(0,σ2I)
即我们要求通过因变量的变换,使得变换过的向量 Y ( λ ) Y^{(\lambda)} Y(λ)与回归自变量具有线性相依关系,误差也服从正态分布,误差各分量是等方差且相互独立。
用极大似然估计来确定 λ \lambda λ。
因为 Y ( λ ) ∼ N ( X β , σ 2 I ) Y^{(\lambda)}\sim N(X\beta,\sigma^2I) Y(λ)∼N(Xβ,σ2I),所以对固定的 λ \lambda λ,可以得到 β \beta β和 σ 2 \sigma^2 σ2的似然函数 L ( β ( λ ) , σ 2 ( λ ) ) L(\beta(\lambda), \sigma^2(\lambda)) L(β(λ),σ2(λ)),进而可以得到 β \beta β和 σ 2 \sigma^2 σ2的极大似然估计 β ^ ( λ ) = ( X ⊤ X ) − 1 X ⊤ Y ( λ ) \hat\beta(\lambda)=(X^\top X)^{-1}X^\top Y^{(\lambda)} β^(λ)=(X⊤X)−1X⊤Y(λ) σ ^ 2 ( λ ) = 1 n Y ( λ ) ( I − X ( X ⊤ X ) − 1 X ⊤ ) Y ( λ ) \hat\sigma^2(\lambda)=\frac{1}{n}Y^{(\lambda)}(I-X(X^\top X)^{-1}X^\top)Y^{(\lambda)} σ^2(λ)=n1Y(λ)(I−X(X⊤X)−1X⊤)Y(λ)
再带入似然函数 L ( β ^ ( λ ) , σ ^ 2 ( λ ) ) L(\hat\beta(\lambda),\hat \sigma^2(\lambda)) L(β^(λ),σ^2(λ)),这是关于 λ \lambda λ的一元函数,通过求它的值来确定 λ \lambda λ
2.9 广义最小二乘估计
误差方差不相等时,即 C o v ( ε ) = σ 2 Σ Cov(\varepsilon)=\sigma^2\bm\Sigma Cov(ε)=σ2Σ,所以此时的回归模型为 Y = X β + ε , E ( ε ) = 0 , C o v ( ε ) = σ 2 Σ (**) Y = X\beta + \varepsilon,E(\varepsilon)=0,Cov(\varepsilon)=\sigma^2\bm\Sigma\tag{**} Y=Xβ+ε,E(ε)=0,Cov(ε)=σ2Σ(**)
为了求参数的估计,需要经过适当变换,把它化成前面讨论过的情形。
既然 Σ \bm\Sigma Σ时正定阵,所以存在 n × n n\times n n×n正交阵 P \bm P P使其对角化 Σ = P ⊤ Λ P \bm\Sigma=\bm P^\top\bm\Lambda\bm P Σ=P⊤ΛP
其中, Λ = d i a g { λ 1 , ⋯ , λ n } \Lambda=diag\{\lambda_1,\cdots,\lambda_n\} Λ=diag{λ1,⋯,λn}, λ i \lambda_i λi是 Σ \bm\Sigma Σ的特征值。
令 Z = Σ − 1 2 Y Z=\bm\Sigma^{-\frac{1}{2}}Y Z=Σ−21Y, U = Σ − 1 2 X U=\bm\Sigma^{-\frac{1}{2}}X U=Σ−21X, e = Σ − 1 2 ε e=\bm\Sigma^{-\frac{1}{2}}\varepsilon e=Σ−21ε,则 C o v ( e ) = σ 2 I Cov(e)=\sigma^2I Cov(e)=σ2I,进而可以得到广义最小二乘估计,也称为Gauss-Markov估计 β ^ ∗ = ( U ⊤ U ) − 1 U ⊤ Z = ( X ⊤ Σ − 1 X ) − 1 X ⊤ Σ − 1 Y \hat\beta^*=(U^\top U)^{-1}U^{\top}Z=(X^\top\bm\Sigma^{-1} X)^{-1}X^{\top}\bm\Sigma^{-1}Y β^∗=(U⊤U)−1U⊤Z=(X⊤Σ−1X)−1X⊤Σ−1Y
2.10 复共线性
对于有些大型线性回归问题,最小二乘估计不总是令人满意。例如,有时某些回归系数的估计值的绝对值异常大,有时回归系数的估计值的符号与问题的实际意义相违背等。原因之一是回归自变量之间存在着近似线性关系,称为复共线性。
2.10.1 均方误差
设 θ \theta θ为 p × 1 p\times1 p×1未知参数向量, θ ~ \tilde \theta θ~是 θ \theta θ的一个估计。定义 θ ~ \tilde\theta θ~的均方误差为 M S E ( θ ~ ) = E ∣ ∣ θ ~ − θ ∣ ∣ 2 = E ( θ ~ − θ ) ⊤ ( θ ~ − θ ) MSE(\tilde\theta)=E||\tilde\theta-\theta||^2=E(\tilde\theta-\theta)^\top(\tilde\theta-\theta) MSE(θ~)=E∣∣θ~−θ∣∣2=E(θ~−θ)⊤(θ~−θ)
它度量了估计 θ ~ \tilde\theta θ~跟未知参数向量 θ \theta θ平均偏离的大小。一个好的估计应该有较小的均方误差。
可以证明,一个估计的均方误差就是由各分量的方差和偏差所决定的,一个好的估计应该由较小的方差和偏差。
2.10.2 复共线性
可以证明,对于最小二乘估计 β ^ \hat\beta β^, M S E ( β ^ ) = σ 2 ∑ i = 1 p 1 λ i MSE(\hat\beta)=\sigma^2\sum\limits_{i=1}^p\frac{1}{\lambda_i} MSE(β^)=σ2i=1∑pλi1
其中, λ i \lambda_i λi为 X ⊤ X X^\top X X⊤X( X X X已经中心化和标准化)的特征根,所以,如果 X ⊤ X X^\top X X⊤X至少有一个特征根非常小,即非常接近 0 0 0,那么 M S E ( β ^ ) MSE(\hat\beta) MSE(β^)就会很大,从均方误差的角度看,这时的最小二乘估计就不是一个好的估计。
而 X ⊤ X X^\top X X⊤X至少有一个特征根非常小意味着 X X X的列之间有近似线性关系,有几个非常小的特征根,就有几个复共线关系。
2.10.3 条件数
定义 k = λ 1 λ p k=\frac{\lambda_1}{\lambda_p} k=λpλ1为方阵 X ⊤ X X^\top X X⊤X的条件数,即 X ⊤ X X^\top X X⊤X的最大特征根与最小特征根之比。
一般来说,若 k < 100 k\lt100 k<100,则认为复共线性的程度很小;若 100 ≤ k ≤ 1000 100\le k\le1000 100≤k≤1000,则认为存在中等强度或较强的复共线性;若 k > 1000 k\gt1000 k>1000,则认为存在严重的复共线性。
2.11 岭估计
当设计阵存在严重的复线性关系时,最小二乘估计并不是一个好的估计,此时可用岭估计,即 β ^ ( k ) = ( X ⊤ X + k I ) − 1 X ⊤ Y \hat\beta(k)=(X^\top X+kI)^{-1}X^\top Y β^(k)=(X⊤X+kI)−1X⊤Y
其中, k > 0 k\gt0 k>0为可选择参数,称为岭参数或偏参数。
岭估计是有偏估计,以牺牲无偏性,换取方差部分的大幅度减少,最终降低其均方误差。
确定 k k k可用 H o e r l − K e n n a r d Hoerl-Kennard Hoerl−Kennard公式和岭迹法(即看 β ^ i ( k ) \hat\beta_i(k) β^i(k)在什么时候趋于稳定)。
三、假设检验与预测
3.1 回归系数的显著性检验
对每个自变量逐一做显著性检验,即对固定的 i , 1 ≤ i ≤ p − 1 i, 1\le i\le p-1 i,1≤i≤p−1,做如下检验 H i : β i = 0 H_i:\beta_i=0 Hi:βi=0
由2.6节可知, β ^ ∼ N ( β , σ 2 ( X ⊤ X ) − 1 ) \hat\beta \sim N(\beta, \sigma^2(X^\top X)^{-1}) β^∼N(β,σ2(X⊤X)−1),不妨令 C = ( X ⊤ X ) − 1 C=(X^\top X)^{-1} C=(X⊤X)−1,其元素为 c i j c_{ij} cij,则有 β ^ i ∼ N ( β i , σ 2 c i i ) \hat\beta_i\sim N(\beta_i,\sigma^2c_{ii}) β^i∼N(βi,σ2cii)
所以当 H i H_i Hi成立时,有 β ^ i σ c i i ∼ N ( 0 , 1 ) \frac{\hat\beta_i}{\sigma\sqrt{c_{ii}}}\sim N(0, 1) σciiβ^i∼N(0,1)
- 当 σ 2 \sigma^2 σ2已知时, β i \beta_i βi的 100 ( 1 − α ) % 100(1-\alpha)\% 100(1−α)%的置信区间为 β ^ i ± u 1 − α / 2 σ c i + 1 , i + 1 , i = 0 , … , p − 1 \hat\beta_i \pm u_{1-\alpha/2}\sigma\sqrt{c_{i+1,i+1}},\ i=0,\dots,p-1 β^i±u1−α/2σci+1,i+1, i=0,…,p−1,所以对于假设 H 0 : β i = 0 v s . H 1 : β i ≠ 0 H_0:\beta_i= 0\ vs.\ H_1:\beta_i\neq 0 H0:βi=0 vs. H1:βi=0,相应的拒绝域为: W = { ∣ β ^ i ∣ > u 1 − α / 2 σ c i + 1 , i + 1 } W = \{|\hat\beta_i|\gt u_{1-\alpha/2}\sigma\sqrt{c_{i+1,i+1}}\} W={∣β^i∣>u1−α/2σci+1,i+1}
- 当 σ 2 \sigma^2 σ2未知时,可用 σ ^ 2 = S e 2 n − p \hat\sigma^2 = \frac{S_e^2}{n-p} σ^2=n−pSe2来估计 σ 2 \sigma^2 σ2,而 S e 2 σ 2 ∼ χ 2 ( n − p ) \frac{S_e^2}{\sigma^2}\sim \chi^2(n-p) σ2Se2∼χ2(n−p),且 S e 2 S_e^2 Se2与 β ^ \hat\beta β^独立,所以 t i = β ^ i σ ^ c i i ∼ t ( n − p ) t_i=\frac{\hat\beta_i}{\hat\sigma\sqrt{c_{ii}}}\sim t(n-p) ti=σ^ciiβ^i∼t(n−p)所以 β i \beta_i βi的 100 ( 1 − α ) % 100(1-\alpha)\% 100(1−α)%的置信区间为 β ^ i ± t 1 − α / 2 ( n − p ) σ ^ c i + 1 , i + 1 , i = 0 , … , p − 1 \hat\beta_i \pm t_{1-\alpha/2}(n-p)\hat\sigma\sqrt{c_{i+1,i+1}}, i=0,\dots,p-1 β^i±t1−α/2(n−p)σ^ci+1,i+1,i=0,…,p−1,所以对于假设 H 0 : β i = 0 v s . H 1 : β i ≠ 0 H_0:\beta_i= 0\ vs.\ H_1:\beta_i\neq 0 H0:βi=0 vs. H1:βi=0,相应的拒绝域为: W = { ∣ β ^ i ∣ > t 1 − α / 2 ( n − p ) σ ^ c i + 1 , i + 1 } W = \{|\hat\beta_i|\gt t_{1-\alpha/2}(n-p)\hat\sigma\sqrt{c_{i+1,i+1}}\} W={∣β^i∣>t1−α/2(n−p)σ^ci+1,i+1}
3.2 模型整体的显著性检验
模型整体的显著性检验就是检验模型的自变量中是否至少有一个变量对因变量有重要的可解释性作用,即检验 H 0 : β 1 = ⋯ = β p − 1 = 0 v s . H 1 : β i ∗ ≠ 0 for some i ∗ ≥ 1 H_0:\beta_1=\dots=\beta_{p-1}=0\ vs.\ H_1: \beta_{i^*}\neq 0\text{ for some }i^*\ge 1 H0:β1=⋯=βp−1=0 vs. H1:βi∗=0 for some i∗≥1
如果原假设成立,那么模型的所有自变量对因变量都没有重要的可解释性作用,亦即所有的自变量都不显著;如果原假设不成立,那么至少有一个模型的自变量对因变量有重要的可解释性作用,要想继续探究是哪个自变量,可继续使用上一小节使用的 t − t- t−检验进行检验。
先说明几个定义:
总平方和(总离差平方和) S T 2 = ∑ i = 1 n ( y i − y ˉ ) 2 S_T^2 = \sum_{i=1}^n(y_i-\bar y)^2 ST2=i=1∑n(yi−yˉ)2
回归平方和 S R 2 = ∑ i = 1 n ( y ^ i − y ˉ ) 2 S_R^2 = \sum_{i=1}^n(\hat y_i-\bar y)^2 SR2=i=1∑n(y^i−yˉ)2
误差平方和(残差平方和) S e 2 = ∑ i = 1 n ( y i − y ^ i ) 2 S_e^2 = \sum_{i=1}^n(y_i-\hat y_i)^2 Se2=i=1∑n(yi−y^i)2
它们之间的关系是 S T 2 = S R 2 + S e 2 S_T^2=S_R^2+S_e^2 ST2=SR2+Se2,可证明如下:
S T 2 = ∑ i = 1 n ( y i − y ˉ ) 2 = ∑ i = 1 n ( y i − y ^ i + y ^ i − y ˉ ) 2 S_T^2=\sum\limits_{i=1}^n (y_i-\bar y)^2 = \sum\limits_{i=1}^n (y_i-\hat y_i+\hat y_i-\bar y)^2 ST2=i=1∑n(yi−yˉ)2=i=1∑n(yi−y^i+y^i−yˉ)2
= ∑ i = 1 n [ ( y i − y ^ i ) 2 + ( y ^ i − y ˉ ) 2 + 2 ( y i − y ^ i ) ( y ^ i − y ˉ ) ] =\sum\limits_{i=1}^n [(y_i-\hat y_i)^2+(\hat y_i-\bar y)^2+2(y_i-\hat y_i)(\hat y_i-\bar y)] =i=1∑n[(yi−y^i)2+(y^i−yˉ)2+2(yi−y^i)(y^i−yˉ)]
= S R 2 + S e 2 + 2 ∑ i = 1 n ( y i − y ^ i ) ( y ^ i − y ˉ ) =S_R^2+S_e^2 +2\sum\limits_{i=1}^n (y_i-\hat y_i)(\hat y_i-\bar y) =SR2+Se2+2i=1∑n(yi−y^i)(y^i−yˉ)
而 ∑ i = 1 n ( y i − y ^ i ) ( y ^ i − y ˉ ) = ϵ ^ ⊤ Y ^ − Y ˉ ∑ i = 1 n ϵ ^ i \sum\limits_{i=1}^n (y_i-\hat y_i)(\hat y_i-\bar y) = \hat\epsilon^\top\hat Y-\bar Y\sum\limits_{i=1}^n \hat\epsilon_i i=1∑n(yi−y^i)(y^i−yˉ)=ϵ^⊤Y^−Yˉi=1∑nϵ^i
= [ ( I n − P ) Y ] ⊤ P Y − Y ˉ ( 1 , 1 , … , 1 ) ( I n − P ) Y = [(I_n-P)Y]^\top PY-\bar Y(1,1,\dots,1) (I_n-P)Y =[(In−P)Y]⊤PY−Yˉ(1,1,…,1)(In−P)Y
= Y ⊤ ( I n − P ) P Y − Y ˉ [ ( 1 , 1 , … , 1 ) − ( 1 , 1 , … , 1 ) P ] Y =Y^\top (I_n-P)PY - \bar Y[(1,1,\dots,1)-(1,1,\dots,1)P]Y =Y⊤(In−P)PY−Yˉ[(1,1,…,1)−(1,1,…,1)P]Y
= 0. = 0. =0.
结论:在模型的假设 ϵ ∼ N ( 0 , σ 2 I n ) \epsilon\sim N(0,\sigma^2I_n) ϵ∼N(0,σ2In)下, S R 2 , S e 2 , y ˉ S_R^2,S_e^2,\bar y SR2,Se2,yˉ相互独立,而且在原假设 H 0 : β 1 = ⋯ = β p − 1 = 0 H_0:\beta_1=\dots=\beta_{p-1}=0 H0:β1=⋯=βp−1=0成立的前提下, S R 2 / σ 2 ∼ χ 2 ( p − 1 ) S_R^2/\sigma^2\sim\chi^2(p-1) SR2/σ2∼χ2(p−1)
下面我们用广义似然比假设检验进行检验。
我们的模型为 Y = X β + ε Y = X\beta + \varepsilon Y=Xβ+ε,根据我们模型的假设 ϵ ∼ N ( 0 , σ 2 I n ) \epsilon\sim N(0,\sigma^2I_n) ϵ∼N(0,σ2In),我们可以知道 Y ∼ N ( X β , σ 2 I n ) Y\sim N(X\beta, \sigma^2I_n) Y∼N(Xβ,σ2In),所以 Y Y Y的似然方程为 L ( β , σ 2 ) = ( 2 π σ 2 ) − n / 2 e − ∣ ∣ Y − X β ∣ ∣ 2 2 σ 2 L(\beta,\sigma^2) = (2\pi \sigma^2)^{-n/2} e^{-\frac{||Y-X\beta||^2}{2\sigma^2}} L(β,σ2)=(2πσ2)−n/2e−2σ2∣∣Y−Xβ∣∣2。
在全假设下,即 Θ = { ( β , σ 2 ) ∣ β ∈ R p , σ 2 > 0 } \Theta=\{(\beta,\sigma^2)|\beta\in \mathbb{R}^p,\sigma^2\gt0\} Θ={(β,σ2)∣β∈Rp,σ2>0},使似然方程 L ( β , σ 2 ) L(\beta,\sigma^2) L(β,σ2)最大化的参数为 β ^ = ( X ⊤ X ) − 1 X ⊤ Y , σ ^ 2 = ∣ ∣ Y − X β ^ ∣ ∣ n = S e 2 n \hat\beta = (X^\top X)^{-1}X^\top Y,\ \hat\sigma^2 = \frac{||Y-X\hat \beta||}{n}=\frac{S_e^2}{n} β^=(X⊤X)−1X⊤Y, σ^2=n∣∣Y−Xβ^∣∣=nSe2。
在原假设下,即 Θ 0 = { ( β , σ 2 ) ∣ β i = 0 , i ≥ 1 , σ 2 > 0 } \Theta_0=\{(\beta,\sigma^2)|\beta_i=0,i\ge 1,\sigma^2\gt0\} Θ0={(β,σ2)∣βi=0,i≥1,σ2>0},使似然方程 L ( β , σ 2 ) L(\beta,\sigma^2) L(β,σ2)最大化的参数为 β ^ ∗ = ( Y ˉ , 0 , … , 0 ) ⊤ , σ ^ ∗ 2 = ∣ ∣ Y − X β ^ ∗ ∣ ∣ n = S T 2 n \hat\beta^* = (\bar Y,0,\dots,0)^\top,\ \hat\sigma^{*2} = \frac{||Y-X\hat \beta^*||}{n}=\frac{S_T^2}{n} β^∗=(Yˉ,0,…,0)⊤, σ^∗2=n∣∣Y−Xβ^∗∣∣=nST2
所以似然比统计量为: λ = sup θ ∈ Θ L ( β , σ 2 ) sup θ ∈ Θ 0 L ( β , σ 2 ) = L ( β ^ , σ ^ 2 ) L ( β ^ ∗ , σ ^ ∗ 2 ) = ( S T 2 S e 2 ) n / 2 = ( 1 + S R 2 S e 2 ) n / 2 (7) \lambda = \frac{\sup_{\theta\in\Theta}L(\beta,\sigma^2)}{\sup_{\theta\in\Theta_0}L(\beta,\sigma^2)} = \frac{L(\hat\beta,\hat\sigma^2)}{L(\hat\beta^*,\hat\sigma^{*2})}= \left(\frac{S_T^2}{S_e^2}\right)^{n/2}= \left(1+\frac{S_R^2}{S_e^2}\right)^{n/2}\tag7 λ=supθ∈Θ0L(β,σ2)supθ∈ΘL(β,σ2)=L(β^∗,σ^∗2)L(β^,σ^2)=(Se2ST2)n/2=(1+Se2SR2)n/2(7)
根据我们以上的结论,原假设成立时, S R 2 / σ 2 ∼ χ 2 ( p − 1 ) S_R^2/\sigma^2\sim\chi^2(p-1) SR2/σ2∼χ2(p−1), ( n − p ) σ ^ 2 σ 2 = S e 2 σ 2 ∼ χ 2 ( n − p ) \frac{(n-p)\hat\sigma^2}{\sigma^2}=\frac{S_e^2}{\sigma^2}\sim \chi^2(n-p) σ2(n−p)σ^2=σ2Se2∼χ2(n−p),所以 F = S R 2 / ( p − 1 ) S e 2 / ( n − p ) ∼ F ( p − 1 , n − p ) F=\frac{S_R^2/(p-1)}{S_e^2/(n-p)}\sim F(p-1,n-p) F=Se2/(n−p)SR2/(p−1)∼F(p−1,n−p)。所以我们把 F F F作为检验统计量,拒绝域为 W = { F > C } W=\{F\gt C\} W={F>C}, C = F 1 − α ( p − 1 , n − p ) C = F_{1-\alpha}(p-1,n-p) C=F1−α(p−1,n−p), sup θ ∈ Θ 0 P θ ( F > C ) = α \sup\limits_{\theta\in\Theta_0}P_\theta(F\gt C)=\alpha θ∈Θ0supPθ(F>C)=α。
进而可以得到方差分析表。
3.3 模型优度
定义判决系数为 R 2 = S R 2 S T 2 = 1 − S e 2 S T 2 R^2 =\frac{S_R^2}{S_T^2}=1-\frac{S_e^2}{S_T^2} R2=ST2SR2=1−ST2Se2,一般认为 R 2 R^2 R2越大,模型越好。但注意到,回归变量越多,残差平方和 S e 2 S_e^2 Se2越小, R 2 R^2 R2自然变大。所以用 R 2 R^2 R2衡量模型拟合程度显然不合适,没有考虑过度拟合的情况。另一方面,随机项的方差的估计 σ ^ 2 = S e 2 / ( n − p ) \hat\sigma^2=S_e^2/(n-p) σ^2=Se2/(n−p)并不一定随着 S e 2 S_e^2 Se2越小而变小,因为变量个数 p p p变大,该估计量的分母变小。
于是,有人提出对 R 2 R^2 R2进行调整,调整的决定系数(adjusted R-squared)为 R ~ 2 = 1 − S e 2 / ( n − p ) S T 2 / ( n − 1 ) = 1 − ( n − 1 ) S e 2 ( n − p ) S T 2 \tilde{R}^2 = 1-\frac{S_e^2/(n-p)}{S_T^2/(n-1)}=1-\frac{(n-1)S_e^2}{(n-p)S_T^2} R~2=1−ST2/(n−1)Se2/(n−p)=1−(n−p)ST2(n−1)Se2。不难看出 R ~ 2 < R 2 \tilde{R}^2\lt R^2 R~2<R2,这是因为 p > 1 p\gt 1 p>1, 调整的决定系数 R ~ 2 \tilde{R}^2 R~2综合考虑了模型变量个数 p p p的影响。
3.4 预测
先说明一下置信区间与预测区间的区别:
置信区间:利用估计的回归方程,对于自变量 x x x的一个给定值 x 0 x_0 x0,求出因变量 y y y的平均值的估计区间。
预测区间:利用估计的回归方程,对于自变量 x x x的一个给定值 x 0 x_0 x0,求出因变量 y y y的一个个别值的估计区间。
而且,预测区间的范围总是要比置信区间的范围要大的。就是说,给定一个 x x x,估计对应 y y y的平均值平均值比估计一个个别值更精确一点。其实也好理解,估计平均值比估计个别值貌似更简单一点嘛。个别值更容易受一些外界因素影响而有差异性,而平均值则相对稳定些。
E [ y n + 1 ] E[y_{n+1}] E[yn+1]的置信区间:
考虑到 y n + 1 = β 0 + β 1 x n + 1 , 1 + ⋯ + β p − 1 x n + 1 , p − 1 + ϵ n + 1 y_{n+1} = \beta_0+\beta_1x_{n+1,1}+\dots+\beta_{p-1}x_{n+1,p-1}+\epsilon_{n+1} yn+1=β0+β1xn+1,1+⋯+βp−1xn+1,p−1+ϵn+1,而在模型 ϵ ∼ N ( 0 , σ 2 I n ) \epsilon\sim N(0,\sigma^2I_n) ϵ∼N(0,σ2In)的假设下,可以知道 y n + 1 = v ⊤ β + ϵ n + 1 ∼ N ( v ⊤ β , σ 2 ) y_{n+1}=v^\top \beta+\epsilon_{n+1}\sim N(v^\top \beta,\sigma^2) yn+1=v⊤β+ϵn+1∼N(v⊤β,σ2),其中, v = ( 1 , x n + 1 , 1 , x n + 1 , 2 , … , x n + 1 , p − 1 ) ⊤ v = (1,x_{n+1,1},x_{n+1,2},\dots,x_{n+1,p-1})^\top v=(1,xn+1,1,xn+1,2,…,xn+1,p−1)⊤。
所以 E [ y n + 1 ] = v ⊤ β E[y_{n+1}]=v^\top \beta E[yn+1]=v⊤β的无偏估计为 y ^ n + 1 = v ⊤ β ^ ∼ N ( v ⊤ β , σ 2 v ⊤ ( X ⊤ X ) − 1 v ) \hat y_{n+1} = v^\top \hat\beta \sim N(v^\top \beta, \sigma^2 v^\top(X^\top X)^{-1}v) y^n+1=v⊤β^∼N(v⊤β,σ2v⊤(X⊤X)−1v),而根据 2.5 2.5 2.5节的几个结论,可以知道, y ^ n + 1 − v ⊤ β σ ^ v ⊤ ( X ⊤ X ) − 1 v ∼ t ( n − p ) \frac{\hat y_{n+1}-v^\top \beta}{\hat\sigma\sqrt{v^\top(X^\top X)^{-1}v}}\sim t(n-p) σ^v⊤(X⊤X)−1vy^n+1−v⊤β∼t(n−p)。所以 E [ y n + 1 ] E[y_{n+1}] E[yn+1]的 100 ( 1 − α ) % 100(1-\alpha)\% 100(1−α)%置信区间为: y ^ n + 1 ± t 1 − α / 2 ( n − p ) σ ^ v ⊤ ( X ⊤ X ) − 1 v (8) \hat y_{n+1}\pm t_{1-\alpha/2}(n-p)\hat{\sigma}\sqrt{v^\top(X^\top X)^{-1}v}\tag8 y^n+1±t1−α/2(n−p)σ^v⊤(X⊤X)−1v(8)
y n + 1 y_{n+1} yn+1的预测区间:
同样的,可以得到 y n + 1 − y ^ n + 1 σ ^ 1 + v ⊤ ( X ⊤ X ) − 1 v ∼ t ( n − p ) \frac{y_{n+1}-\hat y_{n+1}}{\hat{\sigma}\sqrt{1+v^\top(X^\top X)^{-1}v}}\sim t(n-p) σ^1+v⊤(X⊤X)−1vyn+1−y^n+1∼t(n−p),所以 y n + 1 y_{n+1} yn+1的 100 ( 1 − α ) % 100(1-\alpha)\% 100(1−α)%预测区间为: y ^ n + 1 ± t 1 − α / 2 ( n − p ) σ ^ 1 + v ⊤ ( X ⊤ X ) − 1 v (9) \hat y_{n+1}\pm t_{1-\alpha/2}(n-p)\hat{\sigma}\sqrt{1+v^\top(X^\top X)^{-1}v}\tag9 y^n+1±t1−α/2(n−p)σ^1+v⊤(X⊤X)−1v(9)
参考:
更多推荐


所有评论(0)