林木子代試驗中小區(qū)平均值法與轉(zhuǎn)化分析法比較

葉金俊, 葉金水, 包小梅, 黃琳, 毛海淋

1. 遂昌縣生態(tài)林業(yè)發(fā)展中心，浙江 遂昌 323300；2. 遂昌縣自然資源和規(guī)劃局，浙江 遂昌 323300

單株觀察值法常用于處理林木統(tǒng)計學(xué)研究試驗數(shù)據(jù)，不少學(xué)者[1-6]應(yīng)用Monte Carlo（蒙特卡羅）模擬法，在試驗數(shù)據(jù)不平衡或嚴重不平衡時，比較驗證不同分析方法間的分析效果、優(yōu)點和不足之處，分析其與方差分量估計法的優(yōu)劣[1-2,7-9,10-12]。

林木田間試驗由于外在環(huán)境及植物間競爭等因素，導(dǎo)致部分試驗植株死亡，最終得到的是非平衡試驗資料。按照傳統(tǒng)的現(xiàn)代線性模型理論、方差分析等來處理這些非平衡的試驗數(shù)據(jù)，常常獲得負的方差分量，不能給出數(shù)學(xué)解釋，同時也沒有生物學(xué)意義，說明試驗數(shù)據(jù)的分析方法可能有問題。2009年，齊明對林木正交試驗設(shè)計的非平衡數(shù)據(jù)處理，提出了一個轉(zhuǎn)化理論[13]，2014年，齊明采用計算機模擬，論證了轉(zhuǎn)化分析法的統(tǒng)計學(xué)基礎(chǔ)，并得出就同一非平衡數(shù)據(jù)，轉(zhuǎn)化分析法無論是在隨機模型，還是在固定模型上均優(yōu)于Henderson方法I的結(jié)論[14,15]。但小區(qū)平均值法處理林木田間試驗數(shù)據(jù)的現(xiàn)象仍然隨處可見，并且數(shù)據(jù)量越大（如多地點多年度的試驗數(shù)據(jù)），采用小區(qū)平均值法的頻度就越高。知網(wǎng)文獻檢索結(jié)果顯示:至2021年2月采用小區(qū)平均值法來處理林木田間試驗的非平衡數(shù)據(jù)的文獻多達8 943條。 林木遺傳育種中，無性系材料的田間試驗資料，可采用小區(qū)平均值法進行數(shù)據(jù)處理；林木遺傳育種的正交田間試驗數(shù)據(jù)，由于平衡不平衡（缺株）、規(guī)則不規(guī)則（缺區(qū)），宜以單株值進行轉(zhuǎn)化分析法統(tǒng)計分析[13-14]。為系統(tǒng)地評估小區(qū)平均值法在處理林木子代試驗非平衡數(shù)據(jù)中的優(yōu)缺點，借助Monte Carlo模擬法，產(chǎn)生非平衡試驗數(shù)據(jù)，以轉(zhuǎn)化分析法分析結(jié)果為參照對象，評估了林木試驗中小區(qū)平均值法的統(tǒng)計效率，分析精確性、參數(shù)大小和精度，比較研究了這兩種方法在隨機模型和固定模型條件下的分析后果: （1）著重評估數(shù)據(jù)處理方法能否達到田間試驗的目的;（2）其統(tǒng)計分析的效率如何;（3）其統(tǒng)計分析的精確性如何。為選用正確而合適的統(tǒng)計分析方法和促進林木數(shù)量遺傳學(xué)發(fā)蔚縣提供科學(xué)依據(jù)。

1 模型與方法

1.1 平衡試驗設(shè)計及其統(tǒng)計分析模型

從研究的普遍性、計算工作量的大小等方面考慮，采用單因素隨機區(qū)組（RCB）設(shè)計。林木子代測定的田間試驗中,小區(qū)的設(shè)計一般為4～6株小區(qū)、也有8或10株小區(qū)。以下五個試驗，采用Monte Carlo法模擬。為體現(xiàn)研究性狀內(nèi)在的遺傳變異幅度,采用了40個半同胞家系，每個試驗分別產(chǎn)生100套數(shù)據(jù):

（Ⅰ）40個半同胞家系，10個區(qū)組，4株小區(qū)參試；

（Ⅱ）40個半同胞家系，8個區(qū)組，5株小區(qū)參試；

（Ⅲ）40個半同胞家系，8個區(qū)組，6株小區(qū)參試；

（Ⅳ）40個半同胞家系，6個區(qū)組，8株小區(qū)參試；

（Ⅴ）40個半同胞家系，4個區(qū)組，10株小區(qū)參試。

1.2 平衡試驗數(shù)據(jù)的產(chǎn)生

按照林木遺傳育種的經(jīng)驗，假設(shè)半同胞家系間、區(qū)組重復(fù)、小區(qū)株間因子服從正態(tài)分布，指定株間變異的方差Ve、區(qū)組重復(fù)因子的方差Vb、半同胞家系的方差Vf分別為280、15 、20。

按原試驗設(shè)計模式（即Henderson方法I）來模擬試驗數(shù)據(jù)：處理與區(qū)組重復(fù)間的交互作用因子不遵循正態(tài)分布，繪散點圖表明（FB）ij遵從二項分布，因此按二項分布產(chǎn)生其效應(yīng)值。群體平均值為60。采 用Matlab7.X語 言[19,20]編 寫 程 序,以 獲取RCB設(shè)計的試驗平衡資料。每個試驗一次模擬共100套數(shù)據(jù)，保存至EXCEL文件。

1.3 非平衡試驗數(shù)據(jù)的產(chǎn)生及統(tǒng)計分析

由于氣候環(huán)境、造林技術(shù)和株間競爭等多方面因素對造林存活率有所影響，多年來林木測定調(diào)查結(jié)果為造林存活率多數(shù)在75%～95%，多片10多年生杉木試驗林的平均保存率在83%。假設(shè)隨機死亡的保存率為83%，通過MATLAB語言中設(shè)計的刪除語句隨機刪除,獲取非平衡試驗數(shù)據(jù)。

1.4 不平衡試驗數(shù)據(jù)的轉(zhuǎn)化分析法分析

由于是正交試驗，故每套資料先將多株小區(qū)轉(zhuǎn)化為單株小區(qū)[13]，采用轉(zhuǎn)化分析法建立如下線性模型：

其中： i = 1→40； j =1→40～48；k=1或0。u是群體平均效應(yīng)；Fi是第i個家系的效應(yīng)值；Bj是第j個區(qū)組重復(fù)的效應(yīng)值；Eijk是株間變異；Yijk是第i個家系在第j個區(qū)組重復(fù)中第k個觀察值。

在隨機模型條件下，各參試因子都是隨機因子。于是：

方差分析原理參見有關(guān)文獻[13]。期望均方結(jié)構(gòu)見表1。

表中的期望均方的調(diào)節(jié)系數(shù)k計算、非平衡數(shù)據(jù)的處理和參數(shù)估計見參考文獻[13,14]。

1.5 不平衡試驗數(shù)據(jù)的小區(qū)平均值分析法

使用Matlab7.X語言，編寫以獲得以小區(qū)平均值為單位,參與統(tǒng)計分析非平衡數(shù)據(jù)的程序。林木半同胞子代試驗，以小區(qū)平均值進行統(tǒng)計分析時的線性模型為：

以小區(qū)平均值參與計算時，原試驗變成了兩因素?zé)o重復(fù)的方差分析模型，其期望均方結(jié)構(gòu)如下表2。

根據(jù)表2中的結(jié)果，估算隨機模型條件下，因子方差分量和家系遺傳力；然后分析固定模型條件下，家系效應(yīng)值的大小和秩的次序，以轉(zhuǎn)化分析法的結(jié)果為參照物來分析小區(qū)平均值法的家系效應(yīng)值產(chǎn)生的誤差及選擇產(chǎn)生的失誤率大小。

表 1 單因素隨機區(qū)組不平衡數(shù)據(jù)轉(zhuǎn)化分析法的期望均方結(jié)構(gòu)Tab. 1 Expected mean square structure of single factor random block non-equilibrium data in transformation analysis method

表 2 半同胞試驗小區(qū)平均值法的期望方差結(jié)構(gòu)Tab. 2 Expected variance structure of plot average method in half-sib experimental plot

在此模型中，半同胞家系遺傳力為：hf2=σf2/[σf2+(1/b)Ve]

由于缺乏株間變異性信息，無法計算單株遺傳力，許多文獻采用如下公式：

上述公式中Ve為小區(qū)間的誤差，故該計算公式結(jié)果有誤。

1.6 以Monte Carlo模擬數(shù)據(jù)為基礎(chǔ)的小區(qū)平均值法的結(jié)果評價

轉(zhuǎn)化分析法：根據(jù)參試因子的調(diào)節(jié)系數(shù)和參試因子的方差分量，計算半同胞家系遺傳力 Hf2和家系內(nèi)單株遺傳力hi2；小區(qū)平均值法：根據(jù)平衡模型,估計方差分量和家系遺傳力。

根據(jù)100次Monte Carlo模擬結(jié)果，計算出若干參數(shù)的估計值，參數(shù)的偏差（偏差=估計值-參數(shù)真值）、偏差的顯著性檢驗[1-2,7-9]。如果偏差的絕對值被其估計值除得的商，大于5%，可認為該參數(shù)有偏（Graybill & wortham，1956）。均方誤差為MSE=Var（估計值）+偏差2。在相同參數(shù)情況下，均方誤差越小，估計效益和精度越高[2,3-6,7-9,21]。

所有計算分析均在MATLAB7.0[19-20]和Excel 2003平臺上完成。

2 結(jié)果與分析

2.1 隨機模型條件下，小區(qū)平均值法的評估結(jié)果

（1）五個模擬試驗中兩種方法獲得參試因子負方差分量的頻率

五個試驗的非平衡試驗數(shù)據(jù),林木轉(zhuǎn)化分析法和小區(qū)平均值法的分析結(jié)果列于表3.

由表3可見, 五個試驗的非平衡試驗數(shù)據(jù)，林木轉(zhuǎn)化分析法均未獲得負的方差分量；而小區(qū)平均值法中，試驗Ⅵ、Ⅴ中，區(qū)組因子分別有3%、6%的試驗有負的方差分量。結(jié)果支持了轉(zhuǎn)化分析法中參試因子是隨機正態(tài)因子的觀點。只有滿足方差分析的正態(tài)性前提條件，才不會有負的方差分量[13,15]，同時表明轉(zhuǎn)化分析法構(gòu)建的線性模型是科學(xué)的、轉(zhuǎn)化分析法優(yōu)越于小區(qū)平均值法。

（2）幾個參數(shù)效益的比較

在隨機死亡的前提下，采用MATLAB語言的刪除指令，將平衡的試驗數(shù)據(jù)轉(zhuǎn)變成不平衡的試驗數(shù)據(jù)。用轉(zhuǎn)化分析法來處理試驗資料，計算小區(qū)平均值，再采用小區(qū)平均值法來分析試驗數(shù)據(jù)，其結(jié)果列于表4～8。

試驗‖轉(zhuǎn)化分析法和小區(qū)平均值法的比較分析結(jié)果列于表4。

由表4可見，轉(zhuǎn)化分析法的參數(shù)偏差不顯著，小區(qū)平均值法中，誤差方差、家系方差和區(qū)組方差偏性顯著，且其均方誤差也比轉(zhuǎn)化分析法的誤差大。

表 3 五個模擬試驗中兩種方法獲得參試因子負方差分量的頻率（%）Tab. 3 frequencis of negative variance component of the test factors obtained by two methods in five simulated experiments

表 4 試驗Ⅱ轉(zhuǎn)化分析法與小區(qū)平均值法分析結(jié)果的比較Tab. 4 Comparison of analysis results between transformation analysis method and plot average method in experiment Ⅱ

試驗Ⅰ轉(zhuǎn)化分析法和小區(qū)平均值法的比較分析結(jié)果列于表5。

試驗Ⅲ中,轉(zhuǎn)化分析法和小區(qū)平均值法的比較分析結(jié)果列于表6。

試驗Ⅳ中,轉(zhuǎn)化分析法與小區(qū)平均值法分析結(jié)果的比較列于表7。

試驗Ⅴ中,轉(zhuǎn)化分析法與小區(qū)平均值法分析結(jié)果的比較列于表8。

表 5 試驗Ⅰ轉(zhuǎn)化分析法與小區(qū)平均值法分析結(jié)果的比較Tab. 5 Comparison of analysis results between transformation analysis method and plot average method in experiment Ⅰ

表 6 試驗Ⅲ轉(zhuǎn)化分析法與小區(qū)平均值法分析結(jié)果的比較Tab. 6 Comparison of analysis results between transformation analysis method and plot average method in experiment Ⅲ

表 7 試驗Ⅳ轉(zhuǎn)化分析法與小區(qū)平均值法分析結(jié)果的比較Tab. 7 Comparison of analysis results between transformation analysis method and plot average method in experiment Ⅳ

表 8 試驗Ⅴ轉(zhuǎn)化分析法與小區(qū)平均值法分析結(jié)果的比較Tab. 8 Comparison of analysis results between transformation analysis method and plot average method in experiment Ⅴ

從表5～8可見，轉(zhuǎn)化分析法，各參數(shù)都存在一定的偏差，但均未達到顯著水平，參試因子的方差分量估計值可認為無偏的。試驗資料的小區(qū)平均值法，各參數(shù)也存在偏差，但只有表5～6中機誤、表7中機誤和區(qū)組、表8中機誤和區(qū)組兩因子的偏差達到顯著水平，因此其估計值是有偏差的。進一步觀察表5～8可見,參試因子的均方誤差中小區(qū)平均值法都大于轉(zhuǎn)化分析法，因此小區(qū)平均值法的參試因子方差分量的精度要比轉(zhuǎn)化分析法小。

采用蒙特卡羅模擬法產(chǎn)生的非平衡試驗數(shù)據(jù)，按照轉(zhuǎn)化分析法和小區(qū)平均值法進行分析，其結(jié)果（列于表4-8）的變化趨勢完全一致：試驗資料的轉(zhuǎn)化分析法，各參數(shù)都存在一定的偏差，但在這些偏差中，只有個別試驗中，區(qū)組的方差分量偏差達到顯著水平；而小區(qū)平均值法，通常區(qū)組和誤差因子的偏性分析結(jié)果均顯著。家系方差Vf 的均方誤兩者比較接近，對于Ve和Vb，轉(zhuǎn)化分析法的均方誤通常顯著小于小區(qū)平均值法的分析結(jié)果。

比較表4-8中參數(shù)均方誤差的相對大小，除家系遺傳方差兩者十分接近外，其它誤差因子和區(qū)組因子的均方誤差，轉(zhuǎn)化分析法的結(jié)果均明顯小于小區(qū)平均值法的均方誤差，而均方誤差的大小則表明該法參數(shù)估計效益的高低和分析結(jié)果精確性的高低，估計這種結(jié)果與區(qū)組數(shù)太少有關(guān)[2-6]。

（3）小區(qū)平均值法獲得的遺傳參數(shù)評估

試驗林在83%存活率的條件下，根據(jù)Ve、Vb、Vf的大小和調(diào)節(jié)系數(shù)的大小，可計算出理論上遺傳力的大小。以此為參照物，對非平衡試驗數(shù)據(jù)轉(zhuǎn)化前后的實際遺傳力進行評價，見表9。

從表9可見：（1）所有的試驗，轉(zhuǎn)化分析法的統(tǒng)計效率要比小區(qū)平均值法的的統(tǒng)計效率要高，利于逆向選擇和前向選擇；（2）轉(zhuǎn)化分析法的分析精確性要高于小區(qū)平均值法，因為前者的（Mse/bn）^.5比后者要??；（3） 小區(qū)平均值法的家系遺傳力比轉(zhuǎn)化分析法的結(jié)果低，并且遺傳力的相對誤差在-0.36%到-6.91%間變化，不利于逆向選擇。

表 9 轉(zhuǎn)化分析法與小區(qū)平均值法遺傳力的大小和精度比較Tab. 9 Comparison of heritability and accuracy between transformation analysis method and plot average method

綜合以上研究，從參數(shù)的精度和統(tǒng)計效益上證明了林木中小區(qū)平均值法不及轉(zhuǎn)化分析法的分析結(jié)果，在林木子代測定資料的分析中宜優(yōu)先采用轉(zhuǎn)化分析法。

2.2 固定模型條件下,小區(qū)平均值法的評估結(jié)果

在固定模型條件下，同一試驗資料，以轉(zhuǎn)化分析法獲得的家系效應(yīng)值為參照對象，分析了小區(qū)平均值法獲得的家系效應(yīng)值的大小，相對誤差，秩變化及20%的入選率時，優(yōu)良家系選擇的失誤率大小如下:

以轉(zhuǎn)化分析法為參照對象，小區(qū)平均值法的家系效應(yīng)值，將產(chǎn)生1.02%～7.08%的相對誤差；小區(qū)平均數(shù)值法中的許多家系秩與這些家系轉(zhuǎn)化分析獲得的秩不一致；采取20%的入選強度時，逆向選擇會產(chǎn)生1/8—2/8的失誤概率。這一結(jié)果與以前的研究結(jié)果[13]相一致。

3 討論與結(jié)論

林木遺傳育種中，規(guī)則的、平衡的試驗設(shè)計常常獲得的是非平衡試驗資料。由于林木遺傳育種是以單株值作為操作程序，此時采用轉(zhuǎn)化分析法具有優(yōu)越性，小區(qū)平均值法更適合處理平衡規(guī)則的林木無性系田間試驗數(shù)據(jù)，存在缺區(qū)時，還要用最小二乘法補平，這非常麻煩.據(jù)齊明（2009）利用9年生的杉木全同胞子代試驗林資料（實生起源）的研究表明：當出現(xiàn)數(shù)據(jù)缺失時，轉(zhuǎn)化分析法（即直接采用缺株分析），其分析結(jié)果造成的統(tǒng)計誤差，比該資料用小區(qū)平均值法（用最小二乘法補齊數(shù)據(jù)）獲得的分析結(jié)果的誤差小[13]。

蒙特卡羅（Monte Carlo）模擬法，產(chǎn)生的非平衡試驗數(shù)據(jù)（針對實生苗），其分析結(jié)果在4、5、6、8、10株小區(qū)試驗中表明：（1）小區(qū)平均值法適合于平衡規(guī)則的數(shù)據(jù)，在處理林木試驗非平衡資料時，尤其是存在缺區(qū)時，要用最小二乘法補平，這與轉(zhuǎn)化分析法相比，更加麻煩；（2）相對于轉(zhuǎn)化分析法，小區(qū)平均值法的計算量小，但沒法獲得單株的遺傳變異性，即統(tǒng)計效率低，不利于前向選擇；（3）由于小區(qū)平均值法獲得的家系遺傳力，其大小和精度要低于轉(zhuǎn)化分析法，不利于逆向選擇；（4）轉(zhuǎn)化分析法在若干參數(shù)上的偏性大小和均方誤大小，明顯優(yōu)于小區(qū)平均值法的分析結(jié)果，這說明林木中采用小區(qū)平均值法處理非平衡資料時，其參數(shù)的精確性不及轉(zhuǎn)化分析法。

在處理林木正交試驗非平衡數(shù)據(jù)時，由于小區(qū)平均值法的統(tǒng)計效率低，獲得的參數(shù)精確性差，不利于逆向選擇和前向選擇，建議采用轉(zhuǎn)化分析法來處理。