1、1第二章 多元线性回归模型的有偏估计模型的参数估计依赖于观测样本,样本是随机的(至少 Y 是随机的),因此估计量也是随机的,不一定恰好等于被估计参数的真值。但是我们希望多次估计的结果的期望值接近或等于真值,即 2)( ,)(E这就叫无偏估计。无偏估计被认为是一个估计量应有的优良性质。但是在一些场合,满足无偏性的估计量却不具备其它应有的优良性,比如说稳定性、容许性。统计学家提出了一些新的估计方法,它们往往不具备无偏性,但在特定场合综合起来考虑还是解决问题较好的。本章就分别介绍这些特定场合下的有偏估计。第一节 设计矩阵列复共线与岭回归一、设计矩阵列复共线的影响上一章最后一节讲的是设计矩阵列向量完全
2、线性相关,|X X|=0 的情况。实际工作中常遇到的是,设计矩阵的列向量存在近似线性相关(称为复共线 (multicollinearity),|X X|0。此时一般最小二乘方法尽管可以进行,但估计的性质变坏,主要是对观测误差的稳定性变差,严重时估计量可能变得面目全非。例如我们建立二元线性回归模型(2.1.1)21XY有关资料在下面运算过程可以看到。看一看原始资料,它近似满足 Yi=X1i+X2i, 应该估计出。可是我们调用普通最小二乘回归程序,运算结果却是1 , ,021(2.1.2)2156.430.0Y对现有数据拟合的还挺好,两条曲线几乎成了一条曲线 (图 2.1.1.1),F 值为 30
3、3744,但是代入 X1=0, X2=10,预测值却为 15.66,这与原模型应有的预测值 10 相距甚远。-岭回归与岭迹图计算程序, 例 2.1.4 例 214.D 数据文件中, n=8, M=2 要显示原始资料吗? 0=不显示, 1=显示 2.0100 .9900 1.010021.9900 1.0200 .99004.0100 2.0300 1.99005.9900 2.9700 3.01008.0100 3.9600 4.01007.9900 4.0100 3.990010.0100 5.0400 4.990011.9900 6.0500 5.9900正规方程系数矩阵的行列式的值是 2
4、.12162请输入工作参数, 0=普通回归, 1=岭回归, 2= 计算岭迹 (0)现在作线性回归显著性检验, 计算 t,F,R 统计量请输入显著性水平 a, 通常取 a=0.01, 0.05, 0.10, a=? (0.05)-线 性 回 归 分 析 计 算 结 果 样本总数 8 自变量个数 2-回归方程 Y = b0+b1*X1+.+b2*X2Y = .0033 + .4330 X1 + 1.5660 X2 回归系数 b0, b1, b2, ., b2.0033 .4330 1.5660-残差平方和: .00 回归平方和: 93.92误差方差的估计 : .0001 标准差 = .0098-线
5、 性 回 归 显 着 性 检 验 显著性水平 : .050-回归方程整体显著性 F 检验, H0:b0=b1=.=b2=0F 统计量: 303744.5000 F 临界值 F(2, 5) 5.786全相关系数 R : 1.0000-回归系数逐一显著性 t 检验, H0:bi=0, i=1,.,2t 临界值 t( 5) 2.015回归系数 b1-b 2 的 t 值: .0106 .0382-要作回归预测吗? 键入 0=不预测, 1=要预测 (1)现在作回归预测, 请输入自变量, X1-X 2X( 1)= 0 X( 2)= 10-线 性 回 归 预 测: Y 的预测值= 15.6633给定 X1-
6、X 2= .0000 10.0000-3要作回归预测吗? 键入 0=不预测, 1=要预测 (0)要打印拟合数据吗? 0=不打印, 1=打印 (1)Y 的观测值 Y 的拟合值 差值2.0100 2.0136 -.00361.9900 1.9953 -.00534.0100 3.9987 .01135.9900 6.0030 -.01308.0100 7.9977 .01237.9900 7.9881 .001910.0100 10.0001 .009911.9900 12.0035 -.0135计算结束。 -下面显示拟合图像。图 2.1.1.1024681012141 2 3 4 5 6 7 8
7、原 始 数 据拟 合 数 据对此我们可以作如下理论分析。 作为 的估计是否优良,应该考察它与 的接近程度,这可以用 的均方误差 (Mean Square Error)来度量:(2.1.3))()|(|)2EMSE我们来计算线性模型(2.1.4)nIVarXY2)( ,0)( , 的 MSE( )。由于(2.1.5)XXYX 111 )()()()(4故由公式 得)()()( yAVartEyAyE(2.1.6)12222 )()()( XtrXtrXtrMS 进一步,若 N (0, 2In),则因对于对称矩阵 A 有(2.1.7)pjivuvitjijaA1,)( 注意到(2.1.8)ujii
8、jiEitj 或, 0, 3)(4故(2.1.9)2241,211, 241242)( )(3)(trAtaaApjiipiipji jiijjipii 因此 2422)(trAEAVar (2.1.10)于是24224 )()()(2| XtrXXtrar(2.1.11)由于 XX 为正定阵,其特征根皆为正数,设为 1 2 p0,则pitr1)((2.1.12)代入(2.1.6)与(2.1.11)得)|(|)2EMS(2.1.13)5piVar1242)|(| (2.1.14)当设计矩阵 X 的列向量存在复共线关系时, p0, 很大,就使 与p )|(|2E都很大。尽管这时按平均来说, 的无
9、偏估计,但具体在每一次计算,)|(|2Var 是由于均方误差太大,使得 估计值偏差很大,以致前面的数值例子变得面目全非。二、岭回归统计学界由 A. E. Hoerl 在 1962 年提出并和 R. W. Kennard 在 1970 年系统发展的岭回归(Ridge Regression)方法,可以显著改善设计矩阵列复共线时最小二乘估计的均方误差,增强估计的稳定性。这个方法在计算数学称为阻尼最小二乘,出现得较早一些。岭回归方法主要就是在病态的(XX) 中沿主对角线人为地加进正数,从而使 p 稍大一些。我们知道模型(2.1.4)中 的最小二乘估计为 YX1)((2.1.15)则 的岭估计定义为 k
10、kIXkp0 ,)()1(2.1.16)从式子直觉看出,当 k=0 时,它就是最小二乘估计;当 k+, 。于是就要0)(k问 k 究竟取多大值为好? 同时我们也要知道 的统计性质究竟如何。)(k性质 1. 岭估计不再是无偏估计,即 。E因为 1111)()()( ()( XkIXkIXYkEpp无偏性一直被认为是一个好的统计量所必须具有的基本性质,但是在现在所讨论的问题场合,我们只好牺牲无偏性,以改善估计的稳定性。性质 2. 岭估计是线性函数。记 S=XX, Zk= (I+kS-1)-1, 则因 LkZkSIYXSkIY )()( 111 6(2.1.17)可见 不仅是 Y 的线性函数,而且是
11、原来最小二乘估计 的线性函数。)(k L性质 3. Zk 的特征根都在 (0,1)内。设有正交阵 P 与 P使 ),( 11 ppdiagS (2.1.18)则 )(1 )()(1111 kkkkIPSIPSIPZ ppK 故知 Zk 的特征根分别为 ,都在(0,1)内。ki性质 4. 岭估计是压缩估计,即 。|)(|k这是因为由性质 2、性质 3,有 2222 |)(|)(|)(| PkPZkk当然,由于 并不一定互相相等,这种压缩一般不是各方向上的均匀压缩。pii ,1,性质 5. 岭估计的均方误差较小,即 22|)(| EkE(2.1.19)这个性质我们放到下一章与方差分量线性模型的广义
12、岭估计一起证明。三、岭迹分析与岭参数选择因为岭估计 YXkIk1)()(2.1.20)7是 k 的函数,所以在二维坐标平面上若以横轴为 k,纵轴为 ,它将画出一条曲线。这条)(k曲线我们称之为岭迹。前已指出,当 k0 时,岭迹反映了最小二乘估计 的不稳定性。当 k+时,岭迹将2趋于 0。在 k 从 0 到+的变化过程中, 的变化可能比较复杂。)(k关于岭参数选择的问题,已有许多文献讨论。讨论来讨论去,并没有一个公认最优的准则。许多办法含有未知参数,又要对其进行估计。下一章我们将仔细介绍在方差分量模型的广义岭估计中岭参数的选择办法。这里我们只简要介绍几种较有影响的方法和原则。1. 岭迹稳定观察岭
13、迹曲线,原则上应该选取使 稳定的最小 k 值,同时残差平方和也不增加太)(k多。2. 均方误差小岭估计的均方误差 还是 k 的函数,可以证明它能在某处2|)(|)(kEkMS取得最小值。计算并观察 ,开始它将下降,到达最小值后开始上升。取它最小处的 k 作为岭参数。3. 22max/i假设回归模型 Y=X + 的设计阵 X 已中心化,即 0121 pX并设 P 为正交方阵,使(2.1.21)),()(21pdiagP记 =P , 称为典则参数,Z =XP,则原模型变为 ZY(2.1.22)这个形式被称为线性回归的典则形式。此时 的最小二乘估计与岭回归估计为 YXP11)((2.1.23)(2.
14、1.24)kIkIZk )()()1于是 ,(1 12 ZpnYpnn ),(11pYXG8都是可以计算的,从而选取岭参数 。22max/ik4 pjjk122/这是 Bayes 原理推出的法则。假若 为 i的公共先验方差,则从 Bayes 原理出发,2得到的岭估计为 2*/k(2.1.25)而 正是 k*的一个估计量。pjjk122/5. pjj12/直观考虑是,当 XX=I 时,取 可使岭估计具有最小的均方误差。于是pjjk12/以 即得这个估计量。2替 代因为编制岭回归计算程序及自动打印岭迹的程序现在并不困难,所以在现在条件下我们建议主要采取岭迹图标分析。四、广义岭回归前面我们介绍了线性
15、回归模型的典则形式 ZY(2.1.26)其中 =P 称为典则参数,Z= XP称为典则变量,P 为正交方阵使 P(XX) P= 。此时 的岭估计为 YkI1)()(2.1.27)这里主对角线上是统一地加上相同的 k。 如果灵活一些,使主对角线上可以加上不同的 ki, i=1,p, 显然有可能使均方误差进一步下降。而且,原来狭义岭估计是广义岭估计的特例。将这个思想写成式子就是 YXPK1)()9(2.1.28)回到原来参数就是 YXPKPK1)()((2.1.29)这里 ),(21pkdiag(2.1.30)广义岭回归确实能使估计的均方误差进一步下降,岭参数的选择更为复杂一些。这方面的证明与讨论我
16、们放到下一章结合方差分量模型一起进行。算例 2.1.4 岭回归与岭迹图仍然采用本节开头的资料,利用本书所附软件运行如下。从岭迹资料看,当岭参数 k 从0 变到 0.04 时,就已趋于稳定了,到 k=0.2 时,已相当稳定了。我们就取 k=0.20 作正式的拟合回归分析。岭迹图也打印在后面(图 2.1.4.1)。看回归方程系数,一个 0.9816,一个1.0060,相当接近真值。实际回归计算工作中如果遇到列复共线情况,就改用本节程序运行。-岭回归与岭迹图计算程序, 例 2.1.4 正规方程系数矩阵的行列式的值是 2.12162请输入工作参数, 0=普通回归, 1=岭回归, 2= 计算岭迹 (2)
17、下面计算岭迹, 计算 40 个点, 需要决定岭参数起点 K1 与终点 K40例数据文件的岭参数, K1=0, K40=0.6 (下面只打印前 20 个点)序号 岭参数 Beta0 Beta1 Beta2 0 .00 .003 .433 1.5661 .02 .010 .846 1.1492 .03 .013 .911 1.0843 .05 .015 .937 1.0574 .06 .017 .951 1.0425 .08 .019 .960 1.0326 .09 .021 .966 1.0267 .11 .023 .970 1.0218 .12 .025 .973 1.0179 .14 .02
18、6 .976 1.01410 .15 .028 .977 1.01211 .17 .030 .979 1.01012 .18 .031 .980 1.00813 .20 .033 .981 1.00614 .21 .035 .982 1.0051015 .23 .036 .983 1.00416 .24 .038 .983 1.00317 .25 .039 .984 1.00218 .27 .041 .984 1.00119 .29 .042 .985 1.00020 .30 .044 .985 .999计算指定岭参数的岭回归, 请输入岭回归合适的岭参数 k, k=? (0.2)现在作线性回归
19、显著性检验, 计算 t,F,R 统计量请输入显著性水平 a, 通常取 a=0.01, 0.05, 0.10, a=? (0.05)-线 性 回 归 分 析 计 算 结 果 样本总数 8 自变量个数 2-回归方程 Y = b0+b1*X1+.+b2*X2Y = .0336 + .9816 X1 + 1.0060 X2 回归系数 b0, b1, b2, ., b2 .0336 .9816 1.0060-残差平方和: .01 回归平方和: 93.21误差方差的估计 : .0007 标准差 = .0264-线 性 回 归 显 着 性 检 验 显著性水平 : .050-回归方程整体显著性 F 检验, H
20、0:b0=b1=.=b2=0F 统计量: 41914.2000 F 临界值 F(2, 5) 5.786全相关系数 R : 1.0000-回归系数逐一显著性 t 检验, H0:bi=0, i=1,.,2t 临界值 t( 5) 2.015回归系数 b1-b 2 的 t 值: .0242 .0246-要作回归预测吗? 键入 0=不预测, 1=要预测 现在作回归预测, 请输入自变量, X1-X 2X( 1)= 0 X( 2)= 10-线 性 回 归 预 测: Y 的预测值= 10.0932给定 X1-X 2= .0000 10.0000-要作回归预测吗? 键入 0=不预测, 1=要预测 (0)要打印拟合数据吗? 0=不打印, 1=打印 (0)计算结束。