裂區(qū)試驗(yàn)設(shè)計(jì)中統(tǒng)計(jì)分析的SAS線性混合模型研究

董中東 張 寧 趙 磊 任 妍 孫叢葦 陳 鋒

（河南農(nóng)業(yè)大學(xué)農(nóng)學(xué)院，河南鄭州 450046）

裂區(qū)試驗(yàn)設(shè)計(jì)是農(nóng)業(yè)試驗(yàn)中常用的一種試驗(yàn)設(shè)計(jì)方法[1-3]。SAS GLIMM功效分析表明，當(dāng)試驗(yàn)單元誤差相同時(shí)，裂區(qū)試驗(yàn)設(shè)計(jì)較不完全隨機(jī)區(qū)組和完全隨機(jī)區(qū)組有著更高的統(tǒng)計(jì)功效。裂區(qū)試驗(yàn)設(shè)計(jì)可以分為二因素裂區(qū)和三因素裂區(qū)，后者又有其中兩因素組合作主區(qū)或者作副區(qū)的不同裂區(qū)試驗(yàn)。但是，在裂區(qū)試驗(yàn)的設(shè)計(jì)和分析過程中經(jīng)常會(huì)存在一些問題，尤其是其統(tǒng)計(jì)分析。裂區(qū)試驗(yàn)的統(tǒng)計(jì)分析難點(diǎn)主要有2個(gè)方面：一方面，在進(jìn)行方差分析各因素主效應(yīng)和互作效應(yīng)的檢驗(yàn)所用到的誤差項(xiàng)不同；另一方面，最主要的難點(diǎn)是其互作效應(yīng)的多重比較，由于在對(duì)不同的處理組合進(jìn)行比較時(shí)用到的誤差方差及其自由度的計(jì)算非常繁瑣，因而裂區(qū)試驗(yàn)在統(tǒng)計(jì)分析時(shí)經(jīng)常會(huì)存在一些問題[4-5]。

SAS 軟件（SAS Institute，Cary，NC）是一款功能非常強(qiáng)大的統(tǒng)計(jì)分析軟件，尤其是近年來其開發(fā)了線性混合模型（linear mixed model）和廣義線性混合模型（generalized linear mixed model）的應(yīng)用程序 PROC MIXED和PROC GLIMMIX，這些程序均使用了很多近年來的統(tǒng)計(jì)理論和方法，使SAS對(duì)試驗(yàn)數(shù)據(jù)的分析有了很大提高。由于裂區(qū)試驗(yàn)的分析需要估計(jì)2種誤差方差以及互作項(xiàng)的多重比較分母誤差項(xiàng)的矯正等問題，在SAS公司出版的著作SAS for Mixed Models中不建議使用PROC GLM程序進(jìn)行裂區(qū)試驗(yàn)的分析。當(dāng)試驗(yàn)數(shù)據(jù)滿足方差分析的基本假定時(shí)，若涉及隨機(jī)效應(yīng)的估計(jì)，PROC MIXED基本不需要太多額外的程序編寫，用其標(biāo)準(zhǔn)程序語句就可以實(shí)現(xiàn)。PROC MIXED還可以實(shí)現(xiàn)三因素裂區(qū)、時(shí)間裂區(qū)、時(shí)空裂區(qū)、多年多點(diǎn)等試驗(yàn)的統(tǒng)計(jì)分析。這里僅以二因素裂區(qū)試驗(yàn)設(shè)計(jì)的統(tǒng)計(jì)分析介紹SAS mixed程序和分析結(jié)果。

1 裂區(qū)統(tǒng)計(jì)分析的混合線性模型程序及結(jié)果解釋

首先以一個(gè)裂區(qū)試驗(yàn)介紹一般情況下的SAS分析程序：有一水稻不同灌水方式（A）和種植密度（B）試驗(yàn)，主處理為3種灌水方式，分別以A1、A2和A3表示，副處理為3種密度，分別以B1、B2和B3表示，設(shè)置3個(gè)區(qū)組，其小區(qū)產(chǎn)量結(jié)果如表1，本例題和數(shù)據(jù)引自《農(nóng)業(yè)試驗(yàn)設(shè)計(jì)與統(tǒng)計(jì)分析》[6]。

表1 灌水方式和密度二因素裂區(qū)試驗(yàn)產(chǎn)量結(jié)果

分析程序如下：

程序解釋：首先用DATA步建立一個(gè)名為spd_1的數(shù)據(jù)集，數(shù)據(jù)集包含4個(gè)變數(shù)，其名稱分別為blk、A、B 和 yld，blk代表區(qū)組，分別用 1、2、3代表 3個(gè)不同區(qū)組；A代表A因素，A1、A2和A3代表A因素的3個(gè)水平，為字符型，需加$說明其為字符型；B代表B因素，B1、B2和B3代表B因素的3個(gè)水平，也為字符型，需加以說明；yld代表小區(qū)產(chǎn)量。@@符號(hào)表明數(shù)據(jù)連續(xù)讀取，datalines下面的數(shù)據(jù)為要讀取的數(shù)據(jù)，到分號(hào)代表數(shù)據(jù)輸入結(jié)束。

接下來為PROC MIXED程序步，首先PROC MIXED說明調(diào)用的程序名稱，nobound選項(xiàng)的作用是在方差組分的估計(jì)時(shí)，避免SAS把估計(jì)的負(fù)的方差組分當(dāng)作0來輸出和進(jìn)一步計(jì)算，此選項(xiàng)在主區(qū)誤差組分為負(fù)值時(shí)尤其重要，它會(huì)影響到主處理的檢驗(yàn)結(jié)果。model語句等號(hào)左邊yld為要分析的試驗(yàn)指標(biāo)或性狀，右邊為固定效應(yīng)項(xiàng)，注意這里和常用的PROC GLM用法不同，mixed程序中等號(hào)右邊只包含固定效應(yīng)，而不能有隨機(jī)效應(yīng)，A|B，為A、B和A*B這3項(xiàng)的簡寫。后面的選項(xiàng)DDFM，是進(jìn)行計(jì)算分母自由度的矯正方法，有6種方法可供選擇，現(xiàn)在一般推薦使用KENWARDROGER（可簡寫為KENROG或 KR）或 SATTERTHWAITE（可簡寫為 SATTERTH或SAT）方法。random語句是指定隨機(jī)效應(yīng)，其對(duì)于裂區(qū)的統(tǒng)計(jì)分析尤為重要，只有把隨機(jī)效應(yīng)項(xiàng)設(shè)置正確才能得到正確的分析結(jié)果。這里把區(qū)組作為隨機(jī)效應(yīng)進(jìn)行分析，blk*A則是定義主區(qū)誤差項(xiàng)的隨機(jī)效應(yīng)。lsmeans語句定義輸出A、B和A*B項(xiàng)的最小二乘法估計(jì)的均值，后面的選項(xiàng)diff是進(jìn)行A、B和A*B的各個(gè)水平間最小二乘法估計(jì)均值差數(shù)的t檢驗(yàn)（LSD 法），adjust=tukey，輸出 tukey法多重比較的均值差數(shù)的結(jié)果。輸出結(jié)果如表2、表3、表4所示。

表2 各隨機(jī)效應(yīng)項(xiàng)的方差估計(jì)

表3 固定效應(yīng)的檢驗(yàn)結(jié)果

表4 各處理和處理組合均值差數(shù)顯著性部分結(jié)果

表2為各隨機(jī)效應(yīng)項(xiàng)的方差分量估計(jì)值，由于在程序中有nobound選項(xiàng)，所以區(qū)組項(xiàng)的方差為負(fù)值，如果沒有該選項(xiàng)，其值將會(huì)為0。另外2項(xiàng)分別為主區(qū)誤差組分和副區(qū)誤差組分。

表3是各固定效應(yīng)的F檢驗(yàn)結(jié)果，結(jié)果表明A因素各水平間差異不顯著（P=0.471 8），B因素的水平間差異極顯著（P=0.004 0），2個(gè)因素存在極顯著的交互作用（P=0.000 8）。由于互作項(xiàng)存在極顯著差異，主要對(duì)互作項(xiàng)的結(jié)果進(jìn)行解釋和分析。表4為各因素不同水平和處理組合差數(shù)的比較結(jié)果，對(duì)于各因素水平間的差異顯著性結(jié)果，只是看其計(jì)算結(jié)果是否正確，本例題應(yīng)主要看各處理組合的差異顯著性分析。先看同主區(qū)不同副區(qū)的一個(gè)比較結(jié)果，以LSD法的A1B1和A1B2為例，該比較只與副區(qū)誤差有關(guān)，其自由度為副區(qū)誤差自由度12。再看不同主區(qū)的比較，以A1B1和A2B2為例，由于不同主區(qū)比較的誤差是主區(qū)誤差和副區(qū)誤差的合成，所以其既不是主區(qū)誤差自由度，也不是副區(qū)誤差自由度，需要根據(jù)公式矯正，其矯正公式如下：

式中，b為副處理水平，MSEa為主區(qū)誤差均方，MSEb為副區(qū)誤差均方，dfEa為主區(qū)誤差自由度，dfEb為副區(qū)誤差自由度，矯正后的自由度為12.3。從分析結(jié)果可以看出，對(duì)于互作項(xiàng)的多重比較，程序會(huì)自動(dòng)根據(jù)要比較的2個(gè)處理組合算出其所需標(biāo)準(zhǔn)誤和自由度，這種比較使用PROC GLM不能實(shí)現(xiàn)。

通過使用PROC MIXED對(duì)二因素裂區(qū)試驗(yàn)結(jié)果的分析可以看出，分析程序非常簡潔，而且可以得到所要分析的結(jié)果。但是，由于處理組合的比較非常多，需要把結(jié)果進(jìn)一步整理。當(dāng)試驗(yàn)需要對(duì)一些綜合的效應(yīng)進(jìn)行比較時(shí)，比如A1水平下B1+B2的平均效應(yīng)與A1水平下B3的效應(yīng)進(jìn)行比較時(shí)需要使用estimate或者contrast語句才能分析出結(jié)果，本文不再介紹。

2 裂區(qū)分析的一般線性模型程序及結(jié)果

由于SAS mixed程序輸出的方差分析表中沒有各效應(yīng)項(xiàng)的SS和MS，為便于驗(yàn)證mixed程序結(jié)果的正確性，這里給出SAS GLM程序和輸出的方差分析結(jié)果。

GLM程序如下：

方差分析結(jié)果如表5、表6所示。從表5和表6可以看出，A因素的F值（0.91）和對(duì)應(yīng)概率（0.471 8）、B 因素的 F 值（9.06）和對(duì)應(yīng)概率（0.004 0）、A×B 互作的F值（10.22）和對(duì)應(yīng)概率（0.000 8）與mixed程序輸出結(jié)果完全一致。處理組合的多重比較可以在程序中加入語句“l(fā)smeans A*B/diff；”進(jìn)行，在結(jié)果中有需要誤差和自由度合成的結(jié)果都是不正確的。通過GLM和mixed程序的運(yùn)行結(jié)果可以看出，在裂區(qū)分析中mixed程序較glm程序更加準(zhǔn)確、方便和快捷。

表5 方差分析結(jié)果

表6 使用Ⅲ型MS作為blk*A的誤差項(xiàng)的假設(shè)檢驗(yàn)