運用SAS廣義線性混合模型分析裂區(qū)試驗數(shù)據(jù)

張久權(quán), 李彩斌, 凌愛芬, 董建新

(1.中國農(nóng)業(yè)科學(xué)院煙草研究所, 農(nóng)業(yè)農(nóng)村部煙草生物學(xué)與加工重點實驗室, 山東 青島 266101; 2.貴州省煙草公司畢節(jié)市公司, 貴州 畢節(jié) 551700; 3.四川省煙草公司涼山州公司, 四川 西昌 615000)

裂區(qū)設(shè)計由于能夠突破試驗條件限制，為實際操作提供方便等諸多優(yōu)點[1]，受到廣大科研人員的普遍歡迎，在許多研究領(lǐng)域得到廣泛運用[2-11]。然而，由于此類試驗數(shù)據(jù)的統(tǒng)計分析較難，F(xiàn)檢驗的主區(qū)因子和副區(qū)因子所用的誤差項不同[4,6]，主區(qū)因子間、同一主區(qū)因子不同副區(qū)因子間均值比較的t值計算方法也不同[1,12]，使裂區(qū)試驗的數(shù)據(jù)分析具有很大挑戰(zhàn)性。如果僅依賴計算機軟件程序進行統(tǒng)計分析，一般很難處理，有些需要手工計算，但計算量很大[1]。從國內(nèi)外研究進展看，人們對裂區(qū)設(shè)計的統(tǒng)計分析常常存在誤用的現(xiàn)象[2,13]。Yossa等[14]調(diào)查表明，65%的已發(fā)表SCI論文存在試驗設(shè)計、F檢驗或多重比較錯誤。其中，裂區(qū)設(shè)計占比較大，但目前有關(guān)裂區(qū)設(shè)計試驗的統(tǒng)計方法存在不系統(tǒng)或太復(fù)雜等問題，通常需要手工計算，給科技工作者增加了負擔(dān)。

傳統(tǒng)上，學(xué)者們在利用SAS程序進行裂區(qū)設(shè)計的方差分析和均值比較時，主要采用一般線性模型(general linear model，GLM)程序模塊，而該程序模塊是在20世紀80年代開發(fā)的，當(dāng)時沒有考慮到在進行裂區(qū)設(shè)計的F檢驗時，主區(qū)因子與副區(qū)因子需采用不同的誤差均方與自由度，均值比較也需采用不同的誤差均方和自由度，因而計算F值時需要編寫大量的復(fù)雜代碼來進行修正[15]。后來，雖然可以用 “h=……”語句進行彌補，但是進行處理間均值比較時，需要選用合適的、不同的誤差項。然而，對于某些情況，GLM程序模塊沒有可用于定義所需誤差項的選項，所以誤差項無法列出[15]。

近年來，SAS等商業(yè)軟件進行了創(chuàng)新，廣義線性混合模型(general linear mixed model，GLIMMIX)有了很大發(fā)展[15]，其能適用于多種分布的變量[16]。本研究應(yīng)用SAS的GLIMMIX程序模塊對裂區(qū)試驗數(shù)據(jù)進行方差分析和均值比較，通過與傳統(tǒng)的GLM程序模塊和混合模型(MIXED)程序模塊進行對比，介紹GLIMMIX 程序模塊在裂區(qū)試驗數(shù)據(jù)統(tǒng)計中的優(yōu)勢，以供廣大科技工作者參考。

1 材料與方法

1.1 試驗材料與設(shè)計

裂區(qū)試驗的處理設(shè)計一般采用全因子試驗，即各因子的水平全部交叉，構(gòu)成完整的處理組合。關(guān)于主、副區(qū)因子的確定與分配，參照張久權(quán)等[1]方法。下面將采用其中實例進行試驗設(shè)計和原始數(shù)據(jù)分析，實例設(shè)計如下。

第一年在福建省三明市進行烤煙品種比較試驗，單因子4水平隨機區(qū)組設(shè)計，4水平為4個烤煙品種(翠碧1號、KRK26、YN116、CF209)，重復(fù)3次。第二年，為了探明移栽期對各品種產(chǎn)量的影響，對4個品種增加了3個移栽期:1月24日(1)、1月31日(2)、2月8日(3)處理，試驗改為裂區(qū)設(shè)計，品種為主區(qū)因子，移栽期為副區(qū)因子。第二年的產(chǎn)量數(shù)據(jù)見表1。

表1 產(chǎn)量數(shù)據(jù)Excel輸入表Table 1 Yield data entry for Excel

1.2 方差分析模型

對于兩因素裂區(qū)設(shè)計，如果主區(qū)因子的排列為隨機區(qū)組，副區(qū)因子在主區(qū)內(nèi)隨機排列，其線性可加效應(yīng)模型如下。

Yijk=μ+αi+βj+αiβj+rk+wik+eijk

(1)

式中，Yijk為某小區(qū)的觀察值，如小區(qū)產(chǎn)量；μ為總體均數(shù)；αi為主區(qū)因素A第i水平的效應(yīng)，βj為副區(qū)因素B第j水平的效應(yīng)，αiβj為A因素i水平與B因素j水平交互作用效應(yīng)，均為固定效應(yīng)；rk為區(qū)組K的效應(yīng)，wik為主區(qū)誤差，eijk為副區(qū)誤差，均為隨機效應(yīng)。為了便于理解，將主區(qū)元素放在前，副區(qū)元素放在后，將(1)式改寫為(2)式。

Yijk=μ+rk+αi+wik+βj+αiβj+eijk

(2)

式中，區(qū)組效應(yīng)rk，主區(qū)誤差wik和副區(qū)誤差eijk應(yīng)滿足以下條件。

(1)所有的隨機效應(yīng)項rk、wik和eijk相互獨立。

(2)區(qū)組效應(yīng)rk為iidN(0, σR2)，即rk服從同一正態(tài)分布，相互獨立，方差均為σR2。

(3)主區(qū)誤差wik為iidN(0, σw2)。

(4)副區(qū)誤差eijk為iidN(0, σ2)。

從模型可看出，對于這類資料的分析可以分兩步。第一步，將總變異分解為主區(qū)組間變異、主區(qū)因子間變異和主區(qū)誤差；第二步，將整個試驗的總變異分解為主區(qū)變異、副區(qū)因子變異、主副因子之間的交互作用和副區(qū)誤差。

如果主區(qū)因子為完全隨機，需要將公式(1)和(2)中的區(qū)組效應(yīng)rk刪除。其他試驗設(shè)計的線性模型依此類推。

1.3 方差分析和多重比較SAS程序

1.3.1用GLIMMIX模塊進行F檢驗和所有處理組合間的均值比較 運用SAS 9.4軟件的GLIMMIX模塊對表1數(shù)據(jù)進行F檢驗。SAS程序1的代碼如下。

ODS RTF;

ODS graphics on;

PROC IMPORT

OUT= split

DATAFILE= "E:yancao.xls"

DBMS=EXCEL REPLACE;

SHEET="Sheet1$";

GETNAMES=YES;

RUN;

PROC GLIMMIX data = split;

Class block variety date;

Model yield = variety date variety*date / ddfm=satterth;

Random block block*variety;

Lsmeans variety*date/plot=meanplot (sliceby=variety join);

Lsmeans variety date variety*date /diff;

Lsmestimate variety*date ′date1 vs date2 in variety1 ′ 1 -10;

RUN;

ODS graphics off;

ODS RTF close;

程序說明：輸入代碼時注意分號為英文字符。數(shù)據(jù)文件存放地址為E:yancao.xls sheet1，數(shù)據(jù)格式見表1。PROC IMPORT直接讀取EXCEL文件的數(shù)據(jù)。其他數(shù)據(jù)格式如文本格式也可以輸入，但EXCEL比較方便。PROC GLIMMIX調(diào)用GLIMMIX程序，采用限定的最大似然法進行F檢驗和多重比較。Class 語句列出所有的分類變量。Model語句“=”右邊只需列出固定效應(yīng)變量。ddfm=satterth對誤差方差、自由度等進行校正，以保證F值的正確性?！癓smeans variety*date/plot = meanplot (sliceby=variety join)”語句輸出交互作用效果圖，join選項將各點連成曲線?！癓smeans variety*date /diff” 語句對所有處理水平組合的最小二乘均值(least square means，lsmean)進行兩兩比較，計算t值和差異顯著性。在平衡的試驗設(shè)計中，最小二乘均值與處理均值相同。

1.3.2用GLIMMIX模塊進行F檢驗和部分處理組合間的均值比較 SAS程序1中的語句“Lsmeans variety date variety*date /diff;”輸出所有處理組合的最小二乘均值及其差值，本例中共64對，輸出超過30頁。為了有選擇性地比較處理組合間的均值差異，需要將上述語句進行替換。如要比較date1(1月24日)與date2(1月24日)的產(chǎn)量差異，即需要替換為“Lsmestimate variety*date ′date1 vs date2 in variety1′ 1 -10;”；要以variety1 date1(翠碧1號，1月24日移栽)為對照進行產(chǎn)量差異比較，則需替換為“Lsmeans variety*date /diff = control (′1′ ′1′);”；要比較品種1(翠碧1號)與所有其他3個品種平均產(chǎn)量的差異，則需替換為“Estimate ‘variety1 vs others’ variety -3 1 1 1 /divisor =3;”。

1.3.3用MIXED程序模塊進行F檢驗和均值比較 SAS的MIXED也可用來進行裂區(qū)設(shè)計的方差分析。進行F檢驗時，class、model、random、lsmeans語句的用法與GLIMMIX的完全相同，F(xiàn)檢驗和均值比較的程序如下。

PROC MIXED data = split;

Class block variety date;

Model yield = variety date variety*date/ddfm=satterth;

Random block block*variety;

Lsmeans variety date variety*date /diff;RUN;

2 結(jié)果與分析

2.1 運用GLIMMIX模塊進行F檢驗、多重比較及交互作用比較

運行程序1，得 到F檢驗的主要結(jié)果見表2?？梢钥闯?，烤煙產(chǎn)量主區(qū)因子品種間、副區(qū)因子移栽期之間，以及他們的交互作用效果都達到極顯著水平(P<0.001)。因此，有必要分別對各處理組合的產(chǎn)量均值差異進行比較。受篇幅所限，僅列出部分處理組合間的產(chǎn)量均值差異(表3)。其中，第1行數(shù)據(jù)表示variety1(翠碧1號)，在date1(1月24日)移栽的烤煙產(chǎn)量比在date2(1月31日)移栽的極顯著地高出780 kg·km-2(P<0.001)。

表2 F檢驗結(jié)果(Ⅲ型)Table 2 Output of F test (Type Ⅲ SS)

表3 部分處理組合的產(chǎn)量均值差異比較Table 3 Difference comparison of mean yield between part treatment combinations

品種和移栽期的交互作用效果圖見圖1，可以看出，各品種烤煙產(chǎn)量隨不同移栽期處理的變化趨勢明顯不同，variety1、2隨移栽期的推遲，產(chǎn)量下降；variety3、4隨移栽期的推遲，產(chǎn)量增加。在date1(1月24日)移栽時，variety1產(chǎn)量最高，variety4最低；而在date3(2月8日)移栽時，品種1產(chǎn)量最低，品種3最高。如果品種和移栽期處理間的交互作用不存在，這些線將呈現(xiàn)平行趨勢。

綜上，結(jié)合均值比較結(jié)果(圖1)能夠很容易地確定產(chǎn)量最高的最佳處理組合，即最佳移栽期為date1(1月24日)的variety1(翠碧1號)。從表3的第1、2行可以看出，variety1的date1組合與variety1的date2組合，及variety1的date1組合與variety1的date3組合差異均為極顯著水平。因此，variety1(翠碧1號)的最佳移栽期為date1(1月24日)。

圖1 品種與移栽期交互作用Fig.1 Profile of interaction between variety and date

為了進行方法說明，現(xiàn)假定品種與移栽期的交互作用效果不顯著，即假定表2中的variety*datePr>F大于0.05，即品種的產(chǎn)量效果不受移栽期的影響，同時移栽期的產(chǎn)量效果與品種無關(guān)，此時的均值比較僅需對品種和移栽期分別進行。用Lsmeans variety/diff比較品種間的均值差異，結(jié)果見表4?？梢?，綜合3個移栽期的平均產(chǎn)量結(jié)果，variety1、2、3與4的產(chǎn)量差異達到極顯著水平，而其他品種間無顯著差異。用Lsmeans date/diff 比較不同移栽期產(chǎn)量之間的均值差異，結(jié)果見表5?？梢?，綜合4個品種的平均結(jié)果，date 1與date2和3的產(chǎn)量差異均達到極顯著水平，而date2與date 3間無顯著差異。需要特別說明的是，這種均值比較只有在品種與移栽期交互作用效果不顯著的前提下才能進行，否則沒有意義。

表4 各品種間產(chǎn)量差異顯著性(所有移栽期平均)比較示例Table 4 Example of mean difference between varieties across all of the dates

表5 各移栽期間平均產(chǎn)量的差異顯著性示例Table 5 Example of mean difference of yield between dates across all varieties

2.2 部分處理組合間的均值比較結(jié)果

為了節(jié)省輸出文檔的篇幅，避免產(chǎn)生大量不必要的輸出，可以進行某些特定均值間的差異比較，如1.3.2 中語句“Lsmestimate variety*date ′date1 vs date2 in variety1 ′ 1 -10;”的輸出結(jié)果見表6，可以看出， date1(1月24日)與date2(1月31日)的產(chǎn)量差值為780.0 kg·hm-2，與表3中date1(1月24日)與date2(1月31日)的產(chǎn)量差值相同。說明語句“Lsmestimate variety*date ′date1 vs date2 in variety1 ′ 1 -10;”的作用與“Lsmestimate variety*day”是相同的，但該語句的輸出內(nèi)容更簡潔。

表6 Lsmestimate 輸出示例Table 6 Example output of Lsmestimate

3 討論

在進行裂區(qū)設(shè)計的方差分析時，通常采用ANOVA或GLM程序模塊[1, 5, 17-18]。由于主、副區(qū)處理需要采用不同的誤差方和及自由度，ANOVA或GLM缺乏內(nèi)置程序進行自動處理，某些情況下需要手工計算，大大增加了計算量和難度[1]，且容易出錯。本研究采用GLIMMIX程序模塊能恰當(dāng)?shù)亟鉀Q這些問題，即使是非平衡設(shè)計(如各處理組合的重復(fù)次數(shù)不等或缺區(qū))的情況。

本研究采用的實例數(shù)據(jù)與孫久權(quán)等[1]的實例數(shù)據(jù)完全相同，而孫久權(quán)等[1]采用GLM進行統(tǒng)計分析，程序繁雜，如其在計算variety處理效應(yīng)(主區(qū)因子效應(yīng))的F值時，SAS的默認方法有誤，需要通過語句"Test h = replicate variety e = replicate*variety”進行修正，才能得到更準確的、與采用GLIMMIX程序模塊相同的F值(64.92)。MIXED和GLIMMIX能夠區(qū)別對待固定效應(yīng)和隨機效應(yīng)變量，在Model語句中，只需列出固定效應(yīng)變量即可，如本研究中可以采用以下語句。

Model yield = variety date variety*date / ddfm = satterth;

Random block block*variety;

而GLM程序模塊的model語句需要列出所有變量，如孫久權(quán)等[1]的model語句如下(replicate與block等同)。

Model yield = replicate variety date replicate*variety variety*date;

Test h = replicate variety e = replicate*variety;

為了正確計算主區(qū)因子的F值，GLM程序模塊需要使用Test語句。類似地，使用LSMeans進行主區(qū)因子間的差異比較時，必須單獨指明誤差項，即“LSMeans variety/ e = replicate(variety) stderr, tdiff;”；否則GLM程序模塊將默認使用副區(qū)誤差項計算t值。

在MIXED和GLIMMIX程序模塊中，只要在random語句中列出隨機效應(yīng)項，如本研究中的“Random block block*variety”，SAS在計算F和t值等統(tǒng)計量時，就會根據(jù)適當(dāng)?shù)姆讲罱Y(jié)構(gòu)，自動選用恰當(dāng)?shù)恼`差項方和與自由度等，這是GLIMMIX程序本身事先內(nèi)置的功能。GLM缺乏該內(nèi)置功能，需要列出誤差項，即使使用random語句也不能實現(xiàn)。然而，有些情況下，GLM程序模塊缺乏此內(nèi)置誤差項，因而GLM不能自動計算該統(tǒng)計量[15,19-20]。此外，GLIMMIX程序模塊可以方便地輸出交互作用效果圖(圖1)，但GLM和MIXED無此功能。

雖然GLIMMIX模型在分析裂區(qū)設(shè)計試驗數(shù)據(jù)時，優(yōu)勢很多。但是，GLIMMIX模型是基于線性混合模型進行迭代的[20]，不能用于嵌套模型的比較，且允許殘差方差及不同協(xié)方差結(jié)構(gòu)的類型有限[21]，由于PROC GLIMIX對隨機效應(yīng)的標準誤估計有偏差，因而不能用于隨機效應(yīng)的顯著性檢驗。

在GLIMMIX進行裂區(qū)設(shè)計試驗數(shù)據(jù)分析時，仍有一些要點需要注意。本研究程序1用ddfm=satterth對標準誤、自由度等進行校正，以保證F值的正確性。Kenward和Roger[17]用“ddfm = kr”來修正標準誤、自由度、檢驗統(tǒng)計量等的計算，所得結(jié)果的準確性得到進一步提高。如果裂區(qū)設(shè)計為平衡設(shè)計，satterth和kr效果完全相同；如果為非平衡設(shè)計或存在缺區(qū)，尤其是在重復(fù)測量和空間相關(guān)的情況下，二者有本質(zhì)區(qū)別，推薦使用kr。在輸出交互作用圖時，本研究將變量date視為分類變量，所以x軸為等距。如果x軸要按比例顯示，需要采用Gplot程序作圖。如果要將variety設(shè)為橫坐標，可用sliceby= date 代替sliceby=variety 選項。程序1的輸出不含某些均值比較。例如，檢驗品種1與所有其他3個品種平均產(chǎn)量的差異顯著性時，可以通過estimate、contrast或lsmestimate語句后的系數(shù)來實現(xiàn)，語句為“Estimate ‘variety1 vs others’ variety -3 1 1 1 /divisor =3;”。其中，系數(shù)的順序一定要與程序1中class語句的變量順序一致，即variety在date之前。本研究中，各品種的產(chǎn)量與date存在直線或曲線關(guān)系(圖1)。因此，可以以移栽期為自變量，產(chǎn)量為因變量進行回歸分析，建立產(chǎn)量與移栽期的回歸方程，用以預(yù)測產(chǎn)量，具體可采用SAS的REG、MIXED等模塊進行擬合和統(tǒng)計檢驗。

在MIXED模塊中，SAS使用限制最大可能性估計(REML)或ANOVA方式進行方差分量，默認方式為REML。若采用ANOVA進行估計，語句為“PROC MIXED data = split method = type3”。對于平衡設(shè)計，REML與ANOVA結(jié)果相同，但ANOVA能夠輸出更詳細的方差分析結(jié)果。在進行均值比較時，GLIMMIX模塊可采用estimate、contrast或lsmestimate語句。但lsmestimate語句是針對GLIMMIX模塊開發(fā)的，MIXED模塊不能使用。然而，estimate、contrast語句的寫法較繁瑣[19]，SAS公司對其進行了簡化，變?yōu)閘smestimate，寫法更加簡便[20]。