全 文 :
第 26 卷 第 6 期 作 物 学 报 V o l. 26, N o. 6
2000 年 11 月 A CTA A GRONOM ICA S IN ICA N ov. , 2000
数量性状分离分析中分布参数估计的 IECM 算法X
章元明 盖钧镒
(南京农业大学大豆研究所, 农业部国家大豆改良中心, 江苏南京 210095)
提 要 在主基因+ 多基因混合遗传分析中, 随着模型的扩展, 估计成分分布参数的 EM 算法显示其
局限性。本文在ECM 算法和剖分成分分布方差为主基因、多基因及环境三种方差组分基础上, 推演出
S 个CM 步骤的一般迭代公式, 称为迭代 ECM 算法 (简称 IECM 算法)。文中给出利用个别分离世代鉴
定主基因和多基因存在, 以及利用联合多个世代分离分析的 IECM 算法。用 T urbo C+ + 语言编写了所
有计算程序。最后给出的实例说明 IECM 算法比 EM 算法更易收敛。
关键词 IECM 算法; 混合模型; 参数估计
The IECM A lgor ithm for Estima tion of Com ponen t D istr ibution
Param eters in Segrega ting Ana lys is of Quan tita tive Tra its
ZHAN G Yuan2M ing GA I Jun2Y i
(S oy bean R esearch Institu te, N anj ing A g ricu ltu ra l U niversity ; N ationa l Cen ter of S oy bean Imp rovem en t, M inistry of
A g ricu ltu re, N anj ing , 210095)
Abstract Based on bo th the ECM algo rithm and the sp lit of variances of com ponen t
d ist ribu t ion s of the m ix tu re m odel in to their m ajo r gene com ponen t, po lygene com ponen t and
environm en ta l com ponen t, the itera ted ECM ( IECM ) algo rithm w as suggested to est im ate
the d ist ribu t ion param eters in sam p le likelihood funct ion fo r m ajo r gene p lu s po lygene m ixed
inheritance ana lysis in the rep lacem en t of EM algo rithm w h ich w as though t to be no t
sufficien t w hen the genet ic m odels get t ing com p lica ted. T he genera l itera ted fo rm u las in CM i
step s of IECM algo rithm fo r est im at ing the above dist ribu t ion param eters w ere g iven to
iden t ify the ex istence of m ajo r genes and po lygenes fo r sing le segrega t ion popu la t ion, and fo r
jo in t ana lysis of m u lt i2genera t ion popu la t ion s. A ll com pu ter p rogramm es w ere w rit ten in
T u rbo C+ + , and then com p iled and linked in to execu t ive files. T he resu lt from an exam p le
show ed tha t the convergence by u sing IECM algo rithm is easier and bet ter than tha t by u sing
EM algo rithm.
Key words Itera ted ECM algo rithm ; M ix tu re m odel; Param eter est im at ion
数量性状分离分析是建立在混合分布理论基础上的, 它将分离群体分布看作为多个主基
因型受多基因和环境修饰所形成的多个正态分布的混合分布[ 1~ 7 ]。因此, 研究混合分布中成
分分布参数的估计方法十分重要。一般采用极大似然法估计其参数。自D em p ster 等 (1977)
提出 EM 算法后[ 8 ] , 有关算法的研究有较多报道[ 9~ 15 ]。但是, 在数量性状分离分析中常用
X 国家 863 项目和重庆市科委应用基础研究项目
收稿日期: 1998211208, 接受日期: 1999202223
EM 算法估计分布参数[ 1~ 7 ]。这时, 对样本似然函数求偏导常忽略成分分布方差中同时含有
一阶遗传参数和二阶分布参数的项使参数估计相对简化[ 2~ 7, 9, 10 ]。在 2 对主+ 多基因混合遗
传模型中, 这种项占大多数, 将这些项都忽略显然会使参数估计不可靠。若采用 ECM 算法,
也因分离群体成分分布方差不全等使估计环境方差和多基因方差组分时会出现解高达十几次
方的高次方程的情况, 若用 Gau ss2N ew ton 法解这些高次方程也因其解不唯一可能造成迭代
发散。本文根据数量遗传学原理, 在将成分分布方差剖分为环境、多基因和主基因 (由一阶遗
传参数表示)三种方差组分和 ECM 算法的基础上, 推导出了CM i 步骤中估计一阶分布参数、
多基因方差和环境方差组分的迭代公式, 有效地避免了采用 EM 和 ECM 算法的缺陷, 避免
了解高次方程的解不唯一性, 降低了参数估计的维数, 缩短了计算时间, 这称为迭代 ECM
(简称 IECM ) 算法。它有效地解决 2 对主+ 多基因混合遗传分析的主基因存在和多基因存在
的鉴定和联合多个世代分离分析的分布参数的极大似然估计。
1 主+ 多基因混合遗传模型
本文涉及的数量性状主+ 多基因混合遗传的数学模型和有关符号与文献[2~ 5 ]一致。若
x 是分离群体家系数量性状平均数, 则有 e~N (0, R2e ö n) , n 为家系内观测植株数。
2 2 对主基因时分离群体成分分布的方差组分
假定亲本主基因型分别为A A B B (P 1) 和 aabb (P 2) , 则 F 2 和 F 2∶3群体分别为A A B B 、
A A B b、A A bb、A aB B 、A aB b、A abb、aaB B 、aaB b 和 aabb 9 种主基因型及其衍生家系按孟德
尔分离比的混合; B 1 和B 1∶2群体分别为A A B B 、A A B b、A aB B 和A aB b 4 种主基因型及其衍
生家系等比例的混合, B 2 和B 2∶2群体分别为A aB b、A abb、aaB b 和 aabb 4 种主基因型及其衍
生家系等比例的混合。若主基因服从加性2显性2上位性模型, 根据文献[4, 5 ]的假定, F 2 和
F 2∶3群体都为 9 个正态分布的混合, B 1、B 2、B 1∶2和B 2∶2群体都为 4 个正态分布的混合。
若数量性状受主基因控制同时有多基因的修饰, 由数量遗传学原理可知B 1、B 2 和 F 2 群
体的各自成分分布方差分别相等, 可剖分为多基因和环境两种方差组分, 分别记为 R2j0和 R2e;
B 1∶2、B 2∶2和 F 2∶3群体的各自成分分布方差不全等, 多数可剖分为主基因 (由一阶遗传参数表
示)、多基因和环境三种方差组分。记 F 2∶3群体 9 个成分分布的平均数与方差分别为 Lj 和 R2j
( j = 1, ⋯, k 3) , 多基因方差组分为 R2t0, 误差方差为 R2e。假定家系间附加的非遗传变异 E b 为
0, 则 F 2∶3的家系平均数方差分别为:
R21 = R23 = R27 = R29 = R2t0 + R2e ö n
R22 = R2t0 + R2e ö n + [1 ö 2 (d b + i) 2 + 1 ö 4 (hb + j ab) 2 ] ö n
R24 = R2t0 + R2e ö n + [1 ö 2 (d a + i) 2 + 1 ö 4 (h a + j ba) 2 ] ö n
R25 = R2t0 + R2e ö n + [d 2a + d 2b + i2 + (d a + j ab) 2 + (d b + j ba) 2 + (ha + 1 ö 2l) 2
+ (hb + 1 ö 2l) 2 + 1 ö 4l2 ] ö (4n)
R26 = R2t0 + R2e ö n + [1 ö 2 (d a - i) 2 + 1 ö 4 (h a - j ba) 2 ] ö n
R28 = R2t0 + R2e ö n + [1 ö 2 (d b - i) 2 + 1 ö 4 (hb - j ab) 2 ] ö n
(1)
B 1∶2和B 2∶2群体成分分布方差的剖分公式参阅上述公式得到。关于 1 对主基因的情形可参阅
文献[2~ 5 ]。
007 作 物 学 报 26 卷
3 鉴定数量性状主+ 多基因混合遗传模型的 IECM 算法
3. 1 ECM 算法[ 12 ]
ECM 算法是 EM 算法的拓展, 分为 E 步骤和CM 步骤两步。E 步骤与 EM 算法的 E 步
骤是一致的, 见文献[7 ], 这里从略。第 t 次迭代的 CM 步骤是分 S 步进行Q (HßH( t) ) 的极大
化。G = {g s (H) ; s= 1, ⋯, S }是参数估计前选择的 S 个 H的函数, 这里 g 1 (H)是似然函数中除
分布平均数外的参数, g 2 (H) 是除多基因方差组分外的参数, g 3 (H) 是除环境方差组分外的参
数。若多基因方差组分不存在则S = 2, 否则S = 3。在第 t+ 1 次迭代中, 首先进行 E 步骤, 然
后进行 S 个CM 步骤。对于 s= 1, ⋯, S , 在 g s (H) = g s (H{ t+ (s- 1) ö S })和 H参数空间极大化Q (Hß
H( t) ) 以获得 H的条件极大似然估计值 H( t+ s ö S ) , 或者说, ECM 算法第 t 次迭代的第 s 个CM 步
骤是获得 H( t+ sö S )使
Q (H( t+ s ö S ) ßH( t) ) ≥Q (HßH( t) ) H∈ ( (2)
完全资料似然函数的条件期望Q (HßH( t) )的极大值点由下列公式确定:
P( t+ 1)j = ∑
n1
i= 1
w
( t)
j i ö n1 (3)
L (Y ßH) - ∑
k
m = 1
Km rm (L) ö Hs = 0 (s = 1, ⋯, S ) (4)
其中, w ( t)j i 是第 t 次迭代后第 i 个观测值归入第 j 个成分分布的后验概率, rm (L)是分布平均数
间的第m 个约束条件, k 是约束条件个数, n1 是样本容量。取 H( t+ 1) = H( t+ S ö S ) , 以此进行下一
轮循环。
3. 2 数量性状主+ 多基因混合遗传分析的 IECM 算法
数量性状主+ 多基因混合遗传分析的 IECM 算法由 E 步骤和 2~ 3 个CM 步骤组成。其
E 步骤与 ECM 算法的 E 步骤是一致的。CM 1 步骤是在固定多基因方差组分 R2j0 ( j = 4, 5, 6)
和环境方差的条件下用迭代方法求分布平均数的条件极大似然估计; CM 2 是在固定环境方差
和CM 1 步骤中获得的分布平均数的条件下用迭代公式求多基因方差组分的条件极大似然估
计; CM 3 步骤是在固定CM 1 和CM 2 中获得的分布平均数和多基因方差组分条件下用迭代公
式求环境方差的条件极大似然估计。若涉及家系世代, CM 1 步骤中分布平均数条件极大似然
估计可按下列步骤进行: ① 若分布平均数间有约束条件, 由约束条件和平均数公式得到的联
立方程组求L agrange 乘数 Ki; ② 由分布平均数公式求其估计值; ③ 由分布平均数估计值得
到一阶遗传参数估计值及其表示的主基因方差组分, 从而改变了家系群体成分分布方差; ④
重复①~ ③步骤直到平均数变化满足预定的精度为止。
为节省篇幅, 本文的 IECM 算法均略去 E 步骤和估计分布平均数的 CM 1 步骤, 其符号
和样本似然函数参见文献[2~ 5 ], 这里只列出多基因和误差两种方差的迭代公式。
3. 3 利用个别分离世代分离分析的 IECM 算法
3. 3. 1 鉴定主基因存在 IECM 算法的迭代公式 F 2 或 F 2∶3世代鉴定主基因存在的样本似
然函数的形式参见文献[2, 3 ]。只是这里的成分分布是按前述的主基因型来确定的。记 x i、n1
和 k 分别是 F 2 或 F 2∶3群体的第 i 个观测值或家系平均数、样本容量和成分分布个数。CM 2 是
在固定成分分布平均数条件下求 R21 的条件极大似然估计, 其迭代公式为:
1076 期 章元明等: 数量性状分离分析中分布参数估计的 IECM 算法
R21 = ∑
k
j= 1
v
2
j∑
n1
i= 1
w j i (x i - Lj ) 2 ∑
k
j = 1
v j∑
n1
i= 1
w j i (5)
其中, F 2∶3群体有 v j = R21 ö R2j , F 2 群体有 v j = 1。
B 1 和B 2 或B 1∶2和B 2∶2世代鉴定主基因存在的样本似然函数为:
f (Y ßH) = ∑
n1
i= 1
∑
k1
j= 1
P1j f (x 1i; L1j , R21j )∑
n2
i= 1
∑
k2
j= 1
P2j f (x 2i; L2j , R22j ) (6)
其中, x 1i和 x 2i、n1 和 n2、k 1 和 k 2、P1j和 P2j、R21j和 R22j分别是B 1 和B 2 或B 1∶2和B 2∶2群体的第 i
个观测值或家系平均数、样本容量、成分分布个数、第 j 个成分分布的后验概率与方差。CM 2
步骤是在固定成分分布平均数条件下求 R211和 R22k2的条件极大似然估计, 其迭代公式分别为:
R211 = ∑
k1
j= 1
v
2
1j∑
n1
i= 1
w 1j i (x 1i - L1j ) 2 ∑
k1
j= 1
v 1j∑
n1
i= 1
w 1j i (7a)
R22k2 = ∑
k2
j = 1
v
2
2j∑
n2
i= 1
w 2j i (x 2i - L2j ) 2 ∑
k2
j= 1
v 2j∑
n2
i= 1
w 2j i (7b)
其中, B 1∶2和B 2∶2群体有 v 1j = R211 ö R21j , v 2j = R22k2 ö R22j; B 1 和B 2 群体有 v 1j = v 2j = 1. 0 。这里的成
分分布也是按主基因型来确定的。将两回交群体合并分析是为了拓展两对主+ 多基因混合遗
传模型。
3. 3. 2 鉴定多基因存在的 IECM 算法 从文献[2 ]可知, 在个别分离世代基础上增加 P 1、
F 1 和 P 2 三个同质群体估计环境方差以鉴定多基因存在。利用亲本、F 1 和 F 2 或 F 2∶3鉴定多基
因存在的符号和样本似然函数等参见文献[ 2 ]。在构造H 0 时, 多基因不存在包括多基因效应
平均数和多基因方差均为 0 两个方面。由在H 0 和H a 下的最大对数似然函数值L 0 和L a 构造
的似然比统计量 K= 2 ( ln L a- ln L 0)~ x 2df 可鉴定多基因是否存在, 其自由度 df 为两种假设下
相差的遗传参数个数。CM 2 步是在固定误差方差 R2e 和分布平均数, 求H a 条件下 F 2 或 F 2∶3群
体成分分布方差中的多基因方差组分 R240的条件极大似然估计, 其迭代公式为:
R240 = ∑
k
j= 1
v
2
j∑
n4
i= 1
w j i (x 4i - L4j ) 2 ∑
k
j= 1
v j∑
n4
i= 1
w j i - A (8)
其中, F 2 群体有A = R2e 和 v j = 1. 0; F 2∶3群体有A = R2e ö n 和 v j = R241 ö R24j。在H 0 条件下, R240不存
在, 即无这一步骤。CM 3 步是在固定 F 2 或 F 2∶3群体成分分布平均数和方差中的多基因方差
组分 R240, 求误差方差 R2e 的条件极大似然估计, 其迭代公式为:
R2e = ∑
3
j = 1
∑
n j
i= 1
(x j i - Lj ) 2 + A ∑
k
j= 1
v
2
j∑
n4
i= 1
w j i (x 4i - L4j ) 2 ∑
3
i= 1
n i + ∑
k
j= 1
v j∑
n4
i= 1
w j i (9)
其中, F 2∶3群体有 v j = (R2e ö n) ö R24j和A = n , F 2 群体有 v j = R2e ö R24j和A = 1. 0。
若利用亲本、F 1、B 1 和B 2 (或B 1∶2和B 2∶2)鉴定多基因存在, x ti和 n t ( t= 1, 2, 3)的含义同
上, 记 x 4i和 x 5i, 以及 n4 和 n5 分别为B 1 和B 2 或B 1∶2和B 2∶2群体的观测值或家系平均数和样
本容量, 则B 1 和B 2 或B 1∶2和B 2∶2群体分别是 k 1 个N (L4j , R24j )和 k 2 个N (L5j , R25j )的混合。由
此, 样本似然函数为:
f (Y ßH) = ∏
n1
i= 1
f (x 1i; L1, R2e )∏
n2
i= 1
f (x 2i; L2, R2e )∏
n3
i= 1
f (x 3i; L3, R2e )
∏
n4
i= 1
∑
k1
j = 1
P4j f (x 4i; L4j , R24j )∏
n5
i= 1
∑
k5
j = 1
P5j f (x 5i; L5j , R25j ) (10)
207 作 物 学 报 26 卷
H 0、H a 和似然比统计量 K可仿上述内容来构造。CM 2 步是在固定环境方差 R2e 和分布平均数,
求H a 条件下B 1 和B 2 或B 1∶2和B 2∶2群体成分分布方差中的多基因方差组分 R240和 R250的条件极
大似然估计, 其迭代公式为:
R2m 0 = ∑
k
m - 3
j= 1
v
2
m - 3, j∑
n
m
i= 1
w m j i (x m i - Lm j ) 2 ∑
k
m - 3
j= 1
vm - 3, j∑
n4
i= 1
w m j i - A (m = 4, 5) (11)
其中, B 1 和B 2 群体有A = R2e 和 v 1j = v 2j = 1. 0; B 1∶2和B 2∶2群体有A = R2e ö n , v 1j = R241 ö R24j , v 2j =
R25k2 ö R
2
5j。在H 0 条件下, 无这一步骤。CM 3 步是在固定B 1 和B 2 或B 1∶2和B 2∶2群体成分分布方
差中的多基因方差组分 (R240和 R250) 和分布平均数的条件下, 求误差方差 R2e 的条件极大似然估
计, 其迭代公式为:
R2e = ∑
3
j = 1
∑
n j
i= 1
(x j i - Lj ) 2 + A ∑
5
t= 4
∑
k t- 3
j= 1
v
2
t- 3, j∑
n t
i= 1
w tj i (x ti - Ltj ) 2
∑
3
i= 1
n i + ∑
5
t= 4
∑
k t- 3
j= 1
v t- 3, j∑
n t
i= 1
w tj i (12)
其中, B 1∶2和B 2∶2群体有 v 1j = [ R2e ö n ] ö R24j , v 2j = [ R2e ö n ] ö R25j和A = n; B 1 和B 2 群体有 v 1j = R2e ö
R24j、v 2j = R2e ö R25j和A = 1. 0。
3. 4 利用联合多个世代群体分离分析的 IECM 算法
P 1、F 1、P 2、B 1、B 2 和 F 2 6 世代分离分析。其符号、基本假定和似然函数参见文献[ 5 ]。
CM 2 步是在固定环境方差 R2e 和分布平均数, 求B 1, B 2 和 F 2 群体成分分布方差的多基因方差
组分 (R240、R250和 R260)的条件极大似然估计, 其迭代公式为:
R2m 0 = S S m ö nm - R2e S S m = ∑
k
m - 3
j = 1
∑
n
m
i= 1
w m j i (x m i - Lm j ) 2 (m = 4, 5, 6) (13)
CM 3 步是在固定多基因方差组分 (R240、R250和 R260) 和分布平均数, 求误差方差 R2e 的条件极大似
然估计, 其迭代公式为:
R2e = ∑
3
t= 1
∑
n t
i= 1
(x ti - Lt) 2 + ∑
6
t= 4
v
2
tS S t ∑
3
t= 1
n t + ∑
6
t= 4
v tn t (14)
其中, v t= R2e ö (R2t0+ R2e ) , t= 4, 5, 6。
P 1、F 1、P 2、F 2 和 F 2∶3 5 世代分离分析。其符号、基本假定和似然函数参见文献 [ 4 ]。
CM 2 步骤是固定环境方差 R2e 和分布平均数, 求 F 2 和 F 2∶3群体成分分布方差中的多基因方差
组分 (R240和 R250)的条件极大似然估计, 其迭代公式分别为:
R240 = S S 4 ö n4 - R2e S S 4 = ∑
k1
j= 1
∑
n4
i= 1
w 4j i (x 4i - L4j ) 2 (15a)
R250 = ∑
k2
j = 1
v
2
j∑
n5
i= 1
w 5j i (x 5i - L5j ) 2 ∑
k2
j = 1
v j∑
n5
i= 1
w 5j i - R2e ö n (15b)
其中, v j = R251 ö R25j , j = 1, ⋯, k 2。CM 3 步是在固定多基因方差组分 (R240和 R250) 和分布平均数的
条件下, 求环境方差 R2e 的条件极大似然估计, 其迭代公式为:
R2e = ∑
3
t= 1
∑
n t
i= 1
(x ti - Lt) 2 + v 24S S 4 + n∑
k2
j= 1
v
2
5j∑
n5
i= 1
w 5j i (x 5i - L5j ) 2
∑
3
t= 1
n t + v 4n4 + ∑
k2
j = 1
v 5j∑
n5
i= 1
w 5j i (16)
3076 期 章元明等: 数量性状分离分析中分布参数估计的 IECM 算法
其中, v 4= R2e ö R241, v 5j = [ R2e ö n ] ö R25j , S S 4 见 (15a)式。
P 1、P 2、F 1、B 1∶2、B 2∶2和 F 2∶3 6 个家系世代分离分析。关于亲本和 F 1、B 1∶2和B 2∶2群体的
符号与§3. 3. 2 一致, 亲本和 F 1 数量性状平均数 x ti~N (Lt, Re ö n) , t= 1, 2, 3。记 x 6i、n6 和
k 3 分别是 F 2∶3群体的第 i 个家系平均数、样本容量和成分分布数。由§1 和§2 可知, B 1∶2、
B 2∶2和 F 2∶3群体分别为 k 1 个N (L4j , R24j )、k 2 个N (L5j , R25j )和 k 3 个N (L6j , R26j )的混合, 则样本
似然函数为:
f (Y ßH) = ∏
n1
i= 1
f (x 1i; L1, R2e ö n)∏
n2
i= 1
f (x 2i; L2, R2e ö n)∏
n3
i= 1
f (x 3i; L3, R2e ö n) ı
∏
n4
i= 1
∑
k1
j= 1
P4j f (x 4i; L4j , R24j ) ∏
n5
i= 1
∑
k2
j= 1
P5j f (x 5i; L5j , R25j )∏
n6
i= 1
∑
k3
j= 1
P6j f (x 6i; L6j , R26j ) (17)
CM 2 步是在固定环境方差 R2e 和分布平均数条件下, 求B 1∶2、B 2∶2和 F 2∶3群体成分分布方差中
的多基因方差组分 (R240、R250和 R260)的条件极大似然估计, 其迭代公式为:
R2m 0 = ∑
k
m - 3
j= 1
v
2
m j∑
n
m
i= 1
w m j i (x m i - Lm j ) 2 ∑
k
m - 3
j = 1
vm j∑
n
m
i= 1
w m j i - R2e ö n (m = 4, 5, 6) (18)
其中, v 4j = R241 ö R24j , v 5j = R25k2 ö R25j , v 6j = R261 ö R26j。CM 3 步是在固定多基因方差组分 (R240、R250和 R260)
和分布平均数的条件下, 求环境方差 R2e 条件极大似然估计, 其迭代公式为:
R2e = ∑
3
t= 1
∑
n t
i= 1
(x ti - Lt) 2 + ∑
6
t= 4
∑
k t- 3
j= 1
v
2
tj∑
n t
i= 1
w tj i (x ti - Ltj ) 2 ∑
3
t= 1
n t + ∑
6
t= 4
∑
k t- 3
j= 1
v
2
tj∑
n t
i= 1
w tj i (19)
其中, v tj = [R2e ö n ] ö R2tj , t= 4, 5, 6。
以上方法看起来繁杂, 但使用计算机也仅一举手之劳。本文作者应用 T u rbo C+ + 语言编
制了以上各种情况的全套 IECM 算法软件 F 2. EXE (利用 F 2 鉴定主基因的存在) 和 F 2P. EXE
(利用亲本、F 1 和 F 2 鉴定多基因的存在)、F 3. EXE (利用 F 2∶3鉴定主基因的存在)和 F 3P. EXE
(利用亲本、F 1 和 F 2∶3鉴定多基因的存在)、B. EXE (利用B 1 和B 2 鉴定主基因的存在) 和BP.
EXE (利用亲本、F 1、B 1 和B 2 鉴定多基因的存在)、FB. EXE (利用B 1∶2和B 2∶2鉴定主基因的
存在)和 FBP. EXE (利用亲本、F 1、B 1∶2和B 2∶2鉴定多基因的存在)、S IN. EXE (利用亲本、
F 1、B 1、B 2 和 F 2 的联合分析)、M IX. EXE (利用亲本、F 1、F 2 和 F 2∶3的联合分析) 和 FAM I.
EXE (利用亲本和 F 1、B 1∶2、B 2∶2和 F 2∶3的联合分析)。欢迎读者来函联系。
4 应用实例
本文以南京农业大学朱立宏教授提供的南京 6 号×广丛杂交组合 6 个基本世代 (P 1、P 2、
F 1、B 1、B 2 和 F 2) 株高资料为例来说明 IECM 算法。用 EM 和 IECM 算法进行分布参数估计
的结果见表 1。从表 1 可知: IECM 算法的结果比 EM 算法的结果好, 更容易收敛, 极大似然
函数值更大。由此可通过A IC 准则进行模型选择和用适合性检验进行模型检验 (另文报道)。
5 讨论
本文迭代公式中分离群体成分分布数是一般的, 主要是为满足A、B、C、D 和 E 类共 24
种遗传模型的缘故。在多世代联合分离分析的A 和B 两类模型中, 没有多基因方差组分的估
计, 即CM 步骤只分两步进行。迄今为止, 利用个别分离世代鉴定主基因存在一般是先确定
成分分布数目然后鉴定主基因是否存在。确定成分分布数的方法可分为图形方法和统计检验
407 作 物 学 报 26 卷
表 1 用 EM 算法和 IECM 算法计算的 24 种遗传模型
的极大对数似然函数值
Table 1 The max imum log- l ikel ihood function values (ML FV)
of 24 genetic modelsca lculated through the EM and
iterated ECM ( IECM ) algor ithm
模型
极大对数似然值M L FV
EM IECM
模型
极大对数似然值M L FV
EM IECM
A 21 3596. 38 3257. 94 D 3487. 85 3185. 95
A 22 3984. 39 3460. 48 D 21 3512. 90 3185. 30
A 23 3618. 97 3300. 69 D 22 3801. 36 3185. 30
A 24 4122. 60 3631. 70 D 23 3514. 13 3264. 98
B21 3521. 80 3224. 81 D 24 4160. 12 3422. 24
B22 3538. 69 3225. 94 E 3185. 93
B23 3983. 97 3381. 16 E21 3183. 80
B24 4093. 22 3422. 23 E22 3172. 06
B25 3618. 85 3444. 06 E23 3226. 75
B26 - 3389. 06 E24 3254. 51
C 3782. 60 3287. 47 E25 3195. 23
C21 3804. 39 3290. 65 E26 3337. 26
方法两大类[ 2, 3, 17 ]。后者有采用
A IC 准则[ 2, 3, 17 ]、B IC 准则[ 17 ]和自
助 (Boo tst rap )等方法。这些方法常
会导致比如在利用 F 2 世代进行主
基因存在的鉴定中不能有效地区分
1 个位点的 1∶2∶1 的 3 个正态分
布的混合和 2 个位点的 9∶6∶1 的
3 个正态分布的混合; 也不能有效
地区分 3∶1 的完全显性的 2 个正
态分布的混合和 1∶3 的负向完全
显性的 2 个正态分布的混合。从参
数估计角度来看, 因不同模型的分
布平均数间约束条件不同, 致使分
布平均数的迭代公式不同, 因此其
参数估计结果是不一致的。本研究
的利用个别分离世代鉴定主基因存
在是建立在一定的遗传模型基础上
的, 根据特定的遗传模型来确定成分分布数的, 具有特定的遗传背景, 其分布参数有特定的
遗传学意义, 这把分布参数的估计和遗传参数的估计有机地结合在一起, 具体表现在一阶分
布参数估计中考虑了通过一阶遗传参数得到的约束条件和鉴定多基因存在中剖分了分离世代
成分分布方差的数量遗传学构成, 有效地克服了传统分析方法的缺陷[ 2, 3, 18 ]。鉴定多基因存
在中的统计假设不仅考虑了多基因存在与否的分布方差构成, 还应考虑多基因是否存在的分
布平均数。这就是说, 在 H 0 中不仅令多基因方差组分为 0, 还应令多基因效应平均数为 0,
这才比较完整地表示出多基因不存在; 在H a 中两者均应存在。鉴定多基因存在的似然比检
验统计量 K的自由度不是固定为 1 而是在H 0 和H a 中相差的遗传参数个数。
自 EM 算法[ 8 ]提出后, 虽对 EM 算法的改进算法较多, 如用相对简单的 S 个CM 步骤代
替计算复杂的 EM 算法M 步骤的 ECM 算法[ 12 ] , 在CM 步骤中用实际似然函数代替 ECM 算
法中的期望完全资料似然函数的 ECM E 算法[ 10 ] , 用 EM 算法获得近似方差2协方差矩阵的
SEM 算法[ 11 ] , 但在主+ 多基因混合遗传分析的参数估计中一般是采用 EM 算法[ 2~ 7, 9, 16 ]。然
而, 在 2 对主+ 多基因混合遗传分析中, 采用 EM 算法进行参数估计显然是不可靠的, 特别
是有家系世代时。这是因为它忽略的项数较多。若采用 ECM 算法进行参数估计, 也因多基
因和环境两种方差估计时会出现高次方程, 用 Gau ss2N ew ton 法解高次方程时其解可能不唯
一致使迭代可能不收敛。本文以 EM 和 ECM 算法为基础从数量遗传学观点出发将成分分布
方差剖分为主基因、多基因和环境三种方差组分, 推导出CM i 步骤的一般迭代公式。经大量
计算表明, 只要方差组分不为负数, CM i 步骤中的迭代是收敛的。此外, 还有效地降低了参
数估计维数, 节约了计算时间, 在环境方差估计中利用了所有群体的所有分布, 克服了以前
方法在估计误差方差时只考虑了亲本和 F 1 群体分布方差的缺陷。
在参数估计前选择分布参数初值对于 2 对主+ 多基因混合遗传分析是十分重要的。Bake
(1994) 将 EM 算法与N ew ton2R aph son 迭代结合克服初值选择问题[ 19 ] , L a ird (1978) 建议用
5076 期 章元明等: 数量性状分离分析中分布参数估计的 IECM 算法
网点搜索[ 17 ] , M cL ach lan 等 (1988)建议用聚类方法先分组再估计参数作为初值[ 17 ] , 本文就是
采用这一方法, 也有用多组初值的结果进行比较选择的。关于度量参数估计值误差大小的参
数标准误, 有 F isher 信息阵和观察信息阵两种方法[ 14, 17 ]。L it t le 等 (1987) 认为前者不精确,
建议不用; 后者的计算相当复杂。因此, 本文未给出其标准误。在参数估计中停止迭代的准
则有迭代间对数似然函数值变化很小 (如小于 10- 4)和方向导数准则等[ 2, 3, 18 ]。T it tering ton 等
(1985)认为: 停止早迟将会影响参数估计值, 采用前一方法可能会过早停止迭代[ 18 ]。从本文
的结果来看, IECM 算法比 EM 算法更易收敛。
在本文水稻株高性状遗传分析中, 通过 IECM 算法从前四类模型 (表 1) 中只能发现表现
为完全显性的 1 个主基因位点的存在, 与拓展到两对主+ 多基因混合遗传分析结果中的最大
效应主基因是一致的, 此外还有一个表现为次大效应的主基因。进一步的比较分析还发现利
用个别世代和多世代联合的分析结果是一致的, 并且后者优于前者, 利用家系群体的结果优
于利用单株群体的结果, 利用自交群体的结果优于利用回交群体的结果。IECM 算法保证了
参数估计时良好的收敛性, 因而更充分地揭示数据内含的信息。
参 考 文 献
1 盖钧镒, 管荣展, 王建康. 世界科技研究与发展, 1999, 21 (1) : 34~ 40
2 王建康, 盖钧镒. 遗传学报, 1997, 24 (5) : 432~ 440
3 盖钧镒, 王建康. 作物学报, 1998, 24 (5) : 402~ 409
4 王建康, 盖钧镒. 作物学报, 1998, 24 (6) : 651~ 659
5 Gai J Y, J K W ang. T heor. A pp l. Genet. , 1998, 97: 1162~ 1168
6 王建康, 盖钧镒. 生物数学学报, 1995, 10 (4) : 87~ 92
7 王建康, 盖钧镒. 生物数学学报, 1997, 12 (5) : 540~ 548
8 D emp ster A P, N M L aird, D B Rubin. J . R. S ta tist. S oc. B . , 1977, 39: 1~ 38
9 Kao C H , Z B Zeng. B iom etrics, 1997, 53: 653~ 665
10 L iu C, D B Rubin. B iom etrika, 1994, 81 (4) : 633~ 648
11 M eng X L , D B Rubin. J . of the A m erican S ta tist. A ssocia tion, 1991, 86: 899~ 909
12 M eng X L , D B Rubin. B iom etrika, 1993, 80 (2) : 267~ 278
13 Rai S N , E M atthew s. B iom etrics, 1993, 49: 587~ 591
14 Shouk riM M , G J M cL ach lan. B iom etrics, 1994, 50: 128~ 139
15 U im ari P, I Hoeschele. Genetics, 1997, 146: 735~ 743
16 Jansen R C. T heor. A pp l. Genet. , 1992, 85: 252~ 260
17 L eroux B G, M L Puterm an. B iom etrics, 1992, 48: 545~ 558
18 Bohning D , P Sch lattm an. B iom etrics, 1992, 48: 283~ 303
19 M o lenbergh s G, E Goetghebeur. J R S ta tist. S oc. B . , 1997, 59 (2) : 401~ 414
607 作 物 学 报 26 卷