四階巢式設計分組隨區組試驗的非平衡數據處理的理論研究

齊 明，包小梅，何貴平*，張 振

（1.中國林業科學研究院亞熱帶林業研究所，浙江 杭州 311400;2.浙江省遂昌縣生態林業發展中心，浙江 遂昌 323300）

在林木遺傳改良中,為了制定正確有效的選擇育種方案,需要了解群體內不同層次的遺傳變異信息,為此林木育種工作者經常采用多階巢式設計[1-6]。1962年Charles EG 等人[7]發表了平衡數據的四階巢式設計的方差分析原理；2006年El-Saeiti I N 等人[8,9]發表了有缺失數據的四階巢式設計的方差分析原理。這些研究美中不足是沒有考慮區組效應和參試因子與區組重復間的互作,這極大地限制了這些模型在林木遺傳育種中的應用。由于四階巢式設計試驗在林木遺傳改良中還有一定的市場，例如：種源—林分—家系—個體間的變異;或林分—家系—個體—無性系試驗的變異等等。2009年齊明[10]提出了一個轉化理論，建議采用轉化分析法，構建了一個包含區組重復效應和參試因子與區組重復間互作的線性模型，開發了林木遺傳育種中四階巢式設計的方差分析原理，但其缺點是忽視了在方差分析中一級參試因子與二級參試因子，二級參試因子與三級參試因子間不獨立的事實，因此現在有必要對此研究作進一步的完善。另外本研究考慮了在林木遺傳育種中四階巢式設計的田間試驗，參試材料眾多，田間試驗通常采用分組隨機區組設計[1-6]的事實，建議使用環境指數對各亞區組內數據調整后[11]，然后采用轉化分析法的思想[10]進行非平衡（缺株）、非規則（缺區）試驗數據處理。下面以種源—林分—家系—個體四階巢式設計試驗為例，進行方差分析的研究。

假定一個樹種，從其分布區中抽取5～8 個典型種源，每個種源內抽取5～8 個林分，每個林分中選取5～8株優樹，分系采集每個優株的OP 種子，進行育苗造林試驗，采用分組隨機區組設計，同一種源放在一起，占領一個亞區組，不同種源隨機排列，十個重復，4 株行狀小區。采用我們介紹的方法，先用每個重復內的環境指數,對各亞區組參試植株觀察值進行調整[11]，再將每個重復內的k 株小區轉化為k 個重復的單株小區試驗，進而采用轉化分析法構建的隨機模型進行數據處理[10]。

根據我們的研究成果，數據轉化后，在非平衡條件下，所有的參試因子均符合線性，正態的前提條件。種源與重復互作因子肯定是遵從正態分布;而家系與區組重復間因無重復，其互作效應可以不考慮，因此線性模型中不包括此項因子;而林分與區組互作因子，在理論上遵從二項分布，有可能獲得負的方差分量，也有可能獲得正的方差分量，但為了模型的通用性，仍將其包含在模型中。如果采用此模型獲得負的方差分量，可將此項因子刪去，重構線性模型進行分析，這可仿兩因素隨機區組試驗（析因設計）的方法[10]進行。

1 轉化后的線性模型:

以單株觀察值參與統計分析時，平衡不平衡、規則不規則的線性模型如下：

yijklm=u + bi+ pj+ sk/j+fl/k+ (pb)ij+(sb)ik+ eijklm

這里: i=1→a;j=1→p;k=1→s;l=1→f;m=1 或0

對于以上線性模型，在隨機模型條件下，有如下約束條件：

（1）E（bi）=0；E（pj）=0；E（sk/j）=0；E（fj/l）=0；E（pb ）ij）=0；E (sb)ik=0；E（eijklm）=0；

（2）在上述線性模型中,u 是群體平均效應;bi 是重復效應;pj是種源效應;sk/j是第j 個種源內第k 個林分效應；fl/k是第k 個林分中第l 個家系的隨機效應;(Pb)ij是種源與重復的互作效應;(sb)ik是林分與重復間的互作；eijkl是隨機誤差。所有的因子在隨機條件下進行試驗,因此有: bi∽N(0,σb2) ;pj∽N(0,σp2);sk/j∽N(0,σs/p2);fl/k∽N(0,σf/s2);(Pb)ij∽N(0,σpb2);(sb)ik∽N(0,σsb2);eijkl∽N(0,σe2).

2 離差平方和的分解

仿照Henderson CR[12]的做法,對上式取數學期望，因此對以上等式有：所有參試因子兩兩之積的交叉項為0,于是有下式:

即，SST= SSb+ SSp+ SSs/p+ SSf/s+ SSpb+SSsb+SSe

這樣可以得如下離差平方和的分解結果。

3 各因子效應的離差平方和

上式SSe 的展開式太復雜，實際中可用SSe=SST-SSb-SSp-SSs/p-SSf/s-SSpb-SSsb來計算SSe以上各因子效應平方和中的Y 和N 分別表示因子效應的求和及其參試子代樣本數。

4 期望均方結構

從線性模型出發，推導期望均方結構，先對各參試因子的參試子代數加以說明。

nijkl表示第i 個重復內第j 個種源內第k 個林分內第l 個家系的參試子代數；nijkl=1 或0；

ni…表示第i 個重復內的參試子代數；

n.j..表示第j 個種源內的參試子代數；

n..jk.表示第j 個種源內第k 個林分因子內的參試子代數；

n..jkl表示第j 個種源內第k 個林分因子內第l 個家系內的參試子代樣本數；

nij..表示第j 個種源在第i 個重復內的參試子代數；

nijk.表示第j 個種源內第k 個林分在第i 個重復內的參試子代數；

N....表示試驗林中全部的參試子代數；

轉化后的線性模型:yijklm=u+bi+pj+ sk/j+fl/k+ (pb)ij+(sb)ik+eijklm

這里: i=1→a;j=1→p;k=1→s;l=1→f;m=1 或0

根據線性模型,可推導出表1中的結果。

隨機模型條件下，各參試因子的方差系數如表1。

表1 隨機模型條件下，四階不平衡巢式設計各因子的方差分量系數表Table 1 Variance component coefficients of each factor in fourth-order unequal nested design under random model

5 自由度的分解

表2 自由度分解

6 期望均方結構

四階不平衡巢式設計模式轉化分析法的期望均方結構列于表3。

表3 四階不平衡巢式設計模式轉化分析法的期望均方結構Table 3 Expected mean square structure of transformation analysis for four-order unbalanced nested design

這里: 期望均方=期望平方和/自由度，方差調節系數K 值從期望平方和表1中計算而來。

7 主效因子的F 檢驗

第一步，根據期望均方結構，建立聯立線性方程組，求出σb2，σp2,σs/p2，σf/s2，σpb2，σsb2，σe2方差份量。

σe2+K1σsb2+K2σpb2+K3σf/s2+K4σs/p2+K5σp2+K6σb2=MSb

σe2+K7σsb2+K8σpb2+K9σf/s2+K10σs/p2+K11σp2+K12σb2=MSp

σe2+K13σsb2+K14σpb2+K15σf/s2+K16σs/p2+K17σp2+K18σb2=MSs/p

σe2+K19σsb2+K20σpb2+K21σf/s2+K22σs/p2+K23σp2+K24σb2=MSf/s

σe2+K25σsb2+K26σpb2+K27σf/s2+K28σs/p2+K29σp2+K30σb2=MSpb

σe2+K31σsb2+K32σpb2+K33σf/s2+K34σs/p2+K35σp2+K36σb2=MSsb

σe2+K37σsb2+K38σpb2+K39σf/s2+K40σs/p2+K41σp2+K42σb2=MSt

第二步，運用行列式知識求解（可用矩陣法求解，但是為了進行F 檢驗，用行列式法還是必須的）,如下:

于是σb2=Db/D；σp2=Dp/D；σs/p2=Ds/p/D;σf/k2=Df/s/D；σpb2=Dpb/D；σsb2=Dsb/D;σe2=DE/D

為了進行主效因子的F 檢驗，除了將D 計算出行列式的具體數值外，其余的行列式是將MSb，MSp，MSs/p,MSf/s,MSpb，MSsb，MSt示為未知數，進行計算，合并同類項再將σb2，σp2,σs/p2，σf/s2，σpb2，σsb2，σe2展成MSb，MSp，MSs/p;MSf/s,MSpb，MSsb，MSt的線性函數，例：

σb2=β1MSb+β2MSp+β3MSs/p+β4MSf/s+β5MSpb+β6MSsb+β7MSt

σp2=β8MSb+β9MSp+β10MSs/p+β11MSf/s+β12MSpb+β13MSsb+β14MSt

σs/p2=β15MSb+β16MSp+β17MSs/p+β18MSf/s+β19MSpb+β20MSsb+β21MSt

σf/s2=β22MSb+β23MSp+β24MSs/p+β25MSf/s+β26MSpb+β27MSsb+β28MSt

σpb2=β29MSb+β30MSp+β31MSs/p+β32MSf/s+β33MSpb+β34MSsb+β35MSt

σsb2=β36MSb+β37MSp+β38MSs/p+β39MSf/s+β40MSpb+β41MSsb+β42MSt

σe2=β43MSb+β44MSp+β45MSs/p+β46MSf/s+β47MSpb+β48MSsb+β49MSt

上述線性方程中，β1，β2…β49為已知系數,這樣做的目的是為了在進行主效因子的F 檢驗時,配合參比均方,計算出該均方的自由度。

例，在對種源效應作顯著性檢驗時，σe2+K7σsb2+K8σpb2+K9σf/s2+K10σs/p2+K11σp2+K12σb2=MSp

F= MSp/MSx∽遵從

這里，MSx=σe2+K7σsb2+K8σpb2+K9σf/s2+K10σs/p2+K12σb2

將上述σe2，σsb2，σpb2，σf/s2，σs/p2，σb2展成MSb，MSp，MSs/p;MSf/s,MSpb，MSsb，MSt的線性函數，并代入MSx中，再合并同類項，再展成MSb，MSp，MSs/p，MSf/s，MSpb，MSsb，MSt的線性函數，

假設為：MSx=α1MSb+α2MSp+α3MSs/p+α4MSf/s+α5MSpb+α6MSsb+α7MSt

如果MSx的自由度為Nx，則MSx/Nx∽X2（Nx）

因此有下式：

有了Nx，便可以根據查F 表，對假設前提，作出判斷，肯定或否定假設。

其它因子的F 檢驗可依照此例進行，為了節約篇幅，此處從略。

8 因子效應的多重對比

四階不平衡巢式設計，其因子效應的多重對比十分繁復。

這可仿三階巢式設計平衡不平衡試驗資料分析方法中的做法進行[10]。

9 四階不平衡巢式設計遺傳力的計算

種源遺傳力：

hp2=σp2/[（1/K11）σe2+（K7/K11）σsb2+（K8/K11）σpb2+（K9/K11）σf/s2+（K10/K11）σs/p2+σp2+（K12/K11）σb2]；

種源內的林分遺傳力：

hs/p2=K16σs/p2/[σe2+K13σsb2+K14σpb2+K15σf/s2+K16σs/p2+K17σp2+K18σb2] ；也可分子分母同除K16

林分內家系遺傳力：

hf/s2=K21σf/s2/[σe2+K19σsb2+K20σpb2+K21σf/s2+K22σs/p2+K23σp2+K24σb2] ；也可分子分母同除K21

家系內單株遺傳力：

hi2=4σf/s2/[σf/s2+σb2+σe2]或3σf/k2/[σf/s2+σb2+σe2]

10 四階巢式設計分組隨區組試驗的非平衡數據處理軟件

四階巢式設計分組隨區組試驗的非平衡數據處理，可采用《林木遺傳育中平衡不平衡、規則不規則試驗數據處理技巧》中的MLAP 軟件，采用R2016a 平臺進行。