Transcript 第七章个体遗传评定
第七章 个体遗传评定
—BLUP法
第一节 有关基础知识
第二节 BLUP育种值估计
第三节 遗传参数估计的REML方
法
第一节
有关基础知识
矩阵代数基础
纯量、矩阵和向量
纯量(scalar)
• 只有大小的一个数值,也称为标量、数量
或元向量。
• 用数字或经定义的拉丁字母斜体、小写表
示。
• 如a、r和k。
矩阵(matrix)
• 由一定行数和一定列数的纯量,按一定顺序排
列的表。
• 一般用大写粗体字母表示。
• 矩阵的阶数 (order) 或维数(dimension)是指
矩阵的行数 (m) 和列数 (n) ,表示为m n。
• 例如:
向量 (Vector)
• 仅有一列或一行的矩阵,前者称为列向量
(column vector),后者称为行向量(row vector)。
• 通常用小写粗体字母表示。
• 为区别行向量和列向量,通常在字母的右上
角加一撇表示行向量,不加撇表示列向量。
• 行向量的阶数为1j,列向量的阶数为j1。
• 例如:
1
a 2
3
a 1 2 3
一些特殊矩阵
方阵(square matrix):行数与列数相等的矩
阵,如An×n。其他矩阵称为直角阵(rectangular
matrices) 。
a b c
A c a b
b c a
对称阵 (symmetric matrix):元素间满足 aij =
aji 的方阵。
a b c
A b a b
c b a
0.5 0.8
1
B 0.5
1
0.2
0.8 0.2
1
三角阵(triangular matrix):
• 上三角阵:主对角线以
下元素全部为0,即当 j <
i 时, aij = 0 (j < i)
• 下三角阵:主对角线以
上元素全部为0,即当 j
> i 时, aij = 0 (i < j)
a11 a12
0 a22
a1n
a2n
0
0
ann
a11 0
a21 a22
0
0
an1 an 2
ann
对角阵 (diagonal matrix):除 i=j 时的元素
(主对角线元素)外,其它元素均为零的方阵,
即aij = 0 (j i 时)。通常可以用Diag {aj }表示,
其中aj 为该阵的第 i 个对角线元素。
a11
0
0
0
a22
0
0
0
ann
单位矩阵(identity matrix):所有对角线元素为1,
其他元素均为0的矩阵。
1 0
0 1
0 0
0
0
1
分块矩阵(block matrix):用水平和垂直虚线将
矩阵分为若干小块,此时的矩阵称为分块矩阵。
其中的小块称为子阵(sub-matrix)。
2
1
A 0
0
0
1
2
0
0
0
1 0 1
2 3 0
A1
1 0 0
O
0 1 0
0 0 1
A2
I1
0 0
1 0 0
2 1
1 0 1
0 0 O, 0 1 0 I
A
,
A
,
1
2
1 2
2
2
0
0 0
0 0 1
分块对角阵 (block diagonal matrix):主对角线
上的子阵都为方阵,其余子阵都是零阵的分块阵。
如:
Α11
0
Α
0
0
Α 22
0
0
0
Αkk
稀疏矩阵 (sparse matrix):设矩阵Amn中有 s
个非零元素,若 s 远远小于矩阵元素的总数 (即
s<<m×n),则称A为稀疏矩阵。
矩阵的运算
加法(addition):当矩阵A和B同阶时,有:
A B aij bij
例如: A 4 5 3
B 1 0 2
6 0 2
3 4 1
4
1
5
0
3
2
5
5
5
BA
AB
6 3 0 4 2 1 9 4 3
对于 m n 阶矩阵 A、 B和 C,具有如下性质:
• 闭 合 性:A + B 仍然是一个 m n 阶矩阵;
• 结 合 性:(A + B) + C = A + (B + C);
• 交 换 性:A + B = B + A;
• 加性等同:A + 0 = A;
• 加 性 逆: A + (−A) = 0 。
乘法 (multiplication)
纯量与矩阵相乘:一个纯量 与一个矩阵 A
的乘积是用 去乘 A 的每个元素,表示为A。
A Aij
对于m n 阶矩阵 A 和 B 以及纯量 和 ,
具有如下性质:
• 闭 合 性: A 仍然是一个 m n 阶矩阵;
• 结 合 性:() A = (A);
• 分 配 性: (A+B) =A+B;
(+) A=A+A;
• 等 同 性:1A = A。
矩阵相乘:当矩阵A的列数与矩阵B的行数相等
时,A与B可乘,即
A r c Bcl Cr l cij
c
其中, cij aik bkj
k
例如:
矩阵乘法具有如下性质:
• AB BA
• (AB)C = A(BC)
• (A + B)C = AC + BC;C(A + B) = CA + CB
转置(transposition):矩阵的行与列对调,用A´
或AT表示,即:
Αmn
aij Αnm a ji
矩阵的转置有如下性质:
• 当A为对称方阵时, A´ = A;
• (A´) ´ = A
• (AB) ´ = B ´A ´
• (AB ´ C) ´ = C ´BA ´
• (A + B + C) ´= A ´ + B ´ + C ´
矩阵的迹 (trace):一个方阵的迹为其对角线元
素之和,表示为:
tr Α aii
设:
0.51 0.32
A 0.28 0.46
0.21 0.16
i
0.19
0.14
0.33
tr Α 0.51 0.46 0.33 1.30
迹的运算性质:
tr ΑΒ tr ΒΑ
tr ABC tr BCA tr CAB
tr aq tr qa qa
则
范数 (norm):矩阵与其转置矩阵乘积的迹的
平方根,即:
Α tr ΑΑ
0.5
2
aij
i j
范数的性质:
• ||A|| > 0,除非A = 0; ||A|| = 0 ⇐⇒ A=0
• ||kA|| = |k| ||A||(k为一纯量);
• ||A+B||≤||A||+||B||
• ||AB||≤||A|| ||B||
0.5
逆矩阵 (inverse matrix):对于一方阵A,若
存在另一矩阵B,使得AB=BA=I,则称B为A的
逆矩阵,并表示为A–1,即A–1 A=I。
对任意 n 阶矩阵A , 称
1
1
A
A*
A1n
A11
| A|
其中,A*是A的伴随矩阵。
A–1存在的先决条件:
A*
A
n1
ann
为A 的伴随矩阵。其中,Aij
是A 中元素aij的代数余子式。
• A必须是一方阵;
• A的行列式|A|0,即A为非奇异阵。
A–1具有如下性质:
• A–1A=AA–1=I;
非奇异阵也就是满秩矩
阵 :对于方阵 A ,如果
存在一个同阶的方阵B ,
两方阵的积为单位阵,
则称方阵 A 为满秩方阵或
非奇异阵。
• A–1是唯一的;
1
1
• Α
;
Α
• (A–1) –1 =A,因而也是非奇异阵;
•
•
•
如果 n 阶矩阵A
的 行 列 式
(A–1) ´ =(A´) –1;
│ A│≠0 , 则 称
如A为对称阵,则A–1也是对称阵; A是非奇异阵;
否则,称A为奇
若A、B均可逆,则(AB)–1 =B–1 A–1 。 异阵。
广义逆 (generalized inverse):对于任一矩
阵A,若有矩阵G,满足:
AGA = A
则称G为A的广义逆,记为A¯,即
AA –A = A
广义逆的性质:
• 若A为方阵且满秩,则A¯ = A –1 ;
• 对于任意矩阵A, A¯必存在。
Kronecker乘积(或直积, direct product)
设 Α aij 和 Β bij 分别为 mn 和 pq 阶
矩阵,定义 C aij Βmpnq ,称为 A 和 B 的
Kronecker乘积或直积,记为 C A B 。
即:
a11 B a12 B
a21 B a22 B
AB
am1 B am 2 B
a1n B
a2 n B
amn B
直积的有关性质:
• 0A=A0=0
• (A1+A2)B=(A1B)+(A2B)
• A(B1+B2)=(AB1)+(AB2)
• ( A)( B)= (AB) (、 均为纯量)
• (A1B1)(A2B2)=(A1A2) (B1B2) (如A1与A2、
B1与B2可乘)
• (AB)′=A′B′
• (AB)-1=A-1B-1 (如A、B均可逆)
• (x′ y)′=y x′=yx′
Hadamard乘积:两个矩阵A和B的元素间
相乘,要求A和B同阶。用*表示:
Α *Β aij * bij abij
随机变量的期望值和方差
设X为一随机变量,则:
其数学期望: E ( X ) ,且具有如下性质:
E(kX ) kE( X ) k (k 为一常数)
E( X Y ) E( X ) E(Y ) X Y (X、Y相互独立时)
E( XY ) E( X ) E(Y ) X Y
其方差:Var ( X ) E( X )2 , 且具有如下性质:
Var (kX ) k 2Var ( X ) k 2 E( X )2
Var ( X Y ) Var ( X ) Var (Y ) (X、Y相互独立时)
Cov( X , Y ) E( X X ) E(Y Y ) (称为X和Y的协
方差)
将上述内容推广至多个变量。设x1, x2,
xn为 n 个随机变量,则:
x1
x x2
xn
随机向量
1
2
E ( x)
n
12 12
2
Var (x) V 21 2
n2
n1
x的期望
向
量
1n
2n
2
n
V的对角线元素为n
个变量的方差;非
对角线元素为变量
间的协方差。
x 的方差-协
方差矩阵
动物育种中常用的是表型方差-协方差矩
阵V、遗传方差-协方差矩阵G和残差方差-协
方差矩阵R。
p2 p
a2 a
p
a
2
2
p
p
a
a
V p
G a
1
12
2
21
pn1 pn 2
1n
2
pn
2n
e21 e12
2
e2
R e21
en1 en 2
1
12
2
21
an1 an 2
e
e
2
e
1n
2n
n
1n
2
an
2n
G阵的构建需要一个由个体间亲缘相关系数组
成的矩阵A,该矩阵称为加性遗传相关矩阵。由
于它是由Wright计算近交系数公式中的分子计
算而得,故又称为分子血缘相关矩阵。
a11 a12
a21 a22
A
个体 i 和 j 间的 a
an 2
n1
加性遗传相关。
计算公式:
1
aij ( )
2
n1 n2
a1n
a2 n
ann
个体 i 的近交
系数加1。即:
aii 1 fi
(1 f A )
A阵元素计算机计算的递推公式:
aii 1 0.5asi di
aij 0.5(ais j aid j )
式中:si 和 d i 分别为个体 i 的父亲和母亲;
s j和 d j分别为个体 j 的父亲和母亲;
asi di 为 si 和 d i 间的加性遗传或亲缘相关系数。
若个体 i 的一个亲本或双亲未知,
asi di 。
0
ai s j 和 ai d j 分别为个体 i 与 s j 和d j 间的加性
遗传相关。s j 未知时,ai s j 0 ;d j未知时,
aid 0 。
j
线性模型基础
模型 (model)
模型:指描述观察值与影响观察值变异性的各
个因子间关系的方程式。
因子:影响观察值的因子也称为变量
• 变量可分为离散型变量(变异不连续)和连
续型变量;
• 离散型因子分为固定因子和随机因子两类;
• 连续性变量通常作为协变量看待。
固定因子与随机因子:与抽样和目的有关
• 固定因子(fixed factors):抽取因子的若干特定
水平、水平数相对较少、研究目的是要对这些水平
的效应进行估计或比较。
• 随机因子(random factors):因子的各水平是其
所有水平的一个随机样本、水平数相对较多、研究
目的是要对该样本去推断总体。
固定效应与随机效应
• 固定效应 (fixed effects):固定因子各水平对观
察值的效应。
• 随机效应 (random effects):随机因子各水平对
观察值的效应。
线性模型 (linear model)
定义:模型中所包含的各个因子是以相加的
形式影响观察值,即它们与观察值的关系为线性
关系。
组成:一个完整的模型应包括3部分内容:
• 数学方程式(或数学模型式)及其解释;
• 模型中随机变量的数学期望、方差协方差;
• 建立模型时的所有假设和约束条件。
例7.1:
yij ai eij
模型中每个参数的解释
•
•
•
•
yij :第 i 个日龄组的第 j 头肉牛的体重, 为观察值;
:总体均值,是一常量;
ai :第 i 个日龄组的效应,为固定效应;
eij :随机误差或残差效应。
随机变量的期望、方差及协方差
E (eij ) 0
E ( yij ) ai
Var( yij ) Var(eij )
Cov(eij , eij ) Cov(eij , eij ) Cov(eij , eij )
2
约束与假设
•
•
•
•
所有犊牛都来自同一个品种;
母亲年龄对犊牛体重无影响;
犊牛的性别相同或性别对体重无影响;
除日龄组外的其他环境条件相同。
对每一观察值建立方程式
观察值
1
个
日
2
龄
体
组
3
4
日龄组
1
2
3
4
y11=198= + a1
y12=204= + a1
y13=201= + a1
y21=203=
+ a2
y22=206=
+ a2
y23=210=
+ a2
y31=205=
+ a3
y32=212=
+ a3
y33=216=
+ a3
y41=225=
+ a4
y42=220=
+ a4
残
差
+e11
+e12
+e13
+e21
+e22
+e23
+e31
+e32
+e33
+e41
+e42
y11 =198= + a1
+ e11
y12 =204= + a1 a1
+ e12
y13 =201= + a1
+ e13
y21 =203=
+ a2
+ e21
y22 =206=
+ a2 a2
+ e22
y23 =210=
+ a2
+ e23
y31 =205=
+ a3
+ e31
y32 =212=
+ a 3 a3
+ e32
y33 =216=
+ a3
+ e33
y41 =225=
+ a4
+ e41
y41 =220=
+ a4 a4 + e42
y
a
e
关联矩阵:又称设计矩阵或发生矩阵。指示 y 中
的元素与 a 中元素的关联情况。
a1 a2 a3 a4
y11
1 1 0 0 0
y12
1 1 0 0 0
y13
1 1 0 0 0
y 21
1 0 1 0 0
a1
y 22
1 0 1 0 0
a2
y 23
1 0 1 0 0
a3
y 31
1 0 0 1 0
a4
y 32
1 0 0 1 0
y 33
1 0 0 1 0
y 41
1 0 0 0 1
y 42
1 0 0 0 1
y
a
X
模型的矩阵表示:
y = Xa + e
式中,y为观察值向量,a为固定的日龄组向量,
e为随机残差效应向量,X为a的关联矩阵。且有:
E (e) 0
E (y) Xa
Var (y) Var (e) I
2
其中,I 为单位矩阵,σ2为观察值的方差。
线性模型的分类
按模型中各因子的性质分类如下:
固定效应模型 (fixed effect model):模型中除
随机残差外,其余所有效应均为固定效应。
y Xb e
随机效应模型 (random effect model):模型
中除外,其余的所有效应均为随机效应。
y 1μ Zu e
混合效应模型 (mixed effect model):模型中除
和e外,既含有固定效应,也含有随机效应。
y Xb Zu e
式中:y—观察值向量;
b—固定效应(包括)向量;
u—随机效应向量;
e —随机残差效应向量;
X—固定效应的关联矩阵(设计矩阵或发生
矩阵);
Z—随机效应的关联矩阵(设计矩阵或发生
矩阵)
第二节 BLUP育种值估计
背
景
• BLUP法是基于克服传统选择指数法的缺点的。
• 选择指数实质上就是育种值的估计值。
• 选择指数法假定不存在影响观察值的系统环境效应,
或者这些效应已知,可以对观察值进行事先校正,则
选择指数是育种值的最佳无偏估值。但是这一假定几
乎不能成立。
• BLUP的基本思路是在估计育种值的同时对系统环境
效应进行估计和校正。
• 根据这一思路,BLUP法必须基于混合模型。
基本原理
BLUP 的 涵 义 : BLUP 是 Best Linear
Unbiased prediction的首字母缩略词,既最佳
线性无偏预测。其中:
• 最佳 (Best):估计误差方差 Var A Aˆ 最小;
• 线性 (Linear):估计值是观察值的线性函数;
• 无偏 (Unbiased):估计值无偏,即估计值的期望值
就是真值,E A Aˆ 0 ;
• 预测 (prediction):是可以对随机效应进行预测。
通常,对固定效应称估计,对随机效应称预测。
混合模型(Mixed model)
y Xb Zu e
式中,y—观察值向量;b和u分别为固定效应和随
机效应向量;e为随机残差向量;X和Z分别为b和
u的关联矩阵。且有:
这里是假定
固定效应与
Eu 0,Ee 0,Ey Xb;
随机效应间
Var u G,Var e R,Cov u,e 0 无互作。育
Var y Var (Zu e)
种中是假定
Var (Zu) Var (e) 2Cov(u, e)
遗传与环境
ZGZ R
间无互作。
V
Co v y, u Cov((Zu e), u) Cov(Zu, u) Cov(e, u) ZG
混合模型方程组(MME)
1
1
ˆ
XR 1 X
XR Z
X R y
b
(1)
1
1
1
1
Z R X Z R Z G uˆ Z R y
系数矩阵 C
解向量 s
右手项 r
MME的解:
ˆb XV 1 X XV 1y,即bˆ 是一个广义最小二乘估值;
ˆu G ZV 1 y Xbˆ
注意:尽管MME的解中有V-1,但MME中没有,直
接求解MME便可,即 sˆ C1r。
XX
C
C-1可用分块矩阵表示为: C
CZX
1
CXZ
ZZ
C
单性状无重复观察值的动物模型BLUP
模型表达式
r
一个个体时: y b j a e
j 1
n 个个体时:
y Xb Za e
y为个体的观
察值;bj为第
j个固定效应;
a为个体的随
机加性遗传
效应;e为随
机残差效应。
式中,y= 表型观测值向量
b= 为固定环境效应向量
X= 固定环境效应b的关联矩阵
a= s个个体的随机加性遗传效应向量
Z= 随机遗传效应a的关联矩阵,当a中的所
有个体都有观测值时,即s=n时,Z = I
e= 随机残差效应向量
随机效应a和e的数学期望为:
E (a) E(e) 0
随机效应a和e的方差-协方差矩阵分别为:
Var (a) G= A
2
a
Var (e) R I
2
e
其中: a2 为加性遗传方差; e2 为随机残差方差;
A为待估个体间的加性遗传相关矩阵 (additive
genetic relationship matrix),即分子亲缘相关系
数矩阵(numerator relationship matrix)。
混合模型方差组(MME)
'
ˆ
X X
X Z b X y (2)
Z'X Z' Z A 1k aˆ Z' y
'
'
式中,A-1 = A的逆矩阵
X′= X 的转置矩阵
Z′= Z 的转置矩阵
与(1)式相比,
(2)式是在(1)
式两边同除了一个
R-1变换而来。
2
1
/
1 h
k
2
2
2
2
a
a / y
h
2
e
2
a
2
y
2
a
2
a
2
y
MME的求解
-1
'
ˆ
b X X
X
Z
X
y
aˆ Z'X Z'Z A1k Z' y
ˆ 为固定效应的BLUE值; aˆ 为随机效应的
其中,b
'
'
BLUP值。
-1
X X
XZ
Z'X Z'Z A1k 为系数矩阵的逆矩阵。
'
'
也可表示为:
C
CZX
XX
C
CZZ
XZ
个体育种值估计的准确度 (Accuracy, ACC)
用估计育种值与真实育种值之间的相关系数 ra aˆ
i i
来度量,其计算公式为:
当个体 i 为非近交个体时:
2
2
2
d
Cov(ai , aˆi )
aˆ
aˆ
a
a e
ACCi ra aˆ =
1 da k
2
a a
a aˆ a
a
i
i
i
i i
i
i
i
i
式中,dai 为CZZ 中与个体 i 对应的对角线元素。
当个体 i 为近交个体时:
ACCi rai aˆi 1
d ai k
1 fi
( a2i (1 fi ) a2 )
式中,fi 为个体 i 的近交系数。
Relationship
between true BV
and EBV for
three accuracy
values, r= 0.8, 0.5
and 0.3
对育种值估计准确度的进一步理解
• 个体某一性状的真实育种值只有一个,而且永远
不变。
• 个体的估计育种值则与信息来源和信息的多少密
切相关,并随之变化。
• 一般而言,利用的信息越多,准确度越大。
• 准确度还是育种值估计中利用信息多少的一个度
量。它表明了当拥有更多信息时,EBV变化的可
能性。
• 随着准确度的提高, EBV的变化减小。
估计育种值的可靠性 (Reliability)
估计育种值与真实育种值之间的相关系数 ra aˆ 的平
i i
2
方 ra aˆ 称为估计育种值的可靠性。
i i
可靠性是一个决定系数,是对真实育种值变异中
2
2
由估计育种值说明的变异部分( ra aˆ a )的一个度
i i
量。
预测误差方差(Prediction error variance, PEV)
PEV是育种值预测时对误差大小的一个度量。
PEV是对真实育种值变异中未由估计育种值说明
的变异部分 ( (1 raa2ˆ ) a2 ) 的一个度量。
当个体 i 为非近交个体时:
PEVai=Var(ai aˆi )=dai e2 (1 raa2ˆ ) a2 dai k a2
当个体 i 为近交个体时:
PEVai=dai e2 (1- ra2i aˆi )(1 fi ) a2 dai k (1 fi ) a2
预测误差标准差 (Standard error of prediction, SEP)
SEPai= d ai k (1 fi ) a2
育种值估计的准确度可由预测误差方差计算
PEV
ACC 1
(1 f ) a2
单性状有重复观察值的动物模型BLUP
模型表达式—重复力模型 (Repeatability model)
y = Xb + Z1a + Z2p + e
式中,y为观察值向量;b为固定效应向量;a为随
机加性遗传效应向量;p为随机永久环境效应向量;
e为随机残差效应向量。Z1为加性遗传效应的关联
矩阵;Z2为永久环境效应的关联矩阵。且有:
E(a) 0 E(p) 0 E(e) 0 E (y) Xb
2
2
2
Var
(
p
)
I
Var (e) I e
Var(a) A a
p
混合模型方程组
XZ1
XZ 2 bˆ Xy
XX
Z X Z Z A 1k
ˆ
Z1 Z 2 a Z1 y
1 1
1
1
Z2 X
Z2 Z1
Z2 Z 2 Ik 2 pˆ Z2 y
其中:
1 h
k1
2
h
2
e
2
a
2
1 re
k2
2
re h
2
e
2
p
其中:h2为性状的遗传力,re为性状的重复力。
动物模型BLUP的理想性质
• 能最有效地充分利用所有亲属的信息;
• 能校正由于选配所造成的偏差;
• 当利用个体的重复记录时,可将由于淘汰造成的
偏差降到最低;
• 能考虑不同群体及不同世代的遗传差异;
• 能提供个体育种值的最精确的无偏估值。
达到如上理想性质的前提是:
数据特别是系谱正确完整;
所用模型是真实模型;
随机效应方差组分已知,而且是正确的。
多性状动物模型BLUP 育种值估计
模型表达式
以两个性状为例。
X1b1 Z1a1 e1
性状2的模型为: y 2 X2b2 Z2a2 e2
性状1的模型为: y1
y1
a1
e1
令:y y a a e e
2
2
2
X1 0
Z1 0
X
Z
则有:
0 X 2
0 Z 2
y Xb Za e
随机效应的期望和方差-协方差
E (a) E(e) 0
E(y) Xb
A
A a12
Var (a) G=
2
A
A
a
a2
12
I e21 I e12
Var (e) R=
2
I
I
e
e2
12
2
a1
其中, a2 和 a2 分别为性状1和2的加性遗传方差;
2
2
a 为性状1和2间的加性遗传协方差。 e 和 e 分
别为性状1和2的残差方差;
e 为性状1和2间的残差
协方差。 A为分子血缘相关矩阵;I为单位矩阵。
1
2
1
12
12
2
MME及其求解
1
ˆ
XR X
XR Z b XR y
ZR 1 X ZR 1Z G 1 aˆ ZR 1y
1
1
计算G-1和R-1后,代入混合模型方程组求解,便
可得到:
ˆb b1
b 2
a
ˆa 1
a 2
固定效应的BLUE值
加性遗传效应的BLUP值
第三节 遗传参数估计的REML方
法
最大似然法的一般原理
最大似然 (Maximum likelihood, ML) 法 是统计
学参数估计的一个重要方法, 由Fisher提出。
概率函数与密度函数
设f (x, ) 为随机变量X的概率函数(X为离散型变
量)或密度函数(X为连续型变量), 为有关的
已知参数,且X为 的函数。则f (x, )为已知时,
X=x的发生概率(X为离散型变量)或X=x的概率
密度(X为连续型变量)。
似然函数(likelihood function)
假定 为待估计的未知参数,但 X=x 已知,为观
察值。将f (, x) 看作是x一定时 的一个函数,称
之为似然函数,记为L(, x),表示 取不同值时,
X=x发生的可能性或似然性。
最大似然原理
寻找一个使得X=x发生可能性最大的 值ˆ ,即当
X=x时,寻找 ˆ,满足 L(ˆ, x) max{L( , x) : }
其中为参数空间。
如果 ˆ 存在,则称它为 的最大似然估值。
注意事项
: 的估计值必须落在参数空间内。我们
不考虑参数空间外的ML估值。
似然函数通常用对数形式,即ln L或log L。因
为在求最大值过程中需要求对数后求导。
混合模型下参数的ML估计
ML法要求说明数据的分布。对于混合模型,通常
假定呈多元正态分布。似然函数的对数为:
1
1
log L(b, σ y ) C log V ( y Xb)V 1 ( y Xb)
2
2
其中C为常数。根据定义,b和σ的ML估值是参数
空间 中的一个使上式最大化的一个值。
ML法的缺点
对于有限样本,ML估值有偏。有偏来源主要
有两个:
ML估值被严格地限定在参数空间中,不考
虑参数空间之外的具有较大似然值的解;
ML估值不考虑因估计固定效应而造成的自
由度的损失。
约束最大似然法的一般原理
基本思路
用 n r ( X)个线型独立的误差对照 Ky 的约束似然
函数 Lr 代替全似然函数 L 。其中,n为有记录的
个体数,r ( X)为固定效应关联矩阵的秩。
设混合模型为: y Xb Za e
且有:
a ~ N (0, A ) e ~ N (0, I )
2
a
2
2
y ~ N (Xb, ZAZ a I e )
2
e
对于y,总可以找到一个线性函数 Ky ,它满
足 K X 0 ,使得:
E (Ky ) KXb 0
亦即Ky 不受固定效应影响。
若y服从正态分布,则 Ky 也服从正态分布,
即:
Ky ~ N (0, KVK )
Ky 的对数似然函数,也即约束似然函数的
对数为:
1
1
log Lr (σ y ) C log K VK y K (K VK )-1 K y
2
2
其中,C为一常数,且有:
C 0.5(n r (X)) log 2
n 为有记录的个体数,r ( X) 为固定效应关联
矩阵的秩。
求解上述对数约束似然函数关于 a2 和 e2 的
2
2
最大值,便可得到 a 和 e 的REML估值。
于是:
h
2
2
a
a2 e2
REML和ML的区别
思路上:ML是求观察值的似然函数的最大值;
而REML是求观察值的某一特定线性函数的似然
函数的最大值。
结果上:由于 K y不考虑估计固定效应,因此,
REML法克服了ML法中产生有偏估值的一个来
源,即克服了ML法中因估计固定效应而造成的
自由度的损失。
REML和ML的相同点
都具有大样本特性,即小样本数据估值有偏。
REML的有关算法与相关软件
非求导的REML—DFREML
DFREML、MTDFREML、DMU
求一阶偏导数的REML—EMREML
VCE、REMLF90
利用一阶和二阶偏导数的REML—AIREML
DFREML、AIREMLF90、VCE、DMU