全 文 :作物学报 ACTA AGRONOMICA SINICA 2009, 35(11): 1981−1989 http://www.chinacrops.org/zwxb/
ISSN 0496-3490; CODEN TSHPA9 E-mail: xbzw@chinajournal.net.cn
本研究由国家自然科学基金项目(30571072)资助。
第一作者联系方式: E-mail: xiyuanhu@yahoo.com.cn; Tel: 029-87081390
Received(收稿日期): 2009-02-23; Accepted(接受日期): 2009-06-25.
DOI: 10.3724/SP.J.1006.2009.01981
基于协方差阵结构优选的作物品种区域试验分析
胡希远 尤海磊 任长宏 吴 冬 李建平
西北农林科技大学农学院, 陕西杨凌 712100
摘 要: 论述了线性混合模型方差协方差结构与作物品种区域试验分析模型的对应关系, 以我国 2005—2006年东北
华北玉米 8 组区域试验资料为例, 按照线性混合模型分析原理及模型拟合信息量准则与似然比测验, 对区域试验品
种方差协方差的结构特性及不同方差协方差结构模型在品种效应估计与评价的差异状况进行了探讨。结果表明, 在
分析的所有试验中, 环境间品种效应方差协方差均不符合方差分析模型假设的同质性结构, 而是呈现为各种异质性
结构; 产量效应测验差异显著的品种对数目在方差分析模型与最佳方差协方差结构线性混合模型间的一致率平均为
86%, 品种产量效应排序在两种模型间也存在明显不同, 品种产量效应估计的平均误差在最佳方差协方差结构线性
混合模型小于在方差分析模型。
关键词: 作物品种; 区域试验; 线性混合模型; 方差协方差
Analysis of Crop Variety Regional Trials Based on Selection of Covariance
Structures
HU Xi-Yuan, YOU Hai-Lei, REN Chang-Hong, WU Dong, and LI Jian-Ping
College of Agronomy, Northwest A&F University, Yangling 712100, China
Abstract: The method mainly used for analyzing crop variety regional trials is based on analysis of variance (ANOVA), which
requires a homogenous variance-covariance of data. Now, other models different from the ANOVA model are available. However,
the problems that how the models should be assessed and that which model is more suitable for given trial data are not solved and
hence restrict the applicability of the models in practices. This paper tried to solve these problems on the basis of liner mixed
models. Relations between various variance-covariance structures of linear mixed model and models available for analyzing crop
regional trials were discussed. Then on the basis of analyses of the corn regional trials in northeast and north China, using the
information criterion and likelihood-ratio-test, the characteristics of variance-covariance structures of regional trial data and the
performance difference between the ANOVA model and the linear mixed model with optimal variance-covariance structure were
assessed. The results showed that the variance-covariance of variety effect over environments was not homogeneous as defined in
the ANOVA model, but heterogeneous in all the considered trials. The ratio of the same variety contrast with significant difference
between the ANOVA model and the optimal linear mixed model averagely reached 86%. Also, there was obvious difference in the
yield ranking of varieties between the two models. The error of variety effect estimation in the optimal linear mixed model was
smaller than that in the ANOVA model.
Keywords: Crop variety; Regional trial; Linear mixed model; Covariance
作物品种区域试验(简称区试)是对参试品种在
多环境下的试验, 是作物品种审定和为实际生产推
荐优良品种的重要依据。品种性状评价的准确性直
接影响着优良品种推荐工作的可靠性和实际农业生
产的成效。由于多环境下品种与环境互作效应形成
机理的复杂性, 关于区试数据的分析目前还没有形
成一个普遍认可的最佳模型或方法。长期以来, 我
国区试数据的分析主要采用传统的方差分析法
(ANOVA)。该方法在包括作物区试在内的各种试验
的统计分析中发挥了重要作用, 而且在有限的计算
条件下也是实际能够实现的最佳分析方法。然而 ,
方差分析模型具有试验误差同质和不相关及随机效
应方差协方差同质的严格假定条件[1-3], 实际试验不
一定能完全满足这些假设条件, 所以方差分析在许
1982 作 物 学 报 第 35卷
多试验数据分析的应用仅仅是一种一定程度上的近
似方法[4]。在现阶段, 由于计算技术的进步, 为更多
或更完善试验分析方法的常规化应用提供了条件 ,
方差分析法的应用已不再是唯一的最佳选择。例如,
线性模型加权最小二乘法和 Wricke-Shukla 模型可
处理方差不同质的问题 , Eberhart-Russell 模型、
AMMI (additive main effects and multiplicative
interaction)模型和 Finlay-Wilkinson 模型等, 既可处
理方差不同质又可处理协方差不同质的问题[5-10]。不
少文献中的研究结果也证明, 对一些区试, 这些模
型或方法比方差分析法在品种效应估计准确性等方
面有一定的改进[11-13]。但这些模型处理品种-环境互
作效应的方式和假设条件(或者方差协方差的组成)
不尽相同, 各模型中待估参数多少也不同, 不是每
个模型对任何区试资料分析效果都会有改进。对任
一试验资料究竟采用哪种模型分析效果才更好不易
确定, 限制了这些模型在实际区试分析中的广泛应
用。因此, 各种模型在实际区试分析效果的评价与
适用模型的选择是目前区试分析中亟待解决的重要
问题之一。
如今, 线性混合模型理论已相当完善, 统计分
析软件的发展也很好地解决了线性混合模型分析的
复杂计算问题, 使线性混合模型在各研究领域的应
用日益广泛。该模型不受方差分析模型和其他区试
分析模型假定条件的限制, 可借助多种方差协方差
结构来反映各类试验数据或分析模型的特征。而且,
对一特定的试验, 可以通过线性混合模型对试验数
据拟合优劣的指标选用最佳的方差协方差结构, 实
现从多种模型中选取最佳的模型进行分析, 从而提
高分析结论的可靠性。本文拟在介绍线性混合模型
分析原理的基础上, 指出线性混合模型几种主要方
差协方差结构与各种区试分析模型的关系, 然后通
过对实际区试数据的模型拟合与分析, 探讨基于方
差协方差结构选择的线性混合模型在区试分析中的
作用和效果, 以期为区试分析模型选择问题的解决
提供一可行方法和对传统区试分析方法的改进有所
启示。
1 材料与方法
1.1 试验设计
线性混合模型的矩阵形式[14-15]为:
y = Xβ+Zu+e, (1)
式中, y为试验观测值向量, X为固定效应设计矩阵,
β为固定效应向量, Z为随机效应设计矩阵, u为随机
效应向量, 服从均值向量为 0、方差协方差矩阵为 G
的正态分布, 通常表示为 u~N(0, Gp×p), p 为随机效
应向量 u的维数; e为剩余误差向量, 服从均值为 0、
方差协方差矩阵为 R的正态分布, 表示为 e~N(0, Rn
×n), n为观测值向量 y的维数。同时, 假定 u与 e间
不相关, 即 Cov(u, e)=0。这时 y 的方差协方差矩阵
为 Var(y)=ZGZ′ =R=V, y的期望值为 E(y)=Xβ。
与传统 ANOVA 模型相比较, 上述线性混合模
型的根本区别在于它没有对 u 和 e 方差协方差结构
的严格限制, 即 G 和 R 可根据所分析数据的特征呈
现各种方差协方差结构形式。若 R为任意方阵形式,
表示剩余误差方差不同质, 且剩余误差协方差(或相
关系数)不相等; 若 R是主对角线元素不相等而非主
对角线元素相等的方阵, 则表示剩余误差方差不同
质, 但剩余误差协方差相等; 若 R 是主对角线元素
相等非主对角线元素为零的对角方阵, 则表示剩余
误差方差同质协方差为零(即剩余误差不相关)。类似
地, G的结构则反映了随机效应 u的方差协方差特性。
依据线性混合模型分析原理, 模型(1)中固定效
应β的估计为 βˆ =(X′V−1X)−X′V−1y, 其估计误差方
差为 var βˆ =(X′V−1X)−。其中,“−”表示广义逆阵。
如果方差协方差参数值已知, 效应估计为最佳线性
无偏估计(BLUE)。
实际应用中, 方差协方差参数未知, 需要根据
样本数据, 利用限制性极大似然法(REML)等方法进
行估计[16-18], 然后利用估计的 Vˆ 代替 V 进行模型效
应估计和有关统计测验。模型中固定效应线性对比
H0: ˆLβ =0由下面的统计量进行测验[15]。
对两个固定效应的线性对比,
t =
ˆ
ˆvar( )
L
L L
β
β ′
~(tdf)
对多个固定效应的线性对比,
F =
1ˆ ˆ ˆ[ var( ) ]
( )
L L L L
rank L
β β β−′ ′ ′ ~F[rank(L), df]
其中, L 是描述固定效应线性对比的行向量(两个固
定效应的线性对比)或矩阵(多个固定效应的线性对
比); rank(L)表示 L 的秩; df 为测验统计量的自由度
(第二自由度), 通常用 Satterthwaite法计算。
1.2 区试分析模型的方差协方差结构
区试中品种的效应一般可表示为μ i j =μ+α i +
βj+(αβ)ij。其中, μij为第 i个品种在第 j个环境中的
效应, μ为总体平均值, αi为第 i个品种的主效应, βj
第 11期 胡希远等: 基于协方差阵结构优选的作物品种区域试验分析 1983
为第 j个环境的主效应, (αβ)ij为第 i个品种与第 j个
环境的品种-环境互作效应。通常认为αi为固定效应,
βj与 (αβ)ij都为随机效应。假设βj与 (αβ)ij的方差分
别为 2βσ 和 2αβσ 。关于μij的方差协方差在线性混合模
型分析中至少可采用如下 6种结构。
(1) CS结构,
Var
2 2 2 2
1
2 2 2 2
2
2 2
2 2 2 2
j
j
ij
β αβ β β
β β αβ β
β β
β β β αβ
σ σ σ σμ
μ σ σ σ σ
σ σ
μ σ σ σ σ
⎛ ⎞+⎛ ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟+⎜ ⎟ = ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟+⎝ ⎠ ⎝ ⎠
#
# # %
,
该结构表示所有品种的方差相等, Var(μij) = 2βσ + 2αβσ ,
所有品种间的协方差也相等, Cov(μij, μi′j) = 2βσ 。
(2) UN(1)结构,
Var
2 2 2 2
(1)1
2 2 2 2 2
2 (2)
2
2 2 2 2
( )
j
j
ij i
β αβ β β
β β αβ β β
β
β β β αβ
σ σ σ σμ
μ σ σ σ σ σ
σ
μ σ σ σ σ
⎛ ⎞+⎛ ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟+⎜ ⎟ ⎜ ⎟=⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎝ ⎠ +⎝ ⎠
# # % #
,
其中 2 ( )iαβσ 为品种 i与环境互作效应的方差。此结构
与 CS 结构的不同之处在于各品种的方差不同, 即
Var(μij) = 2 2 ( )iβ αβσ σ+ 。
(3) FA1(1)结构,
Var
2 2
1 1 1 2 1
2 2
2 2 1 2 2
2 2
1 2
j i
j i
ij i i i
δ
δ
δ
μ λ σ λ λ λ λ
μ λ λ λ σ λ λ
μ λ λ λ λ λ σ
⎛ ⎞⎛ ⎞ +⎜ ⎟⎜ ⎟ ⎜ ⎟+⎜ ⎟ = ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟⎜ ⎟ ⎜ ⎟+⎝ ⎠ ⎝ ⎠
# # # % #
,
其中λi 为第 i 个品种对环境因子变量的回归系数,
2δσ 为品种-环境互作效应中未被回归解释部分的方
差。此结构表示各品种的方差和各对品种间的协方
差都不同 , 方差为 Var(μij) = 2 2i δλ σ+ , 协方差为
Cov(μij, μi′j) =λiλi′。
(4) AMMI结构,
Var
2 2
1 1 1 2 1
2 2
2 2 1 2 2
2 2
1 2
j i
j i
ij i i i
δ
δ
δ
μ λ σ λ λ λ λ
μ λ λ λ σ λ λ
μ λ λ λ λ λ σ
⎛ ⎞⎛ ⎞ +⎜ ⎟⎜ ⎟ ⎜ ⎟+⎜ ⎟ = ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟⎜ ⎟ ⎜ ⎟+⎝ ⎠ ⎝ ⎠
# # # % #
+ 2i i βσ×J ,
其中 Ji×i为所有元素等于 1的 i阶方阵, 其他符号的
意义同前述。与 FA1(1)结构类似, 此结构也表示各
品种的方差和各对品种间的协方差都不同, 但方差
协方差的构成与 FA1(1)结构的不同。方差协方差都
增加了一个环境效应方差组分 2βσ , 即 Var(μ i j )=
2 2 2
i δ βλ σ σ+ + , Cov(μij μi′j) = 2i i βλ λ σ+ 。
(5) FA(1)结构,
2 2
1 (1) 1 2 11
2 2
2 2 1 2 (2) 2
2 2
1 2 ( )
2
Var
,
ij
j i
ij i i i i
i i
δ
δ
δ
β
λ σ λ λ λ λμ
μ λ λ λ σ λ λ
μ λ λ λ λ λ σ
σ×
⎛ ⎞+⎛ ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟+⎜ ⎟ = ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟+⎝ ⎠ ⎝ ⎠
+
# # # % #
J
其中 2( )iδσ 为第 i 个品种与环境互作效应未被回归解
释部分的方差, 其他符号的意义与 AMMI 结构中的
相同。与 AMMI结构不同的是, 此结构中品种-环境
互作效应未被回归解释部分的方差在各品种间不同,
分别为 2( )iδσ , 品种方差为 Var(μij)= 2 2 2( )i iδ βλ σ σ+ + 。
(6) UN结构,
(11) (12) (1I)1
(21) (22) (1I)2
(1I)
I (I1) (I2) (II)
Var
j
j
j
σ σ σμ
σ σ σμ
σ
μ σ σ σ
⎛ ⎞⎛ ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ = ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠
# # %#
,
该结构中 Var(μij) = 2( )iσ =σii, Cov(μij, μi′j) = σ(ii′), 不
仅各品种的方差和各对品种间的协方差都不同, 而
且对品种的方差协方差没有任何约束条件。
显然, 前 5 种结构都是一定约束(或假设) 条件
下 UN 结构的特殊形式[19], 具有前 5 种方差协方差
结构的线性混合模型依次等价于区试分析 ANOVA
模型、Wricke-Shukla模型[7]、Finlay-Wilkinson模型[8]、
AMMI模型[9, 20]和 Eberhart-Russell模型[9]。
1.3 方差协方差结构评判与测验
如前所述, 线性混合模型具有多种方差协方差
结构可供使用, 究竟哪种结构能更准确地反映试验
数据特征, 即如何进行模型评价与选择的问题, 可
采用模型拟合信息量准则予以解决。常用的 AIC 信
息量准则(Akaike’s information criterion)为[21]:
AIC = − 2 lnL + 2q,
lnL 为模型拟合极大似然值的自然对数值, q 为模型
中待估计方差协方差参数的数目。该式中右边第一
项可解释为模型对试验数据拟合优良度的一个度量,
第二项可解释为对增加模型参数个数的一种平衡。
这样一来, 当有若干个模型可供应用时, 数据拟合
效果好, 而又尽可能节省参数数目的模型, 即 AIC
值最小的模型为最佳模型。
所选方差协方差结构线性混合模型与传统方差
分析模型有无显著差异的问题 , 采用似然比测验
1984 作 物 学 报 第 35卷
(Likelihood-Ratio Test, 简称 LRT)来解决。LRT 统
计量为:
LRT = − 2 (lnL0 − lnLs) (2)
其中, lnL0和 lnLs分别为方差分析模型和所选方差
协方差结构线性混合模型对试验数据拟合的极大似
然值的自然对数。在零假设条件下(模型差异不显著),
LRT 服从χ2分布, 其自由度 df 为所选方差协方差结
构线性混合模型与方差分析模型方差协方差参数数
目之差。
1.4 数据资料与分析
所用数据为中国玉米新品种动态国家级玉米品
种区试报告中 2005 和 2006 年东北华北地区春玉米
组区域试验小区产量(kg)资料。该区每年有 4个区试
组别 (不同组别参试的品种不同 ), 每组试验品种
15~17 个, 所有组别试验设计为随机完全区组设计,
重复 3次, 小区面积 12 m2。所有试验以农大 108为
标准对照品种。2005年每组有 22个试点, 资料为平
衡数据, 2006年每组约有 23个试点, 该年由于有些
品种在一些试点资料不全, 其资料为非平衡数据。
每组区试资料采用前述 6 种方差协方差结构线性混
合模型(其中 CS 结构线性混合模型等价于 ANOVA
模型)进行模型拟合分析, 并分别利用 ANOVA 模型
和所选最佳方差协方差结构线性混合模型进行品种
效应估计和效应差异显著性测验。鉴于每一试验的
数据相当多, 采用所谓两步分析法, 即先对各试点
分别进行方差分析模型的分析, 估计出各试点品种
效应μij 的值和试验误差, 然后对所有试点μij 和误差
的估计值分别进行 ANOVA 模型和方差协方差结构
线性混合模型的联合分析。分析所用的 SAS程序及
有关说明见本文附录。
2 结果与分析
2.1 不同方差协方差结构模型对试验数据拟合
的效果
从表 1中 AIC最小值 (表 1中下划线的数值)对
应的模型可知, 2005年第一组试验的最佳模型为UN
结构模型, 而第二、三和四组的最佳结构模型均为
AMMI 模型; 2006 年第一、二、三和四组的最佳结
构模型均为 FA ( 1 )模型。最佳结构模型与相应
ANOVA 模型比较的 LRT 值[计算方法见式(2)]在
2005 年第一、二、三和四组试验分别为 317.6
(1077.6–760.0)、70.1 (1001.0–930.9)、44.6 (1003.8–
959.2)和 74.4 (1077.2–1002.8), 对应的自由度 df 依
次为 151 (154–3)、16 (18–2)、17 (19–2)和 16 (18–2)。
这些 LRT 值分别大于对应的χ20 .01 (151)= 194.3、
χ20.01(16)= 32.0和χ20.01(17) = 33.4。该χ2测验结果说明,
相应的方差协方差异质性结构线性混合模型极显著
优于方差协方差同质性结构 ANOVA模型。2006年
最佳方差协方差结构模型与相应 ANOVA 模型比较
的 LRT 值在第一、二、三和四组试验分别为 100.7
(1098.6–997.9)、93.9 (1026.8–932.9)、65.5 (968.2–
902.7)和 94.2 (1059.9–965.7), 对应的自由度 df依次
为 32 (35–3)、32(35–3)、26 (29–3)和 31 (34–3)。该
LRT值大于对应的χ20.01(32) = 53.5、χ20.01(26) = 45.6和
χ20.01(31) = 52.2。2006年 4组试验χ2测验的结果同样
说明, 异质性方差协方差结构线性混合模型极显著
优于同质性方差协方差结构 ANOVA 模型, 只是该
年试验方差协方差异质的结构特性与 2005 年试验
的有所不同。2005 年主要为 AMMI 结构模型, 而
2006 年 4 组试验全部为 FA(1)结构模型。这些结果
表明, ANOVA模型中品种方差协方差同质假设不符
合实际试验数据特性的问题和各试验中品种方差协
方差异质特性不尽相同的现象具有一定的普遍性。
表 1 中所示最佳方差协方差结构模型在不同年份间
比同一年份不同组别间变化大, 说明年份较试验组
别对方差协方差特性有着较大的影响。根据所分析
试验数据的特点这一结果是不难理解的, 因为试验
不同组别只是品种不同但试验的环境相同, 而不同
年份间则品种和环境都不同。UN和 UN(1)方差协方
差结构模型对 2006 年几组试验数据拟合未达到收
敛, 说明该方差协方差结构不适用于这几组试验数
据特性的描述。
2.2 方差协方差结构对品种效应估计与测验的
影响
利用 ANOVA 模型和所选最佳方差协方差结构
线性混合模型对各组试验分析, 得到品种效应两两
比较 t-测验达到显著(α = 0.05)和极显著(α = 0.01)的
品种对数目及其在两种模型间比较的一致性状况汇
总于表 2。可以看出, 无论在显著水平还是极显著水
平 , t-测验得到效应差异显著的品种对数目在
ANOVA 模型和最佳方差协方差结构线性混合模型
间都有一定的差异。两种模型测验结果的一致率[两
模型中都显著(极显著)的品种对数目占总的显著(极
显著)品种对数目的百分率]在 α = 0.05 水平时为
77.8%~92.6%, 在 α = 0.01水平时为 70.6%~93.0%。
若以所有试验的平均值来看, 两种模型测验结果相
同率在两个显著性水平上均约为 86%。这也就是说,
被一种模型测验显著或极显著的品种对中平均有
14%在采用另一模型时测验不出显著或极显著。
除了试验品种间效应差异显著性测验外, 试验
第 11期 胡希远等: 基于协方差阵结构优选的作物品种区域试验分析 1985
表 1 方差分析模型(ANOVA)和各类方差协方差结构(UN(1)、FA1(1)、AMMI、FA(1)和 UN)线性混合模型对 2005和 2006年
各组区试数据拟合的有关统计量
Table 1 Related fitting statistics of ANOVA and the liner mixed model with various variance-covariance structures [UN(1),
FA1(1), AMMI, FA(1), and UN] for the trials in 2005 and 2006
2005 2006 组别
Group
模型
Model −2 lnL AIC q −2 lnL AIC q
一组 Group 1 ANOVA 1077.6 1083.6 3 1098.6 1104.6 3
UN(1) 1049.1 1087.1 19 — — —
FA1(1) 1054.2 1092.2 19 1068.2 1106.2 19
AMMI 1042.0 1082.0 20 1045.2 1085.2 20
FA(1) 1016.2 1088.2 36 997.9 1067.9 35
UN 760.0 1068.0 154 — — —
二组 Group 2 ANOVA 1001.0 1005.0 2 1026.8 1032.8 3
UN(1) 955.6 989.6 17 975.5 1013.5 19
FA1(1) 972.4 1006.4 17 989.4 1027.4 19
AMMI 930.9 966.9 18 984.5 1024.5 20
FA(1) 923.8 983.8 30 932.9 1002.9 35
UN 782.9 1024.9 121 — — —
三组 Group 3 ANOVA 1003.8 1007.8 2 968.2 974.2 3
UN(1) 977.1 1013.1 18 — — —
FA1(1) 969.7 1005.7 18 949.8 985.8 18
AMMI 959.2 997.2 19 926.5 962.5 18
FA(1) 944.1 1012.1 34 902.7 960.7 29
UN 763.0 1037.0 137 — — —
四组 Group 4 ANOVA 1077.2 1081.2 2 1059.9 1065.9 3
UN(1) 1027.2 1063.2 18 — — —
FA1(1) 1048.1 1084.1 18 1038.0 1076.0 19
AMMI 1002.8 1038.8 18 1027.1 1067.1 20
FA(1) 983.9 1049.9 33 965.7 1033.7 34
UN 825.2 1099.2 137 — — —
q为模型中方差协方差参数的数目。−2 lnL = −2 × 极大似然估计值的自然对数。AIC为 Akaike信息指数。— 表示拟合过程不收敛而
无结果。下划线的数值为最佳方差协方差结构模型对应的 AIC值。
q = No. of variance-covariance parameter. −2 lnL = −2 × log-likelihood. AIC is the Akaike’s information criterion. — denotes no results be-
cause of convergence problem. AIC values of the optimally-fitting model are underlined.
品种产量排序及其与标准对照品种产量差异显著性
比较也是区试分析的一项重要内容。在此以 2005年
第四组和 2006年第三组 2组区试资料分析的结果为
依据, 进一步探讨 ANOVA 模型和最佳方差协方差
结构线性混合模型在这方面的差异。2005年第四组
和 2006 年第三组的资料分别代表了平衡和非平衡
数据两种情况, 同时也代表了前文分析所呈现的两
种主要的方差协方差结构模型 [ 2 0 0 5 年主要为
AMMI 结构, 2006 年主要为 FA(1)结构]。表 3 是这
两组区试资料经 ANOVA 模型和最佳方差协方差结
构线性混合模型分析得到的、与对照品种具有显著
差异(包括高于和低于对照)的品种及其相应的产量
差异值。就 2005年第四组试验来看, 与对照差异显
著的品种数目在两种模型中虽然相同, 但品种及其
排序有所不同。例如品种丹 4245 和冀农 9934 在线
性混合模型(LMM)中差异显著, 而在 ANOVA 模型
中则不显著, 天泰 16 在线性混合模型中排列第二,
但在 ANOVA 模型中则排列第三。就数据不平衡的
2006 年第三组试验来看, 两种模型分析结果的差异
更突出, 不仅差异显著的品种及其排序有明显不同,
1986 作 物 学 报 第 35卷
表 2 不同模型对各组区试品种 t-测验效应差异达到显著与极显著的品种对数目及两模型间测验结果的一致率
Table 2 The number of significant and very significant variety contrasts of t-test for the trial groups and the test
consistence between two models
2005 2006
α = 0.05 α = 0.01 α = 0.05 α = 0.01
Consistence Consistence Consistence Consistence
试验组
Trial groups
ANOVA LMM
No. %
ANOVA LMM
No. %
ANOVA LMM
No. %
ANOVA LMM
No. %
一组 Group 1 54 58 48 82.8 34 27 24 70.6 53 50 48 90.6 34 35 30 85.7
二组 Group 2 54 52 50 92.6 40 38 35 87.5 49 51 47 92.2 37 37 33 89.2
三组 Group 3 49 51 47 92.2 36 38 35 92.1 49 54 42 77.8 35 39 32 82.1
四组 Group 4 57 63 53 84.1 42 43 40 93.0 66 56 54 81.8 44 48 43 89.6
平均 Average 53.5 56.0 49.5 87.9 38.0 36.5 33.5 85.8 54.3 52.8 47.8 85.6 37.5 39.8 34.5 86.7
ANOVA和 LMM分别表示方差分析模型和最佳方差协方差结构线性混合模型。α = 0.05和 α = 0.01分别为显著和极显著水平。2005
年一组、2006年一、二和四组各有 17个品种,2005年三、四组和 2006年三组各有 16个品种,2005年二组有 15个品种。品种进行两两比
较时,具有 15、16和 17个品种的试验分别有 15×14/2=105、16×15/2=120和 17×16/2=136对品种比较。Consistence:一致性。
ANOVA and LMM denote the model of analysis of variance and the linear mixed model with optimal variance covariance structure, respec-
tively. α = 0.05 and α = 0.01 mean significant and very significant, respectively. Each trial of group 1 in 2005, group 1, 2, and 4 in 2005 has 17 varie-
ties. Each trial of group 3, 4 in 2005, group 3 in 2006 has 16 varieties. The trial of group 2 in 2005 has 15 varieties. The trial with 15, 16, and 17 va-
rieties has 15×14/2=105, 16×15/2=120 and 17×16/2=136 variety contrasts, respectively.
表3 ANOVA模型和LMM模型对两组区试分析得到与对照品种相比t-测验差异显著的品种及其排序、效应差值和误差
Table 3 Name, Rank, difference value and error of the varieties with significant difference to CK for two trial groups
using ANOVA and LMM
ANOVA LMM
较对照 Compared to CK 较对照 Compared to CK 排序
Rank 品种
Variety 差值Difference 误差 Error
品种
Variety 差值 Difference 误差 Error
2005第四组 Group 4
1 鑫丰 1 Xinfeng 1 1.803 0.293 鑫丰 1 Xinfeng 1 1.599 0.260
2 东 3008 Dong 3008 0.783 0.293 天泰 16 Tiantai 16 0.680 0.247
3 天泰 16 Tiantai 16 0.707 0.293 东 3008 Dong 3008 0.652 0.278
4 辽 526 Liao 526 0.665 0.293 丹 4245 Dan 4245 0.492 0.240
5 LD9025 0.582 0.293 冀农 9934 Jinong 9934 −0.594 0.240
6 奥试 3310 Aoshi 3310 −0.648 0.293 奥试 3310 Aoshi 3310 −0.648 0.292
2006第三组 Group 3
1 明玉 2 Mingyu 2 2.540 0.383 明玉 2 Mingyu 2 2.466 0.354
2 辽 526 Liao 526 2.024 0.383 辽 526 Liao 526 2.042 0.302
3 丹 689 Dan 689 1.867 0.383 丹 689 Dan 689 1.731 0.355
4 丹玉 39 Danyu 39 1.827 0.540 沈 1083 Shen 1083 1.521 0.302
5 沈 1083 Shen 1083 1.652 0.383 K29 1.484 0.319
6 LD6037 1.605 0.383 LD6037 1.435 0.330
7 K29 1.548 0.383 DH3726 1.385 0.275
8 郑 958 Zheng 958 1.507 0.517 洛玉 1 Luoyu 1 1.316 0.295
9 洛玉 1 Luoyu 1 1.314 0.383 LD9032 1.232 0.430
10 利民 303 Limin 303 1.304 0.383 丹玉 39 Danyu 39 1.175 0.382
11 LD9032 1.225 0.383 利民 303 Limin 303 1.162 0.302
12 DH3726 1.216 0.383 京科 515 Jingke 515 0.964 0.312
13 京科 515 Jingke 515 1.062 0.383
14 承 362 Cheng 362 0.843 0.383
ANOVA和 LMM分别表示方差分析模型和最佳方差协方差结构线性混合模型。各组试验的对照品种均为农大 108。
ANOVA and LMM denote the model of analysis of variance and the linear mixed model with optimal variance covariance structure, respec-
tively. CK for each trial is Nongda 108.
第 11期 胡希远等: 基于协方差阵结构优选的作物品种区域试验分析 1987
而且差异显著的品种数目也不同。由表 3 还可以看
出, 试验品种与对照品种比较的差异值及其比较的
误差在两种模型分析中有所不同, 即两种模型对品
种效应估计的准确度不同。但一个普遍的现象是 ,
最佳方差协方差结构线性混合模型中品种效应比较
的误差要小于 ANOVA 模型中对应品种效应比较的
误差。对其余试验的分析(结果略)也表明, 最佳方差
协方差结构线性混合模型中品种比较的平均误差不
同程度地小于 ANOVA 模型中的。说明具有最佳方
差协方差结构的线性混合模型比方差分析模型对品
种效应估计更准确。
3 讨论
3.1 关于区试品种效应估计
传统方差分析采用试验观测值的算术均值估计
品种效应, 只有在试验误差方差同质、随机效应不
相关和数据平衡时, 其效应估计才具有普通最小二
乘估计(ordinary least square estimation)的优良特性。
该估计方法相当于 R 为元素值相等的对角矩阵时线
性混合模型分析的固定效应估计。文献中通常述及
的广义最小二乘估计(generalized least square esti-
mation), 或称为加权最小二乘估计 (weighted least
square estimation), 则是针对剩余误差方差不同质时
品系效应的估计, 但仍要求剩余误差和随机效应都
不相关(协方差为零), 相当于 R 和 G 的结构为元素
值不尽相等的对角矩阵时线性混合模型分析的固定
效应估计。具有方差协方差结构的线性混合模型分
析, 则在方差分析或广义最小二乘估计的条件不能
满足时模型固定效应估计也具有最佳线性无偏估计
(BLUE)的优良特性。
实际区试分析中会存在 3 个方面的问题难以满
足方差分析的条件。首先, 各试验点的环境条件(如
气候、土壤等)、甚至试验操作情况会有所不同, 要
保证各试验点的误差同质是很难实现的; 其次, 各
个品种稳定性会有所不同, 存在一些品种比另一些
品种环境间方差大的情况, 各品种方差同质的假设
难以符合实际情况; 此外, 经常存在一些品种对一
种(或多种)环境因子或病虫害因素的反应与另外的
品种不一致的情况(例如抗旱和不抗旱, 抗病和不抗
病等), 使各品种间存在不尽相同的相关性, 各对品
种间协方差同质的假设也难以符合实际情况。本文
对 2 年 8 组区试数据拟合分析的结果也证明, 方差
分析模型普遍不如其他异质性方差协方差结构线性
混合模型对数据拟合效果好。因此, 在区试分析中,
很有必要依据方差协方差结构线性混合模型分析的
原理进行品种效应估计与评价。
3.2 关于品种方差协方差结构
本文采用的几种方差协方差结构, 除 UN 结构
外, 都可在已有的几种区试分析模型中找到对应的
模型形式或假设, 易于理解和接受。被称作自由结
构(或无结构)的 UN方差协方差结构, 没有任何假设
条件, 理论上应是适用性最广泛的结构。但由于该
结构中待估计方差协方差参数的数目较多, 在一定
样本数据条件下参数估计准确性低, 实际应用中不
一定是最佳的选项。其余 4种方差协方差结构(不包
括 CS 结构)与 UN 结构的共同点在于品种方差和协
方差二者之一或全部是异质性结构, 不同点在于对
方差协方差结构异质性的构成(特性)进行了各种简
化或限制, 显著地减少了参数数目, 参数数目少的
模型结构更适用于数据样本有限的试验分析。所以,
在本文分析结果中, UN结构只对 8组区试中的一组
试验资料拟合效果最佳 , 对其余 7 组试验资料则
AMMI 结构或 FA(1)结构拟合效果最佳, 也就是说,
较简单的模型结构有时比复杂模型结构对试验数据
拟合效果要好。除了本文介绍的几种方差协方差结
构外, 还有更多其他的方差协方差结构, 它们都可
看作是 UN 结构的特定形式, 分别适用于不同特性
试验数据的线性混合模型分析。加强对它们在区试
分析应用的研究将有助于区试分析模型的发展。
3.3 关于方差协方差结构的选择
虽然人们已提出了多种模型用于区试分析, 但
由于品种对环境反应的复杂性以及各试验数据样本
大小和试验条件不同等原因, 目前人们还不能从理
论上断定哪种模型更适用一个具体区试资料的分
析。应用方差协方差结构线性混合模型分析的优点
是, 不仅可以应用各种方差协方差结构满足不同区
试数据分析的需要, 而且可以依据线性混合模型分
析的方法对各种结构模型进行评价和选择最佳模
型。
本文分析结果表明, 试验数据方差协方差不仅
不同质 , 而且其不同质的特性在各试验间有差异 ;
不仅同一方差协方差结构模型在各试验数据拟合结
果不同, 而且也没有一种结构模型在所有试验拟合
均最佳。所以不加选择地应用某一种模型, 如方差
分析模型, 去分析一区试资料, 难以保证分析结果
具有较高的准确性。本文分析结果中方差分析模型
与最佳方差协方差结构线性混合模型在品种效应估
1988 作 物 学 报 第 35卷
计及其差异显著性测验等方面呈现的明显差异, 证
明了不同结构模型分析结果的不同。这就需要人们
依据实际可操作的一定指标, 如本文所用模型拟合
信息量指标, 选择更符合试验数据实际特征的方差
协方差结构线性混合模型进行试验分析。
作为模型评价和选择的依据, 文献中有采用交
叉验证的方法。但该方法利用的毕竟是相对较繁杂
的计算机模拟手段, 不适宜一般实际试验分析者的
应用。本文尝试用信息量准则和似然比测验对线性
混合模型方差协方差结构进行评价与选择, 为区试
分析模型评价与选择问题的解决提供了一个可操作
的方法。该方法所需极大似然值统计量, 由线性混
合模型分析中极大似然估计方法直接提供, 即使对
于普通试验分析者应用也极方便。该方法的可靠性
在许多领域已得到大量验证, 但在区试分析领域对
模型评价与选择的可靠性, 还需要进一步模拟研究
的验证。
3.4 关于品种效应假设
本文对通常假设品种是固定效应的情形进行了
线性混合模型分析的探讨。这时品种效应的估计为
最佳线性无偏估计。近年来也不断有研究者提出[22-23],
可将试验品种看作随机效应, 利用线性混合模型的
最佳线性无偏预测来评价区域试验中的品种效应。
实际上在环境效应看作随机时, 品种与环境互作效
应也需要通过最佳线性无偏预测来估计与评价。关
于线性混合模型最佳线性无偏预测在区试分析中的
合理性与准确性作者将另作报道。
4 结论
方差分析法中试验因素效应方差协方差同质的
假设不能广泛准确地体现作物品种区试中品种效应
的特性; 以该模型对品种效应估计、排序和效应差
异性测验的结果与最佳方差协方差结构线性混合模
型分析的结果有着明显的不同; 线性混合模型分析
通过模型方差协方差结构的选择能客观反映区试中
品种效应的特性, 因而用于区试分析时研究结论会
更可靠。
References
[1] Wu T-X(吴天侠), Gai J-Y(盖钧镒), Ma Y-H(马育华). An in-
quiring into the nearest neighbouring analysis for trails with large
number of lines in plant breeding. Acta Agron Sin (作物学报),
1995, 21(1): 9–16 (in Chinese with English abstract)
[2] Mo H-D(莫惠栋). Statistics for Agricultural Experiments (农业
试验统计), 2nd edn. Shanghai: Shanghai Scientific and Technical
Publishers, 1992. pp 260–278 (in Chinese)
[3] Hu X-Y(胡希远), Spilke J. Spatial variability and its statistical
control in field experiment. Acta Agron Sin (作物学报), 2007,
33(4): 620–624 (in Chinese with English abstract)
[4] Kempthorne O. The Design and Analysis of Experiments. New
York: John Wiley & Sons, 1952. pp 86–156
[5] Zhang Q-Y(张群远), Kong F-L(孔繁玲). Models and methods
for estimating variety means in regional crop trials: Comparisons
of arithmetic mean, weighted least squares estimates and BLUPs.
Acta Agron Sin (作物学报), 2003, 29(6): 884–891 (in Chinese
with English abstract)
[6] Wu Y-Q(吴元奇), Pan G-T(潘光堂), Rong T-Z(荣廷昭). Study
progress in crop stability. J Sichuan Agric Univ (四川农业大学
学报), 2005, 23(4): 482–489 (in Chinese with English abstract)
[7] Shukla G K. Some statistical aspects of partitioning genotype
environmental components of variability. Heredity, 1972, 29:
237–245
[8] Finly K W, Wilkinson G N. The analysis of adaptation in a plant
breeding program. Aust J Agric Res, 1963, 14: 742–754
[9] Eberhart S A, Russell W A. Stability parameters for comparing
varieties. Crop Sci, 1966, 6: 36–40
[10] Gauch H G. Model selection and validation for yield trials with
interaction. Biometrics, 1988, 44: 705–715
[11] Wang L(王磊), Yang S-H(杨仕华), Xie Y-X(谢英贤). AMMI
model and its application to the regional crop trial analysis. J
Appl Fund Eng Sci (应用基础与工程科学学报), 1997, 5(1):
39–46 (in Chinese with English abstract)
[12] Oman S D. Multiplicative effects in mixed model analysis of
variance. Biometrika, 1991, 78: 729–739
[13] Zhang Q-Y(张群远), Kong F-L(孔繁玲). Comparison of statisti-
cal models for regional crop trial analysis. Sci Agric Sin (中国农
业科学), 2002, 35(4): 365–371 (in Chinese with English abstract)
[14] Little R C, Milliken G A, Stroup W W, Wolfinger R D. SAS Sys-
tem for Mixed Models. Cary, NC: SAS Institute, Inc. 1996. pp
303–326
[15] Verbeke G, Molenberghs G. Linear Mixed Models in Practice.
New York: Springer-Verlag, 1997. pp 185–210
[16] Hartley H O, Rao C R. Maximum likelihood estimation for the
mixed analysis of variance model. Biometrica, 1967, 54: 93–108
[17] Harville D A. Maximum likelihood approaches to variance com-
ponent estimation and to related problems. J Am Stat Assoc, 1977,
72: 320–338
[18] Zhu J(朱军). Principles of Linear Model Analysis (线性模型分
析原理). Beijing: Science Press, 2000. pp 160–166 (in Chinese)
[19] Wolfinger R. Covariance structure selection in general mixed
models. Commun Stat Simula Comput, 1993, 22: 1079–1106
[20] Piepho H P. Stability analysis using the SAS system. Agron J,
1999, 91: 154–160
[21] Bozdogan H. Akaike’s information criterion and recent develop-
ments in information complexity. J Math Psychol, 2000, 44:
62–91
[22] Piyepho H P, Mohring J. Best linear unbiased prediction of culti-
var effects for subdivided target region. Crop Sci, 2005, 45:
1151–1159
[23] Smith A B, Cullis B R, Thompson R. The analysis of crop culti-
var breeding and evaluation trials: An overview of current mixed
model approaches. J Agric Sci, 2005, 143: 449–462
第 11期 胡希远等: 基于协方差阵结构优选的作物品种区域试验分析 1989
附录:
假定用 location、variety、block和 yield分别表示试点、
品种、区组和产量 4 个变量,则利用 SAS 对区域试验进行
CS方差协方差结构线性混合模型,分析的程序如下:
/* 区试数据建立为 SAS数据集的程序段 */
data a;
input location$ variety$ block yield;
cards;
dalian nongda108 1 12.3
dalian nongda108 2 11.8
dalian nongda108 3 10.7
dalian dan689 1 10.8
dalian dan689 2 9.9
dalian dan689 3 9.8
dadong nongda108 1 11.0
……
; run;
/* 数据预处理程序段 */
proc sort data=a out=a_s;
by location; run;
proc mixed data= a_s;
class variety block;
model yield= variety block/ ;
lsmeans variety;
ods output lsmeans=lsm;
by location; run;
data lsm; set lsm;
y=estimate;
var_mean=Stderr**2;
w=1/var_mean; run;
/* CS方差协方差结构线性混合模型分析的程序段 */
proc mixed data = lsm;
class variety location;
weight w;
model y = variety / ddfm = satterthwaite;
random variety / subject = location type = cs;
lsmeans variety /pdiff; run;
前两个程序段中各语句与试验分析常规 SAS 程序中
语句的作用和意义相同,其中数据预处理程序段的作用是
对各试点分别进行方差分析,求出各点的品种效应估计值
和误差。第三个程序段中,语句 weight w; 的作用是要求
在数据分析时以w为权重处理各试点误差的异质性。ddfm
= satterthwaite指定计算自由度的方法为 satterthwaite法。
语句 random variety / subject = location type = cs; 的作用
是指出(试点)环境效应与品种-环境互作效应为随机效应
定义方差协方差的结构类型。语句 lsmeans variety/ pdiff
要求给出所有品种效应估计值及两两差异显著性 t-测验
的结果。proc mixed程序默认的方差协方差参数估计方法
是限制性极大似然法(REML)。此处 CS 结构的程序分析
结果相当于加权最小二乘估计的结果。如果将该程序中的
语句 weight w; 去掉,则该程序就是区试数据普通方差分
析的程序。
若要进行其他方差协方差结构线性混合模型的分析,
只要将上述第三程序段中的语句 random variety / subject =
location type = cs; 利用下面各种方差协方差结构对应的一
个或两个 random 语句替代即可。各程序间的区别仅在于利
用一个或两个 random 语句定义各自对应的方差协方差结
构。根据各种模型得到的 AIC值(程序运行后自动给出)选出
一个 AIC值最小模型的分析结果作为最优的分析结果。
/* UN(1)结构 */
random variety / subject = location type = un(1);
random variety / subject = location;
/* FA1(1)结构 */
random variety / subject = location type = fa1(1);
/* AMMI结构 *
random variety / subject = location type = fa1(1);
random variety / subject = location;
/* FA(1)结构 */
random variety / subject = location type = fa(1);
random variety / subject = location;
/* UN结构 */
random variety / subject = location type = un;