安徽农业大学理学院：《高级生物统计》课程教学资源（讲义）第六讲 方差分析最优线性无偏估计法 BLUP OF ANOVA

第六讲方差分析的最优线性无偏估计法 BLUP OF ANOVA (本讲不讲) 1974年美国康奈大学C·R· Henderson博士及其协作者首先提出应用线性模型来评价种公牛育种值的 方法。他们将育种值作为观察值的线性函数采用稍做改良的最小二乘分析法,使其估计值的误差为最小。 并定名为最优线性无偏估计( Best liner unbiased predication简称BLUP法)。由于该法具有对估计效应值 的误差最小、精度最优及无偏的特点,并可由计算机进行大量的繁琐计算,故BLUP法在八十年代就被 些先进国家的学者应用于生产实践中,并在方差分析中得到具体的应用 第一节参数的最小二乘估计 、线性模型 假定观察的随机变量Y与若干个因素X1,X2,…,Xp间存在线性关系,其线性模型为: Y=B1x1+B2X2+…+BpXp+ (5-1) 其中β1,β2…βp为待估的参数,ε为随机误差。若对Ⅺ,X2…Xp,Y作n次观察(实验),则(5 1)式可表示为: y1=x1B1+x12B2 XIpBp 1B1+x22B2 2pB. +a2 y=xmB.+x,B, B 将(5-2)式改写成矩阵形式为: Y=BX+E 式中Y为观察值向量,β为参数向量,Ⅹ为设计矩阵或结构矩阵,ε为误差向量。通常假定ε;~N(o, 02)i=1,2,…,n且相互独立,因此,对于(5-3)式可进一步记为(Y,XB,o21),意为观察值向 量Y有 Er=XB Vary=o/ 对线性模型(5-4),主要对β或β的函数及σ2作出估计并进行假设检验 、效应的最小二乘估计 对于β的估计,根据(5-4)式有: B=X′Y 称(5-5)式为正规方程,它的解为β的最小二乘估计,记作B 1、当(X′ⅹ)为可逆(满秩)矩阵时

55 第六讲 方差分析的最优线性无偏估计法 BLUP OF ANOVA （本讲不讲） 1974 年美国康奈大学 C·R·Henderson 博士及其协作者首先提出应用线性模型来评价种公牛育种值的 方法。他们将育种值作为观察值的线性函数采用稍做改良的最小二乘分析法，使其估计值的误差为最小。 并定名为最优线性无偏估计（Best Liner Unbiased Predication 简称 BLUP 法）。由于该法具有对估计效应值 的误差最小、精度最优及无偏的特点，并可由计算机进行大量的繁琐计算，故 BLUP 法在八十年代就被一 些先进国家的学者应用于生产实践中，并在方差分析中得到具体的应用。 第一节 参数的最小二乘估计 一、线性模型 假定观察的随机变量 Y 与若干个因素 X1，X2，…，XP间存在线性关系，其线性模型为： Y=β1X1+β2X2+…+βpXp+ε （5—1） 其中β1，β2…βp 为待估的参数，ε为随机误差。若对 X1，X2…Xp，Y 作 n 次观察（实验），则（5 —1）式可表示为：        = + + + + = + + + + = + + + + n n n np p n p p p p y x x x y x x x y x x x                 1 1 2 2 2 21 1 22 2 2 2 1 11 1 12 2 1 1 （5—2） 将（5—2）式改写成矩阵形式为： Y=βX+ε （5—3） 式中 Y 为观察值向量，β为参数向量，X 为设计矩阵或结构矩阵，ε为误差向量。通常假定εi～N（o， σ2） i=1，2，…，n 且相互独立，因此，对于（5—3）式可进一步记为（Y，Xβ，σ2 I），意为观察值向 量 Y 有    = = VarY I EY X 2   (5—4) 对线性模型（5—4），主要对β或β的函数及σ2 作出估计并进行假设检验。 二、效应的最小二乘估计 对于β的估计，根据（5—4）式有： (X＇X)β=X＇Y （5—5） 称（5—5））式为正规方程，它的解为β的最小二乘估计，记作  ˆ 。 1、当（X＇X）为可逆（满秩）矩阵时

B=( 若要估计β1,β2,…βp的线性函数 P=C1B1+C2B2+…+CpBp (5-7) 则P=c'B 称为P的最小二乘估计。 对于02的估计记 R2=(Y-XB)(-B)=(-BX)Y-B)由(5-6)式有 Ro=YY-YXB=YY-BXr (5-8) 称R2为残差平方和,记r为X的秩,则 R6/ 为σ2的无偏估计 最小二乘估计的性质有: E(B)=EICXX-XY=(XX)XEY=(XX)XXB=B (5-10) Var(B)=Var(X X)XY=(XX)Xo IX(XX) (5-11) (r-X-YY(XX-1=0(Xr) 类似地 E(CB)=CB Var(CP=oIC(XX)] (5-13) 即:B、C'B分别为β、Cβ的唯一最小方差线性无偏估计,简称BLUE( Best linear unbiased estimator 2、当(X′X)不可逆时,(5-5)式无唯一解,B有无穷多解,此时常对β加上一个约束条件 HB=O (5-14) 使B的各分量不独立,从(5-14)可解出B与其中独立的那部分分量间的关系: B=T B(O) (5-15) βo)是B中在约束(5-14)下,独立的那些分量构成的向量(维数<P),于是原模型(5-4)变成 「EY=(XTB0 (5-16) vAry=oI 对于β0而言,又成为无约束问题了,所以上述(5-5)至(5-13)仍可适用,只须将其中X改为 XT,B改为B0),于是有正规方程 (XT)′(XT)B0)=(XT)′Y (5-17) 解出Bo,并由(5-15)式得到B,这时B、CB是在HB=0约束条件下的唯一的BLUE,记Bm0。 、对于参数的假设检验,即检验:H1B=0 Ro/ox(n-r (5-18) 其中r=rank(X),即X阵的秩 2、在约束H1B=0下 RH=(Y-XBH(Y-XBH) (5-19) 可以证明 R.a2~x2(d1) (5-20)

56  ˆ =(X＇X)—1X＇Y （5—6） 若要估计β1，β2，…βp 的线性函数 P=C1β1+C2β2+…+Cpβp （5—7） 则 ˆ  ˆ P = c  C＇=（C1，C2，…Cp） 称为 P 的最小二乘估计。 对于σ2 的估计记： ) ˆ )( ˆ ) ( ˆ ) ( ˆ ( 2 R0 = Y − X  Y − X = Y − X  Y − X 由（5—6）式有 R = YY −YX = YY − X Y 2 ˆ ˆ 0 （5—8） 称 2 R0 为残差平方和，记 r 为 X 的秩，则 ˆ = R / n − r 2 0 2  （5—9） 为σ2 的无偏估计。 最小二乘估计的性质有：  =   =   =    =  − − − E E X X X Y X X X EY X X X X 1 1 1 ) [( ) ] ( ) ( ) ˆ ( （5—10） 1 1 2 1 2 1 1 1 2 1 ( ) ( ) ( ) ) [( ) ] ( ) ( ) ˆ ( − − − − − − − =   =  =   =    X X X X X X I X X Var Var X X X Y X X X IX X X     （5—11） 类似地： E(C ˆ ) = C （5—12） ) [ ( ) ] ˆ ( 2 1 Var C C X X C −  =    （5—13） 即：  ˆ 、C＇ ˆ 分别为β、C＇β的唯一最小方差线性无偏估计，简称 BLUE（Best Linear Unbiased Estimator） 2、当（X＇X）不可逆时，（5—5）式无唯一解，  ˆ 有无穷多解，此时常对β加上一个约束条件 Hβ=0 （5—14） 使β的各分量不独立，从（5—14）可解出β与其中独立的那部分分量间的关系： β=Tβ(0) （5—15） β(0)是β中在约束（5—14）下，独立的那些分量构成的向量(维数<P)，于是原模型（5—4）变成    = = VarY I EY XT 2 (0) ( )   （5—16） 对于β(0)而言，又成为无约束问题了，所以上述（5—5）至（5—13）仍可适用，只须将其中 X 改为 XT，β改为β(0)，于是有正规方程 （XT）＇（XT）β(0) =（XT）＇Y （5—17） 解出β(0)，并由（5—15）式得到  ˆ ，这时  ˆ 、C  ˆ  是在 Hβ=0 约束条件下的唯一的 BLUE，记 0 ˆ  H 。 三、对于参数的假设检验，即检验：H1β=0 1、 2 2 RO / ～ ( ) 2  n − r （5—18） 其中 r=rank(X)，即 X 阵的秩。 2、在约束 H1β=0 下 ) ˆ ) ( ˆ ( 1 1 1 2 RH Y X H Y − X H = −  （5—19） 可以证明 2 2 / 1 RH  ～ ( ) 1 2  d （5—20）

d1=n-[r(x)-r(H1 (R2-R)a2~x2(d2) d,=r(X)-[()-r(HDI 3、在假设H1B=0成立时(5-18)与(5-20)式独立 (R.-R)/d2 (d2,n-r) (5-21) R2/ 在方差分析中虽然常有约束,如各效应之和为0,即H=(1……1),这是个理想的约束。它的特点是秩 ramk()=P=rmk(X)+mmk(H)即满秩,如 21 X 201 H 101 这时有R=R2,仍可用(5-21)检验H1B=0是否成立。 第二节单因素方差分析的最小二乘估计 、数学模型 yj=μ+at+eii=1,2,…,Pj=1,2,…,n 其中y为A因素第i个水平的第j次观察值,a为A的第i个水平的效应,e是随机误差,服从N(O, 2,u为总体平均数。(1A式可写成 (2A) 其中 Y"=(y1…y1,nl…yplr…yp.np),B′=( a 2 e′=(e1…el,nl…"ep…"ep,mp) 设计矩阵X 00 100 011 、效应的最小二乘估计 正规方程组(XX)B=X7 nn (3A) 其中N= . Y y 常把(3A)写成表格形式

57 [ ( ) ( )] 1 1 H1 d n r r X = − H − 2 2 0 2 ( )/ 1 RH − R  ～ ( ) 2 2  d ( ) [ ( ) ( )] 2 1 H1 d r X r r X = − H − 3、在假设 H1β=0 成立时（5—18）与（5—20）式独立 R n r R R d F H − − = / ( )/ 2 0 2 2 0 2 1 ～ ( , ) 2 F d n − r （5—21） 在方差分析中虽然常有约束，如各效应之和为 0，即 H=（1……1），这是个理想的约束。它的特点是秩 rank( ) P rank(X) rank(H) X H = = + 即满秩，如 3 1 0 1 2 0 1 1 2 1 (1 0 1) 2 0 1 1 2 1 3 3 =           = =         =         =  P rank P H X X H rank 这时有 2 0 2 RH = R ，仍可用（5—21）检验 H1β=0 是否成立。 第二节 单因素方差分析的最小二乘估计 一、数学模型 yij=μ+αi+eij i=1，2，…，P j=1，2，…，n； （1A） 其中 yij 为 A 因素第 i 个水平的第 j 次观察值，αi 为 A 的第 i 个水平的效应，eij 是随机误差，服从 N（0， 2  e ），μ为总体平均数。(1A)式可写成 Y=Xβ+ε （2A） 其中 Y＇=（y11…y1, n1…yp1…yp, np）， β＇=（μ，α1，α2…αp） ε＇=（e11…e1, n1…ep1…ep, np） 设计矩阵 p N X +                       = ( 1) 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 1 1 1 1                                  二、效应的最小二乘估计 正规方程组 X X  = X Y ˆ ( )                 =                                   . 2 1 2 1 2 2 1 1 1 2 .. 0 0 0 0 0 0 P P P P Y Y Y Y n np n n n n N n n n               （3A） 其中 i p i N n =1 =  ij n j i Y y i =1  =  （i=1，2，…，P） ij n j p i i p i Y y y i 1 =1 =1  =  =  =   常把（3A）写成表格形式 （4A）

u a a2 a RHM H: N Y ar1:nnt0… 00 (RHM: Right hand member意为右边数) (3A)式的系数矩阵(X′X)是不可逆的,其秩为P。这时须对a1,a2,…ap加上一个约束条件 ∑a1=0(i=1,2,…,P) (5A) (6A) 于是参数B的估计可用以下三种方法 第一种逆矩阵法 该法的最大优点在于求解求逆过程中的有关数值可直接用于假设检验。 为使系数矩阵可逆,利用ap=-(a1+ax+…+ap-1)可化(3A)或(4A)式为同解方程式 如(4A)式中第一个方程和最后一个方程可作变化 Nu+n a 1+n2 a 2++np-la p-1+np(-ai-a 即:Np+(n-n)a计+(n-n)a2+…+(n-1-np)ap-1=Y 而 即;np以 从第二个方程起均减去最后一个方程(7A)式,得 (8A) 以表格示之的正规方程组 &I n -np n tnp np n2 tnp np (8A)的系数矩阵S为可逆矩阵,记S1,则(8A)式有唯一解。 P-12 则B=S-(X)为约束(5A)下的BLUE。实际上不管用何种方法,解得B都是B的BLUE。 第二种简捷法 逆矩阵法须将(4A)化为(8A),其过程较为繁琐。故可直接从(4A)解得β。从该式中的a;方程知

58  ˆ 1  ˆ 2  ˆ …………  P ˆ RHM μ： N n1 n2 ……………np Y．． 1： n1 n1 0 ……………0 Y1．  2： n2 0 n2 ……………0 Y2． ∶ ∶ ∶ ∶ ∶ ∶ ∶ ∶ ∶ ∶  p ： np 0 0 ……………np Yp． （RHM：Right Hand Member 意为右边数） （3A）式的系数矩阵（X＇X）是不可逆的，其秩为 P。这时须对α1，α2，…αp 加上一个约束条件 0 ( 1,2, , ) 1 i i p p i  = =  =  （5A） 即 i p i ap  1 1 −  =  （6A） 于是参数β的估计可用以下三种方法。 第一种 逆矩阵法 该法的最大优点在于求解求逆过程中的有关数值可直接用于假设检验。 为使系数矩阵可逆，利用αp=－（α1+α2+…+αp－1）可化（3A）或（4A）式为同解方程式。 如（4A）式中第一个方程和最后一个方程可作变化 Nμ+n1α1+n2α2+…+np－1αp－1+np（－α1－α2―…―αp－1）=Y．． 即：Nμ+（n1－np）α1+（n2－np）α2+…+（np―1―np）αp－1=Y．． 而 npμ+0α1+0α2+…+np（－α1－α2－…－αp－1）=Yp． 即：npμ－npα1－npα2―…―npαp－1=Yp． （7A） 从第二个方程起均减去最后一个方程（7A）式，得 （8A） 以表格示之的正规方程组  ˆ 1  ˆ 2  ˆ … … 1 ˆ  p− RHM 1 1 1 1. . 2 2 2 2. . 1 1 1 1. . 1 2 1 : : : : .. p p p p p p p p p p p p p p p p p p p p p p p p p n n n n n n n Y Y n n n n n n n Y Y n n n n n n Y Y N n n n n n n Y − − − − + − − + − − − − − − − − −                   （8A）的系数矩阵 S 为可逆矩阵，记 S ―1，则（8A）式有唯一解。               = − − − − − − − − 1 0 1 1 1 2 1 1 10 11 12 1 1 00 01 02 0 1 1 p p p p p p p C C C C C C C C C C C C S         （9A） 则 ( ) ˆ 1 = s X Y −  为约束（5A）下的 BLUE。实际上不管用何种方法，解得  ˆ 都是β的 BLUE。 第二种 简捷法 逆矩阵法须将（4A）化为（8A），其过程较为繁琐。故可直接从（4A）解得  ˆ 。从该式中的αi 方程知

a=y (10A) 其中y=Y/n,为平均数,因∑a=0,故相加得 =y+2 点(P-1)-Y2-…-n 回代(10A) -+(p-)2-…-F (12A) an=[-F-F2…+(P-1) 将A,a1,…,an-1式中的各平均数前的系数用矩阵K表示 P-1-1-1 (13A) 则由该简捷法相配合的S-1为 S-l=KDK′ (14A) 其中D为 第三种参数变换法 对(4A)方程组作如下变换 a (16A) 则方程(4A)变为 比原式少最后一列和一行的 (17A) u aa HM

59        + = + = + =    P YP Y Y       ˆ ˆ ˆ ˆ ˆ ˆ 2 2 1 1  （10A） 其中 Yi Yi ni /  =  为平均数，因∑αi=0，故相加得 ˆ 1 ( ) = Y1 + Y2 + + YP p   （11A） 回代（10A）        = − − − + − = − + − − − = − − − −          ˆ [ ( 1) ] ˆ [ ( 1) ] ˆ [( 1) ] 1 2 1 1 2 1 2 1 2 1 1 p p p p p P p Y Y p Y Y p Y Y p Y Y Y        （12A） 将 1 1 ˆ , ˆ , , ˆ     p− 式中的各平均数前的系数用矩阵 K 表示 P P P P P P K                  − − − − − − − − − − − − − = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1              （13A） 则由该简捷法相配合的 S －1 为 S －1=KD－1K＇ （14A） 其中 D 为 p p np n n D                  =   2 1 （15A） 第三种 参数变换法 对（4A）方程组作如下变换             =  = −  = −  = −  = + − − ˆ 0 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ 1 1 2 2 1 1 p p p p p p p               （16A） 则方程（4A）变为： 比原式少最后一列和一行的 （17A）  ˆ 1  ˆ 2  ˆ ……… 1 ˆ −   p RHM