四階巢式設(shè)計(jì)分組隨區(qū)組試驗(yàn)的非平衡數(shù)據(jù)處理的理論研究

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

（1.中國林業(yè)科學(xué)研究院亞熱帶林業(yè)研究所，浙江 杭州 311400;2.浙江省遂昌縣生態(tài)林業(yè)發(fā)展中心，浙江 遂昌 323300）

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

假定一個(gè)樹種，從其分布區(qū)中抽取5～8 個(gè)典型種源，每個(gè)種源內(nèi)抽取5～8 個(gè)林分，每個(gè)林分中選取5～8株優(yōu)樹，分系采集每個(gè)優(yōu)株的OP 種子，進(jìn)行育苗造林試驗(yàn)，采用分組隨機(jī)區(qū)組設(shè)計(jì)，同一種源放在一起，占領(lǐng)一個(gè)亞區(qū)組，不同種源隨機(jī)排列，十個(gè)重復(fù)，4 株行狀小區(qū)。采用我們介紹的方法，先用每個(gè)重復(fù)內(nèi)的環(huán)境指數(shù),對(duì)各亞區(qū)組參試植株觀察值進(jìn)行調(diào)整[11]，再將每個(gè)重復(fù)內(nèi)的k 株小區(qū)轉(zhuǎn)化為k 個(gè)重復(fù)的單株小區(qū)試驗(yàn)，進(jìn)而采用轉(zhuǎn)化分析法構(gòu)建的隨機(jī)模型進(jìn)行數(shù)據(jù)處理[10]。

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

1 轉(zhuǎn)化后的線性模型:

以單株觀察值參與統(tǒng)計(jì)分析時(shí)，平衡不平衡、規(guī)則不規(guī)則的線性模型如下：

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

對(duì)于以上線性模型，在隨機(jī)模型條件下，有如下約束條件：

（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 是群體平均效應(yīng);bi 是重復(fù)效應(yīng);pj是種源效應(yīng);sk/j是第j 個(gè)種源內(nèi)第k 個(gè)林分效應(yīng)；fl/k是第k 個(gè)林分中第l 個(gè)家系的隨機(jī)效應(yīng);(Pb)ij是種源與重復(fù)的互作效應(yīng);(sb)ik是林分與重復(fù)間的互作；eijkl是隨機(jī)誤差。所有的因子在隨機(jī)條件下進(jìn)行試驗(yàn),因此有: 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]的做法,對(duì)上式取數(shù)學(xué)期望，因此對(duì)以上等式有：所有參試因子兩兩之積的交叉項(xiàng)為0,于是有下式:

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

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

3 各因子效應(yīng)的離差平方和

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

4 期望均方結(jié)構(gòu)

從線性模型出發(fā)，推導(dǎo)期望均方結(jié)構(gòu)，先對(duì)各參試因子的參試子代數(shù)加以說明。

nijkl表示第i 個(gè)重復(fù)內(nèi)第j 個(gè)種源內(nèi)第k 個(gè)林分內(nèi)第l 個(gè)家系的參試子代數(shù)；nijkl=1 或0；

ni…表示第i 個(gè)重復(fù)內(nèi)的參試子代數(shù)；

n.j..表示第j 個(gè)種源內(nèi)的參試子代數(shù)；

n..jk.表示第j 個(gè)種源內(nèi)第k 個(gè)林分因子內(nèi)的參試子代數(shù)；

n..jkl表示第j 個(gè)種源內(nèi)第k 個(gè)林分因子內(nèi)第l 個(gè)家系內(nèi)的參試子代樣本數(shù)；

nij..表示第j 個(gè)種源在第i 個(gè)重復(fù)內(nèi)的參試子代數(shù)；

nijk.表示第j 個(gè)種源內(nèi)第k 個(gè)林分在第i 個(gè)重復(fù)內(nèi)的參試子代數(shù)；

N....表示試驗(yàn)林中全部的參試子代數(shù)；

轉(zhuǎ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

根據(jù)線性模型,可推導(dǎo)出表1中的結(jié)果。

隨機(jī)模型條件下，各參試因子的方差系數(shù)如表1。

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

5 自由度的分解

表2 自由度分解

6 期望均方結(jié)構(gòu)

四階不平衡巢式設(shè)計(jì)模式轉(zhuǎn)化分析法的期望均方結(jié)構(gòu)列于表3。

表3 四階不平衡巢式設(shè)計(jì)模式轉(zhuǎn)化分析法的期望均方結(jié)構(gòu)Table 3 Expected mean square structure of transformation analysis for four-order unbalanced nested design

這里: 期望均方=期望平方和/自由度，方差調(diào)節(jié)系數(shù)K 值從期望平方和表1中計(jì)算而來。

7 主效因子的F 檢驗(yàn)

第一步，根據(jù)期望均方結(jié)構(gòu)，建立聯(lián)立線性方程組，求出σ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

第二步，運(yùn)用行列式知識(shí)求解（可用矩陣法求解，但是為了進(jìn)行F 檢驗(yàn)，用行列式法還是必須的）,如下:

于是σ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

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

σ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為已知系數(shù),這樣做的目的是為了在進(jìn)行主效因子的F 檢驗(yàn)時(shí),配合參比均方,計(jì)算出該均方的自由度。

例，在對(duì)種源效應(yīng)作顯著性檢驗(yàn)時(shí)，σ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的線性函數(shù)，并代入MSx中，再合并同類項(xiàng)，再展成MSb，MSp，MSs/p，MSf/s，MSpb，MSsb，MSt的線性函數(shù)，

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

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

因此有下式：

有了Nx，便可以根據(jù)查F 表，對(duì)假設(shè)前提，作出判斷，肯定或否定假設(shè)。

其它因子的F 檢驗(yàn)可依照此例進(jìn)行，為了節(jié)約篇幅，此處從略。

8 因子效應(yīng)的多重對(duì)比

四階不平衡巢式設(shè)計(jì)，其因子效應(yīng)的多重對(duì)比十分繁復(fù)。

這可仿三階巢式設(shè)計(jì)平衡不平衡試驗(yàn)資料分析方法中的做法進(jìn)行[10]。

9 四階不平衡巢式設(shè)計(jì)遺傳力的計(jì)算

種源遺傳力：

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

種源內(nèi)的林分遺傳力：

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

林分內(nèi)家系遺傳力：

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

家系內(nèi)單株遺傳力：

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

10 四階巢式設(shè)計(jì)分組隨區(qū)組試驗(yàn)的非平衡數(shù)據(jù)處理軟件

四階巢式設(shè)計(jì)分組隨區(qū)組試驗(yàn)的非平衡數(shù)據(jù)處理，可采用《林木遺傳育中平衡不平衡、規(guī)則不規(guī)則試驗(yàn)數(shù)據(jù)處理技巧》中的MLAP 軟件，采用R2016a 平臺(tái)進(jìn)行。