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