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

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

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

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

1 模型與方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表 1 單因素隨機(jī)區(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 半同胞試驗(yàn)小區(qū)平均值法的期望方差結(jié)構(gòu)Tab. 2 Expected variance structure of plot average method in half-sib experimental plot

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

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

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

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

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

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

所有計(jì)算分析均在MATLAB7.0[19-20]和Excel 2003平臺(tái)上完成。

2 結(jié)果與分析

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

（1）五個(gè)模擬試驗(yàn)中兩種方法獲得參試因子負(fù)方差分量的頻率

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

從表9可見：（1）所有的試驗(yàn)，轉(zhuǎn)化分析法的統(tǒng)計(jì)效率要比小區(qū)平均值法的的統(tǒng)計(jì)效率要高，利于逆向選擇和前向選擇；（2）轉(zhuǎn)化分析法的分析精確性要高于小區(qū)平均值法，因?yàn)榍罢叩模∕se/bn）^.5比后者要??；（3） 小區(qū)平均值法的家系遺傳力比轉(zhuǎn)化分析法的結(jié)果低，并且遺傳力的相對(duì)誤差在-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)計(jì)效益上證明了林木中小區(qū)平均值法不及轉(zhuǎn)化分析法的分析結(jié)果，在林木子代測(cè)定資料的分析中宜優(yōu)先采用轉(zhuǎn)化分析法。

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

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

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

3 討論與結(jié)論

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

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

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