第四章 方差分量线性回归模型

第四章 方差分量线性回归模型 本章考虑的线性模型不仅有固定效应、随机误差，而且有随机效应。我们先从随机效应角度理解回归概念，导出方差分量模型，然后研究模型三种主要解法。最后本章介绍关于方差分量模型的两个前沿研究成果，是作者近期在《应用数学学报》与国际数学杂志《Communications in Statistics》上发表的。 第一节 随机效应与方差分量模型 一、随机效应回归模型 前面所介绍的回归模型不仅都是线性的，而且自变量看作是固定效应。我们从资料对 出发建立回归模型，过去一直是把Y看作随机的，X1，…，Xp看作非随机的。但是实际上，自变量也经常是随机的，而并不是我们可以事先设计好的设计矩阵。我们把自变量也是随机变量的回归模型称为随机效应回归模型。 究竟一个回归模型的自变量是随机的还是非随机的，要视具体情况而定。比如一般情况下消费函数可写为 （4.1.1） 这里X是居民收入，T是税收，C0是生存基本消费，b是待估系数。加上随机扰动项，就是一元线性回归模型 （4.1.2） 那么自变量到底是固定效应还是随机效应?那要看你采样情况。如果你是按一定收入的家庭去调查他的消费，那是取设计矩阵，固定效应。如果你是随机抽取一些家庭，不管他收入如何都登记他的收入与消费，那就是随机效应。 对于随机效应的回归模型，我们可以从条件期望的角度推导出与最小二乘法则等价的回归函数。 我们希望通过X预测Y，也就是要寻找一个函数 ，当X的观察值为x时，这个预测的误差平均起来应达到最小，即 （4.1.3） 这里min是对一切X的可测函数L(X)取极小。由于当 （4.1.4） 时，容易证明 （4.1.5） 故当 时， （4.1.6） 要使上式左边极小，只有取 。 这个结果告诉我们，预测函数取作条件期望E(Y|X)时，可使预测误差最小。我们还可以证明，此时M(X)=E(Y|X)与Y具有最大相关，即 （4.1.7） 这里ρ表示相关系数。 这是因为当 时，易证 ,同时 ,于是 等号当且仅当 （4.1.8） 时成立，此时L(X)是M(X)的线性函数。 (4.1.3)与(4.1.7)表达了 的极好性质，我们称 （4.1.9） 为Y关于X的回归曲线。 上面的L(X)可取一切函数。如果限定L(X)是X的线性函数，即要限定 （4.1.10） 这里 是对X的一切线性函数取极小，则称满足上式的线性函数为Y关于X的回归直线。我们可以求出 的解。记 ，则 （4.1.11） 这里 （4.1.12） （4.1.13） （4.1.14） 对L(β0，β)求微分(矩阵微商 公式 小学单位换算公式大全免费下载公式下载行测公式大全下载excel公式下载逻辑回归公式下载 )得： （4.1.15） 解得 （4.1.16） 这里当然假定 存在，否则使用广义逆。 此时的预测误差方差是 （4.1.17） （4.1.18） 为复相关系数。它指出了Y与多元变量 之间的线性相关程度，是一元相关系数 （4.1.19） 的推广。 从条件期望角度我们导出的随机效应回归模型的回归直线表达式，与从最小二乘角度导出的固定效应的回归方程，表达式是等价的，所以从计算角度，我们不怎么区分。 二、方差分量模型概念 上段我们建立了随机效应概念，将自变量也视作随机变量，这就可以导出方差分量模型。方差分量模型研究工作的奠基人是我国最早的统计学家许宝驭马录先生。 还是刚才提到的消费函数回归模型，我们作随机抽样。考虑居民按职业的分类，如工人、教师、医生、律师、店员等等，记为 ，我们从这些职业中随机抽取了n个样本，则模型可写为 （4.1.20） 这里Xi可看作是第i种职业对收入的效应。如果我们事先安排好取哪个职业的，当然Xi是固定效应。可是我们现在对职业选取是随机的，而且我们还想研究职业效应的方差，这就导入了方差分量模型，因为现在Cij的方差由两部分组成： （4.1.21） 为了数学符号统一，我们将经济学中的符号改过来，刚才建立的模型是 （4.1.22） 它有一项固定效应μ，一项随机效应ξ1，一项随机误差ε。如果还要考虑地区因素对消费的影响，还可以加进第二个随机效应ξ2，于是可得模型 （4.1.23） 这次我们省掉了取值的标记，Y的方差由三项组成。 一般地，我们建立方差分量模型如下： （4.1.24） 这里有固定效应向量β，随机效应向量 （4.1.25） 并且将随机误差项ε也并入了随机效应向量去。设计矩阵X以及 （4.1.26） 都是已知的。对于随机效应 ,合理的假定是 （4.1.27） 当然以后有时还可以考虑ξi是向量的情况，不过这里假定每个ξi是一维变量。记 ， （4.1.28） 则方差分量模型可记为 （4.1.29） 模型的主要任务是要估计固定效应向量β与方差分量 。和一般的多元线性回归模型相比，就是待估的方差多了。 通过这些介绍，我们就可以方便地将各种经济方面的普通线性回归模型改造成方差分量模型，当然要根据实际。 第二节 方差分量模型的解法 对于方差分量模型 （4.2.1） 一般都采用二步估计法，首先估计方差分量 ，然后再估计固定效应β。按照广义最小二乘 （4.2.2） 其中 （4.2.3） 所以方差分量模型解法的关键是估计方差分量。以下介绍的方法，也都是针对方差分量估计方法而言的。 一、方差 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 法 先从一个简单的模型结合数据结构形象地说明方法。考虑模型 （4.2.4） β0为总平均，是固定效应，ξ1,…，ξm是随机效应， ， 。对于随机误差 。这个模型如果记作方差分量模型的标准形式是 （4.2.5） 其中设计阵X=(1,1,…，1)′，随机效应矩阵为 EMBED Equation.3 （4.2.6） 我们手中资料只有 我们采用(4.2.4)记法方便一些，将资料Y排成表 j i 1 2 … k 组内平均 1 Y11 Y12 … Y1k 2 Y21 Y22 … Y2k m Ym1 Ym2 … Ymk 方差分析主要掌握三点，一是计算组内差、组间差，二是作平方和分解，三是计算各自的自由度。 先计算总平均： （4.2.7） 总变差(全体资料与总平均的偏差平方和)： （4.2.8） 各组平均(各组资料横向相加并平均) （4.2.9） 组间差(各组平均数与总平均数的偏差平方和) （4.2.10） 组内差(各组数据与本组平均数的偏差平方和) （4.2.11） 则必有平方和分解 （4.2.12） 将各平方和除以各自的自由度。ST有一个约束 (4.2.7)，自由度为 ；SA有m组差，1个约束，自由度为m －1；Se有mk组差，m个约束，自由度为mk-m。注意有自由度分解： （4.2.13） 于是算出均方： （4.2.14） （4.2.15） （4.2.16） 因为假定为随机效应，可以算出各均方的均值： （4.2.17） （4.2.18） 以 代者 ， 代替 ，得方程组： （4.2.19） 解得 （4.2.20） 这样就作好了方差分量的估计，然后可以按(4.2.2)作出β的估计。因为这里的方差分量是由方差分析法作出的，故称为方差分析法。 推广到一般的方差分量模型时，基本原则是类似的。我们不妨考虑方差分量模型 （4.2.21） 先对总平方和Y′Y作平方和分解 （4.2.22） 其中Sβ是在模型Y=Xβ+ε中，β的回归平方和： （4.2.23） 是在模型 中，消去β影响后ξ1的平方和 （4.2.24） 类似地， 是在模型 中消去β和ξ1影响后，ξ2的平方和： （4.2.25） 最后的Sε为残差平方和 （4.2.26） 可以验证 （4.2.27） （4.2.28） （4.2.29） （4.2.30） 这里 （4.2.31） （4.2.32） （4.2.33） 这里P*表示关于*的投影阵。 下面计算各平方和的均值。 （4.2.34） 因为 ，所以上式第一项为0。在第三项中， （4.2.35） 在第六项中 （4.3.36） 所以最后有 （4.3.37） 其中 （4.2.38） （4.2.39） （4.2.40） （4.2.41） 类似还可以求得 （4.2.42） （4.2.43） （4.2.44） 于是我们得到方程组 （4.2.45） 解此方程组，就可以得到 的估计。然后进入二步估计的第二步，就可以得到关于固定效应的估计。 算例4.2.1 市场收益率与股利和换手率的关系 考虑一个随机效应的多元线性模型 U的形式如同(4.2.6)。 问题的实际背景是，观测对象被分成了m组，可能存在一个随机效应向量对各组资料有不同的作用。模型也可以写作 数据结构及具体数值如下表所示，m=6,k=6。这些资料采自《'96上海股票市场资料总汇》。 我们研究目的一是看过去一年的股利收入与当年换手率对当年市场收益率有何影响，二是想知道是否存在一个潜在的尚未观测到的随机效应，对行业有明显影响。当然这种情况采用方差分量模型比较合适。 要注意本例是两个方差量，上一章第二节模型(3.2.10)也是两个待估的方差量。它们的随机效应作用范围不一样，不是一回事。 表4.2.1 1996年股市资料 类别 股号 股名 1996年收益率% 1995年股利% 1996年日换手率 商 业 类 628 新世界 64.769 20 3.12 631 中百一店 46.845 11.8 1.68 632 华联商厦 41.958 11.3 1.81 655 豫园商城 16.195 11.2 1.10 682 新百公司 79.911 5.2 3.36 694 大连商场 91.388 5.8 4.26 电 子 类 602 真空电子 33.112 10 3.52 651 飞乐音响 8.108 0 1.95 800 天津磁卡 271.763 5 3.74 839 四川长虹 381.686 60 4.41 850 华东计算机 14.638 1 3.27 870 厦华电子 68.579 3.9 4.20 化 工 类 617 联合化纤 -21.871 0.5 1.53 618 卤碱化工 22.370 2 2.63 672 广东化纤 11.860 0 4.65 688 上海石化 179.817 3.5 4.38 886 湖北兴化 236.328 52.3 4.45 889 南京化纤 -33.122 2.4 4.64 4.11医 药 类 664 哈医药 191.666 5 4.32 671 天目药业 111.135 16 4.11 812 华北制药 152.015 8.4 4.11 842 中西药业 13.821 2.6 1.71 849 四药股份 17.892 2.5 1.68 779 四川制药 -24.744 0 12.28 钢 铁 类 608 异型钢管 8.389 10 2.24 665 沪昌特钢 75.39 1.87 4.01 674 四川峨铁 35.932 3 3.64 808 马钢股份 86.528 0.5 4.52 845 钢管股份 -25.170 0 1.61 894 广钢股份 51.371 2.7 7.08 机 械 类 604 二纺机 14.4 10 2.31 605 轻工机械 6.122 0 3.50 610 中纺机 0.701 0 2.27 806 昆明机床 41.852 1.1 4.22 841 上柴股份 66.981 20 2.20 862 南通股份 41.093 20 1.36 首先我们作普通最小二乘回归，得到 ，然后计算 。此时的 已消除固定效应影响，我们将它排成平面表，以作方差分析，计算 与 。计算过程从(4.2.7)至(4.2.20)。 从下面计算过程可以看到，平方和分解式是满足的： 即(ST=SA+Sε)。对于本例资料，随机误差4055.6远大于随机效应方差，组内差远大于组间差，可以认为随机效应不明显，即行业差别不明显。对于选定的方差分量模型，回归结果是 它的标准差很小，为1.0084，这正是采用方差分量模型广义最小二乘意义所在。拟合效果图(图4.2.1.1)令人满意。 ------------------------------------------------------------------------------------------------------------------------- 方差分量模型方差分析法计算程序, 例 4.2.1. 第一列为 Y, 以后各列为 X 例421.D 数据文件中, m=6, k=6, p=2 要显示原始资料吗? 0=不显示, 1=显示 （0） 先作普通最小二乘, 并打印结果: 现在作线性回归显著性检验, 计算t,F,R 统计量 请输入显著性水平a, 通常取a=0.01, 0.05, 0.10, a=? （0.05） ----------------------------------------------------- 线 性 回 归 分 析 计 算 结 果 样本总数 36 自变量个数 2 ----------------------------------------------------- 回归方程 Y = b0+b1*X1+...+b2*X2 Y = 5.3188 + 4.5656 X1 + 6.0128 X2 回归系数 b0, b1, b2, ..., b2 5.3188 4.5656 6.0128 ----------------------------------------------------- 残差平方和: 149400.60 回归平方和: 131759.30 误差方差的估计 : 4150.0170 标准差 = 64.4206 ----------------------------------------------------- 线 性 回 归 显 着 性 检 验 显著性水平 : .050 ----------------------------------------------------- 回归方程整体显著性F检验, H0:b0=b1=...=b2=0 F统计量: 14.5517 F临界值F(2, 33) 3.285 全相关系数 R : .6846 ----------------------------------------------------- 回归系数逐一显著性t检验, H0:bi=0, i=1,...,2 t 临界值 t( 33) 1.6924 回归系数b1-b 2的t值: 5.6318 1.1073 ----------------------------------------------------- 打印方差分析资料 Y(I,J) -50.6210 -39.0275 -38.6721 137.5441 -56.0541 -50.4640 -22.4490 -8.9357 -7.8936 8.0546 37.4223 -20.2416 -25.8348 221.1285 -21.4183 83.6328 -4.9701 -18.2668 -46.8722 75.9169 132.1827 -13.6501 51.7486 6.1371 30.6483 -14.9082 -34.5268 -8.9422 -40.1694 -42.8772 33.9744 20.2007 -77.2975 -103.9000 -8.8455 -81.7529 计算各种平均: 总平均YYba: .0000 各组平均Yba: -16.216 -2.341 39.045 34.244 -18.463 -36.270 计算各种变差: 总变差SST, 组间差SSA, 组内差SSE SST= 149400.6000 SSA= 27731.9200 SSE= 121668.7000 打印方差分量的估计 SIGMAA= 248.4600 SIGMAE= 4055.6240 下面计算协方差阵及其逆的一块,并作分解 计算广义最小二乘估计, 模型转换Y=PY, X=PX, 下面打印矩阵 P 的一块, P’P=∑的逆 SIG1= .1239 -.0007 -.0007 -.0007 -.0007 -.0007 SIG1= .0000 .1238 -.0008 -.0008 -.0008 -.0008 SIG1= .0000 .0000 .1237 -.0008 -.0008 -.0008 SIG1= .0000 .0000 .0000 .1237 -.0008 -.0008 SIG1= .0000 .0000 .0000 .0000 .1236 -.0009 SIG1= .0000 .0000 .0000 .0000 .0000 .1235 下面打印的是广义最小二乘的统计结果: 现在作线性回归显著性检验, 计算t,F,R 统计量 请输入显著性水平a, 通常取a=0.01, 0.05, 0.10, a=? (0.05) ----------------------------------------------------- 线 性 回 归 分 析 计 算 结 果 样本总数 36 自变量个数 2 ----------------------------------------------------- 回归方程 Y = b0+b1*X1+...+b2*X2 Y = .7309 + 4.5682 X1 + 5.7726 X2 回归系数 b0, b1, b2, ..., b2 .7309 4.5682 5.7726 ----------------------------------------------------- 残差平方和: 2290.69 回归平方和: 2014.61 误差方差的估计 : 63.6303 标准差 = 7.9769 ----------------------------------------------------- 线 性 回 归 显 着 性 检 验 显著性水平 : .050 ----------------------------------------------------- 回归方程整体显著性F检验, H0:b0=b1=...=b2=0 F统计量: 14.5113 F临界值F(2, 33) 3.285 全相关系数 R : .6841 ----------------------------------------------------- 回归系数逐一显著性t检验, H0:bi=0, i=1,...,2 t 临界值 t( 33) 1.6924 回归系数b1-b 2的t值: 5.6411 1.0683 ----------------------------------------------------- 比较残差平方和: 普通最小二乘的: 149400.6000 : 广义最小二乘的: 2290.6910 下面打印的是利用广义最小二乘的回归系数去计算原始资料的回归拟合结果: 要作回归预测吗? 键入 0=不预测, 1=要预测 (0) 回归系数: .7309 4.5682 5.7726 要打印拟合数据吗? 0=不打印, 1=打印 (0) 计算结束。 ------------------------------------------------------------------------------------------------------------------------- 计算中需要一些分析与技巧。因为 需要用它及(4.2.2)计算β* 注意UU′是一个分块对角阵，共有m个对角块，每块是一k×k方阵，元素皆为1。这样Σ的计算就容易了。恰好X也是分成m块，每块为k×p阵。于是 而 ，这就大大简化了计算。如果真的要计算36×36的矩阵Σ的逆，一般微机是不可能的。 这个程序开始时调用了普通多元线性回归程序对原始资料回归。得到的回归方程为 拟合效果也很好(见图4.2.1.2)。但是，这个模型的残差向量的标准差为64.4206，远大于方差分量模型的1.0084。这可以对比说明方差分量模型的作用。 二、最小范数二次无偏估计方法 方差分量模型中方差分量的最小范数二次无偏估计(Minimum Norm Quadratic Unbiased Estimator, MINQUE)是C.R.Rao提出的，思路类似于他处理奇异广义线性模型的分块逆矩阵法，先提出估计应满足的性质，根据这些性质去求解(一般是一个极小值问题)，若能解出，当然就有那些设定的性质。 考虑一般的方差分量模型(4.2.1)，我们要估计方差分量 及其线性函数 （4.2.46） 首先考虑φ的估计应具有的形式，因为是估计方差，可考虑采用二次型Y′AY的形式，即 （4.2.47） A为待估对称矩阵。再来考虑Y’AY应具有的性质： (1)关于参数β的平移不变性。 若参数β有平移： （4.2.48） 则方差分量的估计应该不变。此时原模型变为 （4.2.49） 其二次型估计变为 ，应该有 （4.2.50） 这等价于 （4.2.51） 即平移不变性(4.2.50)应满足的充要条件是 (2)估计量的无偏性。 因为 （4.2.53） 今要求对一切σ2成立： （4.2.54） 则其充要条件应为 （4.2.55） （3)最小范数准则。 经过研究，满足平移不变性与无偏性，尚不能唯一确定待估对称矩阵A。于是可以再加一个优良性质。 如果随机效应向量ξi,i=1,…，m是已知的，则φ=c’σ2的估计应该为 （4.2.56） 这里 （4.2.57） 现在用Y’AY去估计φ=c’σ2，在满足不变性条件AX=0时， （4.2.58） 我们自然希望(4.2.56)与(4.2.58)之间相差很小，这只要求矩阵Δ与U’AU之间相差很小。我们选用矩阵范数‖U’AU-Δ‖来度量Δ与U’AU之间的差异，即应选A使 （4.2.59） 范数可选欧氏范数 （4.2.60） 若记V=UU’，则因Δ是分块对角阵及无偏性要求： （4.2.61） Δ是已知的，则极小化范数 等价于极小化 。 总结上述三项优良性要求，求φ=c’σ2的最小范数无偏估计的问题，归结为求下述极值问题： 因为目标函数是矩阵的迹，所以称为最小迹问题。 要解决极值问题(L1)的解，可以先对其简化。因为 ，V正定而 存在。令 （4.2.62） 则模型(L1)可变为 定理4.2.1 极值问题L2的解为 （4.2.63） 其中 为方程组 （4.2.64） 的解， 为L(Z)的正交补空间L⊥(Z)上的投影阵： （4.2.65） 证明 先证方程组(4.2.64)兼容。 设B0满足极值问题(L2)的约束条件，即 （4.2.66） 则 。因是 往 上的投影阵，故 ，记 ，则 （4.2.67） 引进拉直算符， 表示将Wj按列拉成一个n2×1的向量， 表示将B0按列拉直成一个n2×1的向量。定义n2×m矩阵G， （4.2.68） L(G)表示按G的各列展成的子空间，是n2维。定义μ1为向量 在L(G)上的投影 （4.2.69） μ2为向量 在L⊥(G)上的投影 （4.2.70） 具体形式 ,则 （4.2.71） 这个情形与最小二乘的基本原理图标(图1.2.1)一样，那里是将向量Y向子空间L(X)投影，于是存在常数β1，…，βp,使Y=Xβ=X1β1+…+Xpβp，即将Y表为X列向量的线性组合。现在是必存在常数 ，可将 表示为G的列向量的线性组合，当然也可以将 在L(G)的投影μ1表成线性组合： （4.2.72） 于是 （4.2.73） 现在看(4.2.67)，由公式 得 （4.2.74） 这就证明了 是方程组(4.2.64)的一组解。 再证明(4.2.63)的B*是问题(L2)的解。 由于 故B*满足的约束条件 是从方程组tr(BWi) =Ci,i=1,…，m中解出来的，当然满足这个约束条件。余下看B*是否是trB2的极小值解。设B为任一满足(L2)约束条件的解，记D=B－B*，则D对称，DZ=0,tr(DWj)=0, ,于是 （4.2.75） 因此 （4.2.76） 这就证明了B*是极小值解。 现在回到原问题(L1)。由于(4.2.62)所定义的三个变换都是可逆变换，故(L1)与(L2)等价。于是(L1)的解存在， （4.2.77） 对于最小范数二次无偏估计，我们介绍了它的原理、算法、解的存在性。深入的讨论还可以引出一些问题，如解的非负性，计算的复杂性等，我们这里不再讨论了。 三、极大似然法 在假定方差分量模型随机效应服从正态的情况下，可以使用极大似然法求出参数估计。 设模型为 （4.2.78） （4.2.79） 其它记号仍同以前。则 （4.2.80） Y的条件密度为 （4.2.81） 取对数，略去常数，得似然函数方程 （4.2.82） 用矩阵微商公式 （4.2.83） （4.2.84） （4.2.85） 得似然方程组 （4.2.86） 由于似然方程组没有显式解，统计学家提出了一些迭代算法。但是在数值计算技术发达的今天似乎没有必要再介绍这些迭代算法，甚至连求导的似然方程组都没有必要。对于样本 ，以及已知设计阵 ，代入到似然函数(4.2.82)中，调用本软件所附的计算极值的程序(由Sargent改进的Powell算法)，就可以计算出β，σ2的估计。 第三节 方差分量模型参数的广义岭估计 在这一节我们先将方差分量模型的方差分量化为派生模型的均值参数，分别作出其相对于LSE和BLUE的广义岭估计，再根据二步估计法作出原模型均值参数的广义二乘估计及其进一步的岭估计，证明了这样不仅使方差分量估计的均方误差减少，而且使原模型均值参数估计的均方误差也不增加和进一步减少。我们还找到了岭参数仅仅依赖于样本的估计。这样既将岭估计方法推进至方差分量模型，也改进了方差分量模型参数的离差均值对应方法。 一、方差分量岭估计的构造与性质 岭估计和Stein估计是减少均方误差的行之有效的方法。已有文献将岭估计推广到多元线性回归模型和设计阵是列降秩的情况。本节则试图将它们推广到方差分量模型，以期改进离差－均值对应方法。 设有一般方差分量模型 （4.3.1） 这里 为已知设计阵， 为固定效应向量， 为随机效应向量。假定 记 则模型(4.3.1)可记为 （4.3.2） 按照两步估计法我们先求方差分量σ的估计，为此采取离差均值对应。记 ，于是有模型 （4.3.3） 再记 ，因为 = ，可得派生模型 （4.3.4） 此时原模型的方差分量成了派生模型的均值参数。 对派生模型(4.3.4)用最小二乘法求出均值参数σ的最小二乘估计(LSE)： （4.3.5） 是σ的无偏估计，但由于派生模型是广义线性模型，它未必是σ的最佳线性无偏估计(BLUE) 。当F已知时， (4.3.6) 只是在一些特殊情形研究证明了 。 不难举例说明 可能出现负值分量，也不难看出D’D很容易呈现病态。当X，V1，…，Vm中有两个接近相等时，D’D的列向量就呈复共线。因此需要对 加以改进。对 亦然。 考虑原模型中心化与派生模型(4.3.4)中心化的关系，有 引理4.3.1 若模型(4.3.1)设计矩阵X，U1，…，Um均已中心化，则派生模型(4.3.4)也已中心化。 证明 已知 这里 ,则 故 证毕 我们不妨设D已经中心化、标准化。 当F未知时，对模型(4.3.4)定义σ的关于 的广义岭估计为 （4.3.7） 这里K=diag(k1,…,km)，ki>0，i=1,…,m。P是正交方阵致P’(D’D)P=diag(λ1,…,λm)=Λ。当k1=k2=…=km时，就成为狭义岭估计。我们将看到狭义岭估计虽然能改进 的均方误差，但不一定能保证第二步估计时 的均方误差不增加。 (K)的第i个分量记作 。于是 （4.3.8） 当F已知时，定义σ的关于 的广义岭估计为 (4.3.9) （4.3.10） 这里P’(D’F-1D)P=Λ。自然这里K、P、Λ与(4.3.7)不同。 引理4.3.2 若 、 存在，则 证明 仅证后式，由 定义，注意K>0, 证毕 引理4.3.2说明 、 分别是对 和 的压缩。但这种压缩是模长的压缩，有的分量可能还有伸长。能否使各分量都压缩呢?我们有 引理4.3.3 设 存在，00使(4.3.13)成立，同时有 （4.2.24） 或使(4.3.14)成立，同时有 （4.2.25） 证明 因为 ，故 由引理4.3.3及定理4.3.1，存在K>0，使(4.3.13)成立，且有 。于是 即(4.3.24)成立。同理可证(4.3.14)、(4.3.25)可以同时成立。 证毕 如果我们找不到引理4.3.3的方法，得不到定理4.3.4的结论，那么我们的两步估计法对岭估计而言就将失去意义。如果发生了 ，那么再进一步对β采用岭估计我们就无法判断最终β的估计的均方误差到底是减少还是增大了。现在有了定理4.3.4，我们就可以对广义线性模型(4.3.2)继续采用广义岭估计： (4.3.26) 这里 。 是正交方阵致 。则又存在 ，使均方误差进一步缩小： （4.3.27） 且存在 ， （4.3.28） 利用这种形式，可以作出不依赖于未知参数的 的估计。对 亦然，不再赘述。 第四节 方差分量模型参数经验Bayes估计 在这一节我们将构造误差正态的方差分量模型参数的经验Bayes估计，定义逆拉直运算，利用作者自己曾经证明的多参数指数族EB估计收敛结果，证明所构造的经验Bayes估计的收敛性。 要利用的多参数指数族参数经验Bayes估计收敛速度的结果概述如下： 设t=(t1,…,tp)以及t (1),…,t(n)为来自多元密度f(t)的当前样本和历史样本，它们都是独立同分布的，f(t)和它的偏导数f(r)(t)都是局部有界的，r+2是构造密度核估计正交多项式空间的维数。密度f(t|θ)的参数θ有先验分布G(θ)，它的估计的Bayes风险为R(G)，经验Bayes风险为Rn(G)。我们证明了： 定理1 设1/2≤δ≤1，且 则 一、方差分量模型参数经验Bayes估计的构造 设有方差分量模型 （4.4.1） 这里X是k×p0矩阵，Ui,i=1,…,m是k×pi矩阵。β是p0×1固定效应向量，ξi,i=1,…,m是pi×1随机效应向量。假定ξi互相独立并且 令 这里Vi,i=1,…, m以及Σ都是k×k矩阵，令 这里U是 矩阵，ξ是 向量，则模型(5.4.1)成为 且有 这里Λ是对角阵： 令向量 则模型的基本问题是要估计固定效应β和方差分量σ2。 以下我们要将待估参数β和σ2化为多参数指数族的参数。因为已经假定 ，所以y的条件密度为 由于 我们可以定义参数变换为： 这里θ是p×1向量，p=k(k+1)。取资料变换为 这里t是p×1向量，令 则t的条件密度化成了多参数指数族： 如果我们利用样本变换数据 作出了参数θ的估计，则我们可以根据参数变换的逆变换作出β和Σ的估计，进而作出β和σ2的估计。 至于如何利用 作出θ的经验Bayes估计，可以参看第七章第三节与第十章第四节，这里从略。 二、方差分量模型参数经验Bayes估计的收敛性 对一般的线性模型，当σ2是一个数的时候，作者曾得到β2、σ2的联立经验Bayes估计的收敛速度。现在我们对方差分量模型，σ2是一个向量，无法得到EB估计的收敛速度，但可以推出收敛性。 令 为参数θ的EB估计， 是θ的Bayes估计，由定理1知 于是有 由于均方收敛必定导致依概率收敛，我们就有 现在我们定义拉直运算的逆运算。由于逆拉直运算结果可能不唯一，所以须标注脚标。如果对拉直运算有 这里B是将A的列向量顺次排列而成，则我们定义逆拉直运算为 这里A的列向量是将长向量B截成n等分并顺次排列。如 当有多重运算复合时，我们规定从里层做起。当逆拉直成方阵时，我们省去脚标。如B为n2×1向量，则 表示先将B逆拉直成n×n方阵，再求逆矩阵，再拉直成n2×1的列向量。显然，拉直运算与逆拉直运算结果都是唯一确定的。 对 施行逆拉直运算，得 于是 M的元素 。令 则 是含有k2个方程的矛盾方程组，其最小二乘解为 因为我们已经得到 ，所以 同样我们有 上面两个函数分别对于 和 连续，由于已经有 ，于是有 （4.4.3） 再考虑 的收敛性。注意θ2=Σ-1Xβ和 ，我们有 这是含有k个方程的矛盾方程组，其最小二乘解为 因为 和 已经获得，所以 因为 和 已经获得，所以 同样由于上面两式对θ1和θ2分别连续，并且已经证明有 和 ，我们有 （4.4.4） 这样我们从式(4.4.3)和(4.4.4)就完成了方差分量模型参数EB估计收敛性的证明。 1 35 _1066823475.unknown _1067171623.unknown _1067235352.unknown _1067239812.unknown _1070105604.unknown _1070117438.unknown _1070127769.unknown _1070129427.unknown _1070130069.unknown _1070132538.unknown _1070134497.unknown _1070134813.unknown _1070134332.unknown _1070134159.unknown _1070130542.unknown _1070132480.unknown _1070130491.unknown _1070129827.unknown _1070129890.unknown _1070129613.unknown _1070129297.unknown _1070129374.unknown _1070128101.unknown _1070118352.unknown _1070126448.unknown _1070126898.unknown _1070118987.unknown _1070118037.unknown _1070118244.unknown _1070117945.unknown _1070117200.unknown _1070117332.unknown _1070109304.xls 图表1 64.769 110.1049 46.845 64.3333 41.958 62.7997 16.195 58.2443 79.911 43.8813 91.388 51.8175 33.112 66.7322 8.108 11.9874 271.763 45.1612 381.686 300.2789 14.638 24.1754 68.579 42.7916 -21.871 11.847 22.37 25.0491 11.86 27.5734 179.817 42.0034 236.328 265.3347 -33.122 38.4793 191.666 48.5093 111.135 97.5471 152.015 62.8289 13.821 22.4792 17.892 21.8493 -24.744 71.6181 8.389 59.3433 75.39 32.4214 35.932 35.4476 86.528 29.107 -25.17 10.0247 51.371 53.9348 14.4 59.7473 6.122 20.9349 0.701 13.8346 41.852 30.1161 66.981 104.7942 41.093 117.2629 原始数据 拟合数据 圖4.2.1.1 Sheet1 64.769 110.1049 46.845 64.3333 41.958 62.7997 16.195 58.2443 79.911 43.8813 91.388 51.8175 33.112 66.7322 8.108 11.9874 271.763 45.1612 381.686 300.2789 14.638 24.1754 68.579 42.7916 -21.871 11.847 22.37 25.0491 11.86 27.5734 179.817 42.0034 236.328 265.3347 -33.122 38.4793 191.666 48.5093 111.135 97.5471 152.015 62.8289 13.821 22.4792 17.892 21.8493 -24.744 71.6181 8.389 59.3433 75.39 32.4214 35.932 35.4476 86.528 29.107 -25.17 10.0247 51.371 53.9348 14.4 59.7473 6.122 20.9349 0.701 13.8346 41.852 30.1161 66.981 104.7942 41.093 117.2629 _1070109721.xls 图表1 64.769 115.39 46.845 69.294 41.958 67.7928 16.195 63.0672 79.911 49.2627 91.388 57.4136 33.112 72.1395 8.108 17.0437 271.763 50.6345 381.686 305.7691 14.638 29.5462 68.579 48.3783 -21.871 16.8011 22.37 30.2636 11.86 33.2783 179.817 47.6343 236.328 270.8548 -33.122 44.1755 191.666 54.1219 111.135 103.0804 152.015 68.3822 13.821 27.4711 17.892 26.8342 -24.744 79.156 8.389 64.4431 75.39 37.9677 35.932 40.9021 86.528 34.7794 -25.17 14.9994 51.371 60.2165 14.4 64.864 6.122 26.3636 0.701 18.9678 41.852 35.7149 66.981 109.8582 41.093 122.8459 原始数据 拟合数据 圖4.2.1.2 Sheet1 64.769 115.39 46.845 69.294 41.958 67.7928 16.195 63.0672 79.911 49.2627 91.388 57.4136 33.112 72.1395 8.108 17.0437 271.763 50.6345 381.686 305.7691 14.638 29.5462 68.579 48.3783 -21.871 16.8011 22.37 30.2636 11.86 33.2783 179.817 47.6343 236.328 270.8548 -33.122 44.1755 191.666 54.1219 111.135 103.0804 152.015 68.3822 13.821 27.4711 17.892 26.8342 -24.744 79.156 8.389 64.4431 75.39 37.9677 35.932 40.9021 86.528 34.7794 -25.17 14.9994 51.371 60.2165 14.4 64.864 6.122 26.3636 0.701 18.9678 41.852 35.7149 66.981 109.8582 41.093 122.8459 _1070105647.unknown _1067244009.unknown _1070086070.unknown _1070105485.unknown _1070105559.unknown _1070086703.unknown _1067244429.unknown _1067244635.unknown _1067244724.unknown _1067245017.unknown _1067245101.unknown _1067244898.unknown _1067244680.unknown _1067244465.unknown _1067244234.unknown _1067244392.unknown _1067244149.unknown _1067242083.unknown _1067243873.unknown _1067243956.unknown _1067243813.unknown _1067241422.unknown _1067242054.unknown _1067241307.unknown _1067237963.unknown _1067238552.unknown _1067239518.unknown _1067239724.unknown _1067239761.unknown _1067239607.unknown _1067239272.unknown _1067239413.unknown _1067238621.unknown _1067238267.unknown _1067238363.unknown _1067238450.unknown _1067238324.unknown _1067238126.unknown _1067238189.unknown _1067238071.unknown _1067236153.unknown _1067237484.unknown _1067237700.unknown _1067237849.unknown _1067237660.unknown _1067236803.unknown _1067237089.unknown _1067236549.unknown _1067235878.unknown _1067236015.unknown _1067236057.unknown _1067235961.unknown _1067235496.unknown _1067235811.unknown _1067235438.unknown _1067230432.unknown _1067231783.unknown _1067233265.unknown _1067233880.unknown _1067234096.unknown _1067235041.unknown _1067233935.unknown _1067233534.unknown _1067233656.unknown _1067233425.unknown _1067233377.unknown _1067233402.unknown _1067232461.unknown _1067233214.unknown _1067233237.unknown _1067233006.unknown _1067231926.unknown _1067232220.unknown _1067231822.unknown _1067231362.unknown _1067231551.unknown _1067231614.unknown _1067231687.unknown _1067231598.unknown _1067231487.unknown _1067231520.unknown _1067231452.unknown _1067230712.unknown _1067231145.unknown _1067231246.unknown _1067230974.unknown _1067230596.unknown _1067230627.unknown _1067230527.unknown _1067174743.unknown _1067177097.unknown _1067177506.unknown _1067177790.unknown _1067230356.unknown _1067177605.unknown _1067177219.unknown _1067177368.unknown _1067177129.unknown _1067176623.unknown _1067176962.unknown _1067177008.unknown _1067176921.unknown _1067176212.unknown _1067176587.unknown _1067175178.unknown _1067173960.unknown _1067174407.unknown _1067174601.unknown _1067174638.unknown _1067174465.unknown _1067174301.unknown _1067174349.unknown _1067174262.unknown _1067172462.unknown _1067173350.unknown _1067173803.unknown _1067173237.unknown _1067171777.unknown _1067171831.unknown _1067171697.unknown _1067160268.unknown _1067167179.unknown _1067169110.unknown _1067170791.unknown _1067171066.unknown _1067171480.unknown _1067171557.unknown _1067171327.unknown _1067170924.unknown _1067170965.unknown _1067170897.unknown _1067169675.unknown _1067169788.unknown _1067169978.unknown _1067169731.unknown _1067169187.unknown _1067169635.unknown _1067169111.unknown _1067167888.unknown _1067168455.unknown _1067168542.unknown _1067168722.unknown _1067168532.unknown _1067168300.unknown _1067168372.unknown _1067168019.unknown _1067167404.unknown _1067167691.unknown _1067167755.unknown _1067167468.unknown _1067167304.unknown _1067167350.unknown _1067167277.unknown _1067165515.unknown _1067166453.unknown _1067166870.unknown _1067167038.unknown _1067167125.unknown _1067166776.unknown _1067166814.unknown _1067166840.unknown _1067166685.unknown _1067166104.unknown _1067166174.unknown _1067166285.unknown _1067166153.unknown _1067165665.unknown _1067165860.unknown _1067165633.unknown _1067162823.unknown _1067163529.unknown _1067165428.unknown _1067165481.unknown _1067165130.unknown _1067163096.unknown _1067163143.unknown _1067163013.unknown _1067160817.unknown _1067161062.unknown _1067162739.unknown _1067162675.unknown _1067162711.unknown _1067161544.unknown _1067160984.unknown _1067160661.unknown _1067160730.unknown _1067160376.unknown _1067160514.unknown _1066829304.unknown _1067068043.unknown _1067159493.unknown _1067160035.unknown _1067160191.unknown _1067160229.unknown _1067160116.unknown _1067159825.unknown _1067159890.unknown _1067159802.unknown _1067068457.unknown _1067159282.unknown _1067159451.unknown _1067070065.unknown _1067068230.unknown _1067068336.unknown _1067068148.unknown _1067060518.unknown _1067067464.unknown _1067067700.unknown _1067067908.unknown _1067067504.unknown _1067067300.unknown _1067067332.unknown _1067067225.unknown _1066829674.unknown _1067058869.unknown _1067059122.unknown _1066829803.unknown _1066829553.unknown _1066829630.unknown _1066829526.unknown _1066825947.unknown _1066826660.unknown _1066827569.unknown _1066829174.unknown _1066829269.unknown _1066827700.unknown _1066827373.unknown _1066827519.unknown _1066826690.unknown _1066827371.unknown _1066826283.unknown _1066826464.unknown _1066826517.unknown _1066826385.unknown _1066826149.unknown _1066826238.unknown _1066826120.unknown _1066824378.unknown _1066825078.unknown _1066825339.unknown _1066825485.unknown _1066825300.unknown _1066824679.unknown _1066824925.unknown _1066824512.unknown _1066823739.unknown _1066823906.unknown _1066823976.unknown _1066823832.unknown _1066823607.unknown _1066823674.unknown _1066823563.unknown _1066815317.unknown _1066818631.unknown _1066821403.unknown _1066822025.unknown _1066822606.unknown _1066822896.unknown _1066823331.unknown _1066822742.unknown _1066822420.unknown _1066822454.unknown _1066822103.unknown _1066821796.unknown _1066821875.unknown _1066821936.unknown _1066821814.unknown _1066821697.unknown _1066821741.unknown _1066821667.unknown _1066820363.unknown _1066820693.unknown _1066821140.unknown _1066821347.unknown _1066820893.unknown _1066820567.unknown _1066820640.unknown _1066820404.unknown _1066820007.unknown _1066820291.unknown _1066820330.unknown _1066820175.unknown _1066819500.unknown _1066820000.unknown _1066818824.unknown _1066816734.unknown _1066817805.unknown _1066818199.unknown _1066818523.unknown _1066818568.unknown _1066818344.unknown _1066817923.unknown _1066818036.unknown _1066817859.unknown _1066817336.unknown _1066817510.unknown _1066817755.unknown _1066817432.unknown _1066817048.unknown _1066817235.unknown _1066816838.unknown _1066816207.unknown _1066816647.unknown _1066816687.unknown _1066816699.unknown _1066816670.unknown _1066816387.unknown _1066816428.unknown _1066816226.unknown _1066815926.unknown _1066816161.unknown _1066816182.unknown _1066816025.unknown _1066815657.unknown _1066815784.unknown _1066815418.unknown _1066804602.unknown _1066806401.unknown _1066807592.unknown _1066814663.unknown _1066815156.unknown _1066815187.unknown _1066815094.unknown _1066815122.unknown _1066814915.unknown _1066814509.unknown _1066814562.unknown _1066809427.unknown _1066807155.unknown _1066807452.unknown _1066807531.unknown _1066807308.unknown _1066806978.unknown _1066807131.unknown _1066806613.unknown _1066805543.unknown _1066805893.unknown _1066806151.unknown _1066806318.unknown _1066805966.unknown _1066805712.unknown _1066805852.unknown _1066805640.unknown _1066805021.unknown _1066805247.unknown _1066805465.unknown _1066805076.unknown _1066804748.unknown _1066804840.unknown _1066804641.unknown _1066802380.unknown _1066803072.unknown _1066803816.unknown _1066803956.unknown _1066804065.unknown _1066804535.unknown _1066803874.unknown _1066803304.unknown _1066803732.unknown _1066803243.unknown _1066802760.unknown _1066802919.unknown _1066802938.unknown _1066802814.unknown _1066802513.unknown _1066802721.unknown _1066802469.unknown _1066716575.unknown _1066716957.unknown _1066717088.unknown _1066802232.unknown _1066717020.unknown _1066716764.unknown _1066716930.unknown _1066716685.unknown _1066716153.unknown _1066716456.unknown _1066716525.unknown _1066716309.unknown _1066715995.unknown _1066716026.unknown _1066715851.unknown