裂區試驗設計中統計分析的SAS線性混合模型研究

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

（河南農業大學農學院，河南鄭州 450046）

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

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

1 裂區統計分析的混合線性模型程序及結果解釋

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

表1 灌水方式和密度二因素裂區試驗產量結果

分析程序如下：

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

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

表2 各隨機效應項的方差估計

表3 固定效應的檢驗結果

表4 各處理和處理組合均值差數顯著性部分結果

表2為各隨機效應項的方差分量估計值，由于在程序中有nobound選項，所以區組項的方差為負值，如果沒有該選項，其值將會為0。另外2項分別為主區誤差組分和副區誤差組分。

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

式中，b為副處理水平，MSEa為主區誤差均方，MSEb為副區誤差均方，dfEa為主區誤差自由度，dfEb為副區誤差自由度，矯正后的自由度為12.3。從分析結果可以看出，對于互作項的多重比較，程序會自動根據要比較的2個處理組合算出其所需標準誤和自由度，這種比較使用PROC GLM不能實現。

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

2 裂區分析的一般線性模型程序及結果

由于SAS mixed程序輸出的方差分析表中沒有各效應項的SS和MS，為便于驗證mixed程序結果的正確性，這里給出SAS GLM程序和輸出的方差分析結果。

GLM程序如下：

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

表5 方差分析結果

表6 使用Ⅲ型MS作為blk*A的誤差項的假設檢驗