運用廣義線性混合模型分析隨機區(qū)組重復(fù)測量的試驗資料

張久權(quán) 閆慧峰 褚繼登 李彩斌

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

農(nóng)業(yè)科研常常需要進行長期定位試驗和多點試驗, 需要對同一對象(小區(qū)、某株作物等)進行重復(fù)測量, 例如, 為了探明烤煙移栽后的生長速度, 需要定點定株對烤煙的株高進行連續(xù)測量。這種對同一受試對象(subject)進行多次測量的試驗即為重復(fù)測量試驗[1-3]。在農(nóng)業(yè)試驗中, 主要包括3種類型的重復(fù)測量: (1)不同時間點對同一對象測量; (2)對同一對象不同部位的測量, 如烤煙上、中、下3個部位;(3)不同地點的測量, 如作物品種區(qū)試。時間、部位、地點均可以視為重復(fù)測量因子。由于個體間的差異,觀測值開始就高的個體, 后面也高, 反之一樣, 例如, 定點測量煙株的高度, 第 1次測量較高的煙株,后面測量時也會高。因此, 同一對象不同觀察值間往往存在時間或空間上的自相關(guān)性, 觀測點間隔越近, 往往相關(guān)性越大[3-6], 自相關(guān)是重復(fù)測量數(shù)據(jù)最大的特點。由于存在自相關(guān)性, 數(shù)據(jù)間相互獨立的基本要求不能滿足, 不能采用常規(guī)的統(tǒng)計方法進行統(tǒng)計分析[4-5,7-8]。

重復(fù)測量試驗在醫(yī)學領(lǐng)域十分普遍, 在農(nóng)業(yè)研究領(lǐng)域也越來越受到人們的重視, 特別是長期定位試驗。過去, 主要采用裂區(qū)設(shè)計或多變量統(tǒng)計方法分析這類數(shù)據(jù)[2,6,9], 不僅忽視了分析的前提條件[10],可能得出錯誤的結(jié)論, 也大大浪費了重復(fù)測量因素的時間或空間信息, 降低了統(tǒng)計分析的效能, 導致高層次、高效率的試驗得出低質(zhì)量、低水平的研究結(jié)論, 造成不必要的浪費與損失[4,11-12]。國外期刊從2004年開始就已經(jīng)對重復(fù)測量數(shù)據(jù)的分析進行了嚴格要求[1]。目前, 我國農(nóng)業(yè)領(lǐng)域, 有關(guān)對重復(fù)測量試驗數(shù)據(jù)進行科學統(tǒng)計分析方法的報道鮮見。

近年來, SAS等統(tǒng)計軟件進行了創(chuàng)新和發(fā)展。運用混合模型(mixed)和廣義線性混合模型(Generalized Linear Mixed Models, GLIMMIX)能很好地進行重復(fù)測量的數(shù)據(jù)分析[13]。重復(fù)測量試驗可以與許多試驗設(shè)計結(jié)合, 如完全隨機、隨機區(qū)組、裂區(qū)、拉丁方等等。隨機區(qū)組試驗是農(nóng)業(yè)試驗中普通采用的試驗設(shè)計, 尤其是田間試驗。因此, 本文以兩因素隨機區(qū)組試驗數(shù)據(jù)為例, 系統(tǒng)介紹應(yīng)用GLIMMIX進行重復(fù)測量數(shù)據(jù)方差分析和均值比較的方法, 并闡明傳統(tǒng)統(tǒng)計方法的不足, 供廣大農(nóng)業(yè)科研工作者參考。

1 材料與方法

1.1 試驗設(shè)計

1.1.1 土柱試驗 該試驗為室內(nèi)模擬試驗二因子(降水強度、N肥水平)三重復(fù)全因子設(shè)計, 小區(qū)排列方式為隨機區(qū)組。土壤為中國農(nóng)業(yè)科學院西南試驗基地(四川西昌市)收集的紫色土。降水強度分別為0.44、0.68、1.20 mm min-1, 降水量均為 50 mm 次-1。設(shè) N1 (常規(guī)施肥, CK)、N2 (減量 30%)、N3 (減量30%+生物炭)共3個氮肥水平。對土柱進行了12次模擬降雨, 第 6次開始收集到下滲淋洗液, 一直到第12次降雨完成, 每個土柱收集到7次淋洗液, 時間間隔相等。采用堿性過硫酸鉀消解紫外分光光度法測定淋洗液中的全氮[14]。7次淋洗被當作重復(fù)測量因素。采用SAS 9.4的GLIMMIX進行統(tǒng)計分析,數(shù)據(jù)輸入示例見表1。

表1 數(shù)據(jù)Microsoft Excel輸入表示例Table 1 Data entry example for Microsoft Excel

1.1.2 田間試驗 2018年4月開始, 在貴州畢節(jié)威寧縣黑石鎮(zhèn)開展煙田生物炭長期定位試驗, 采用單因子(生物炭用量)三水平(0、15、40 t hm-2)三重復(fù)隨機區(qū)組設(shè)計, 小區(qū)面積 67 m2, 供試作物為烤煙,分別于2018年的6月13日、7月17日、9月20日,2019年的 6月 11日、7月 24日、9月 26日采集0～20 cm 耕層土壤測定土壤速效鉀含量。為了進行說明, 設(shè)置2018年7月17日取樣的處理15 t hm-2重復(fù)1數(shù)據(jù)缺失, 按缺區(qū)處理。

1.2 分析比較

圖1顯示了土柱試驗3種不同施肥方式下總N淋失量隨淋次數(shù)的變化, 淋洗次數(shù)為重復(fù)測量的時間因素。我們可能需要確定如下處理組合間的差異是否顯著: (1)特定時間兩處理水平間的差異,如第6次淋洗時, 全氮淋失量 N1與N2處理的差異; (2)特定處理隨時間的變化, 如 N2處理, 全氮淋失量第6次與第8次淋洗的差異; (3)兩處理水平在各特定時間點的差異, 以便確定最優(yōu)處理組合,如N1第7次淋洗與N2第9次淋洗全氮淋失量的差異顯著性; (4)隨時間的變化(所有處理平均), 如3種氮處理平均, 全氮淋失量第 8次與第 10次淋洗的差異; (5)處理間差異(所有時間點平均), 如N1與 N2處理在所有時間點平均, 全氮淋失量的差異; (6)其他差異比較。

田間試驗與土柱試驗類似, 我們可以比較 2018年6月13日取樣時, 不同生物炭用量間土壤速效鉀含量的差異; 施用生物炭15 t hm-2后, 2018年6月與2019年6月土壤速效鉀含量的差異等。

1.3 方差分析模型

對于三因素(2個處理因素和1個時間因素)隨機區(qū)組設(shè)計(土柱試驗), 其線性可加效應(yīng)模型如下:

式中,Yabij為某小區(qū)的觀察值, 如全N淋失量;μ為總體均值;αa為因素 A第a水平的效應(yīng);βb為因素B第b水平的效應(yīng); (αβ)ab為因素A與因素B在ab水平上的交互作用效應(yīng);Bk為第k個區(qū)組效應(yīng);δi(ab)為第i個對象(小區(qū))在ab水平上的效應(yīng), 即對象間誤差, 滿足i.i.d. N(0,σ2主間);rj為重復(fù)測量因素在時間點j的效應(yīng); (αr)aj、(βr)bj分別為因素 A、B 與時間點的交互作用; (αβr)abj為三因素交互作用;eabij為誤差項, 滿足i.i.d. N(0,σe2)。Bk、δi(ab)、eabij為隨機效應(yīng),其余各項為固定效應(yīng)。δi(ab)、eabij相互獨立。

對于一因素三水平隨機區(qū)組設(shè)計(田間試驗),其線性可加效應(yīng)模型如下:

式中,αa為處理因素A第a水平的效應(yīng);Bk為第k個區(qū)組效應(yīng);δia為對象間誤差, 滿足i.i.d. N(0,σ2主間);rj為重復(fù)測量因素在時間點j的效應(yīng); (αr)aj為因素A與時間點的交互作用;eaij為誤差項, 滿足i.i.d. N(0,σe2)。

1.4 協(xié)方差結(jié)構(gòu)的篩選

在重復(fù)測量試驗中, 某一測量時間點上測定值變異的大小構(gòu)成方差, 2個不同時間點上測定值相互變異的大小用協(xié)方差來度量。如果時間點I上的取值不影響時間點II上的取值, 則協(xié)方差為0, 反之則大于0。運用GLIMMIX程序分析重復(fù)測量的數(shù)據(jù)主要分兩步, 首先是挑選協(xié)方差結(jié)構(gòu), 其次是進行方差分析和均值比較。對于土柱試驗, 土柱(core)為受試對象(subject)。運用GLIMMIX挑選最恰當?shù)膮f(xié)方差結(jié)構(gòu)的參考SAS程序如下。

程序說明: 輸入代碼時注意分號為英文字符。數(shù)據(jù)存放在E:數(shù)據(jù)N.xls sheet1里, 格式見表1。程序(1)讀取 Excel文件的數(shù)據(jù)。程序(2)～(8)內(nèi)容基本相同, 僅在“type =”后的協(xié)方差結(jié)構(gòu)選項不同。對某些協(xié)方差結(jié)構(gòu), 包括一階自回歸[autoregressive(1),AR(1)]、循環(huán)相關(guān)(toeplitz, TOEP)、一階前依賴[antedependence, ANTE(1)]等協(xié)方差結(jié)構(gòu), 需要考慮對象間誤差, 因此需要增加random block core(Rain N)語句[15]。程序(2)～(8)調(diào)用 GLIMMIX 過程, 分別采用方差分量結(jié)構(gòu)(variance components, VC)、復(fù)合對稱結(jié)構(gòu)(compound symmetry, CS)、不規(guī)則結(jié)構(gòu)(unstructured, UN)、空間冪相關(guān)結(jié)構(gòu)[space power,SP(POW)]、一階自回歸[AR(1)]、循環(huán)相關(guān)結(jié)構(gòu)(TOEP)、一階前依賴結(jié)構(gòu)[ANTE(1)]模型進行方差分析。Class語句列出所有的分類變量。Model語句根據(jù)上文(1)式編寫, Model語句后僅需要列出固定效應(yīng)[4], 注意此與 mixed的不同。Rain|N|times表示 3個因素的主效應(yīng)、2級和 3級交互作用效應(yīng)的線性組合。采用 KR法對標準誤和自由度進行修正[16]。random_residual_語句相當于 mixed模型中的repeated語句[5]。選項type指定協(xié)方差結(jié)構(gòu)類型。選項sub指定數(shù)據(jù)集中的受試對象(subject), 本例中為土柱(core), 即小區(qū), 若為單因素試驗, 直接指定對象名稱; 若為多因素試驗, 則在對象名稱后加因素名稱, 并加“()”。ods output語句輸出模型擬合信息fitstatistics。程序(9)建立宏rename, 對數(shù)據(jù)集中已有變量名value更名, 否則原有數(shù)據(jù)列value的數(shù)據(jù)會被覆蓋[17]。程序(10)對不同協(xié)方差結(jié)構(gòu)擬合參數(shù)進行合并, 以便程序(11)輸出擬合結(jié)果, 為挑選最佳協(xié)方差結(jié)構(gòu)提供依據(jù)。

對于田間試驗, SAS代碼與上述類似, 現(xiàn)就VC結(jié)構(gòu)的代碼列出如右上:

這里, rate為生物炭用量, block為區(qū)組, plot為小區(qū), 即試驗對象, time為取樣時間。由于本田間試驗取樣時間間隔不相等, 對于 SP(POW)協(xié)方差結(jié)構(gòu),需要計算時間間隔變量, 本列以天數(shù) days表示, 并增加語句 “random _residual_/type=SP(POW)(days)sub=plot;”。對于需要增加 random 語句的協(xié)方差結(jié)構(gòu), 包括 SP(POW)、ANTE(1)、AR(1)+RE、TOEP,random 語句寫法為“random plot block block*time;”。

1.5 方差分析和均值比較

挑選出最恰當?shù)膮f(xié)方差結(jié)構(gòu)后, 就可以運用GLIMMIX進行方差分析和均值比較了。土柱試驗以一階前依賴結(jié)構(gòu) ANTE(1)模型為例進行說明, 參考SAS程序如下。

程序說明: 程序的開始部分讀入數(shù)據(jù), 代碼同上文的程序I (1), 讀者可拷貝粘貼。程序(2)以一階前依賴協(xié)方差結(jié)構(gòu)ANTE(1)模型進行F檢驗, 代碼基本上同程序I (8)。輸出結(jié)果主要包括擬合情況和F檢驗結(jié)果(表4)。語句(3)輸出所有因子主效應(yīng)及交互作用最小二乘均值及其差值的t檢驗結(jié)果, 由于輸出內(nèi)容龐大, 本例為83頁文檔, 許多內(nèi)容我們不用關(guān)注。為了減少篇幅, 突出重點, 建議將該語句注釋掉。利用特定的 lsmestimate語句進行均值比較。用法可參閱文獻[5]、文獻[18]和查閱SAS幫助文檔和文獻[4]。

語句(4)檢驗第6次淋洗N1與N2處理全N淋失量差異顯著性, 即進行特定時間兩處理水平間的比較。語句(4a)和(4b)作用完全相同, 前者采用傳統(tǒng)的定位式(positional)寫法, 較為復(fù)雜; 后者采用新式的非定位式(non-positional)寫法[18], 更為簡單, 讀者可以任選其一。語句(5)檢驗N2處理第6次淋洗與第8次淋洗全N淋失量差異顯著性, 即特定處理隨時間的變化差異; 語句(6)檢驗N1處理第7次淋洗與N2處理第9次淋洗全氮淋失量差異顯著性, 即兩處理水平在各特定時間點的差異, 由此可以找出最優(yōu)處理組合; 語句(7)檢驗 3種氮水平處理平均全氮淋失量, 第 8次淋洗與第 10次淋洗的差異顯著性,即所有處理平均隨時間的變化差異是否顯著; 語句(8)檢驗N1與N2在所有時間點平均, 全氮淋失量的差異顯著性, 即兩處理水平在所有時間點平均的差異, 輸出結(jié)果見表5。讀者可以根據(jù)需要, 進行其他處理組合間的差異顯著性檢驗。

語句(9)以淋洗次數(shù) times為橫坐標, 總氮淋失量為縱坐標, 輸出類似圖 1的折線圖。“sliceby=N”表示按N水平進行分類。Join將各點連成線。cl表示輸出95%的置信區(qū)間(confident interval)。圖1中誤差線為標準誤, 該語句輸出的誤差線為 95%置信區(qū)間。NLOPTIONS語句是為了將 PROC GLIMMIX的優(yōu)化技術(shù)設(shè)置為與 mixed 程序(Newton-Raphson算法)相同的技術(shù)[5]。

田間試驗由于處理組合較少, 我們可以所有處理組合均值的差異輸出, 不需要lsmestimate, 以VC結(jié)構(gòu)為例, 僅需要在代碼(12)的run之前增加如下語句:由于采用極大似然分析, 缺區(qū)時自動按不等重復(fù)進行, 如本例中的缺區(qū)處理按 2個重復(fù)進行計算,代碼不需要改動。

2 結(jié)果與分析

2.1 不同協(xié)方差結(jié)構(gòu)模型引起的F檢驗結(jié)果差異

SAS程序I輸出了F檢驗結(jié)果, 表 2是對土柱試驗數(shù)據(jù)7種協(xié)方差結(jié)構(gòu)模型F檢驗結(jié)果的匯總?？梢钥闯? 協(xié)方差結(jié)構(gòu)模型對F檢驗的結(jié)果影響較大。降雨強度與施肥方式之間的交互作用(Rain*N)雖然P值都大于0.05, 但差異較大, 最低的為0.2375,最高的為0.3161, 相差33.09%。降雨強度與淋洗次數(shù)之間的交互作用效應(yīng)(Rain*Times) UN和ANTE(1)不顯著, 而其余4種模型顯示極顯著。因此, 選用恰當?shù)膮f(xié)方差結(jié)構(gòu)模型進行F檢驗很關(guān)鍵, 否則有可能得出錯誤的結(jié)論。

2.2 協(xié)方差結(jié)構(gòu)模型的優(yōu)選

根據(jù)統(tǒng)計學理論, 選用協(xié)方差結(jié)構(gòu)模型時, 可參考赤池信息準則(akaike information criterion,AIC)、為小樣本修正的赤池信息準則(akaike information corrected criterion, AICC)、修正的赤池信息準則(corrected akaike information criterion, CAIC)、貝葉斯信息準則(bayesian information criterion, BIC)、漢南-奎因信息準則(hannan–quinn information criterion, HQIC)、-2殘差對數(shù)似然值準則(–2 res log likelihood, -2logL)[1,4-5,15,17]等, 值越小表示擬合性越好, 如果相近, 可通過 χ2檢驗[17]并結(jié)合試驗本身的特點進行判斷。另外, 協(xié)方差結(jié)構(gòu)模型需要估算的參數(shù)個數(shù)越少越好。從表3可以看出, 土柱試驗UN和ANTE(1)模型各準則值明顯低于其余5種協(xié)方差模型的值, 應(yīng)優(yōu)先考慮。UN模型協(xié)方差矩陣中需要估算的參數(shù)個數(shù)在所有模型中最多, 為n(n–1)/2,n為重復(fù)測量次數(shù), 計算時可能無法收斂, 估算無法完成[5]。本例中n=7, 參數(shù)個數(shù)為21個; ANTE(1)模型需要估算的參數(shù)為2n-1, 本例中為15。綜合考慮,選用ANTE(1)模型進行F檢驗和均值比較。

田間試驗結(jié)果顯示, UN結(jié)構(gòu)無法收斂, 無法計算擬合參數(shù), 其余 6種協(xié)方差結(jié)構(gòu)模型均收斂。綜合比較, VC結(jié)構(gòu)模型最合適。

表2 采用不同協(xié)方差結(jié)構(gòu)模型時F檢驗P值(III型檢驗)Table 2 P-value for F test with various covariance structures (III)

表3 不同協(xié)方差結(jié)構(gòu)模型擬合性(土柱試驗)Table 3 Model fit statistics with various covariance structures (soil column experiment)

2.3 F檢驗結(jié)果

選用ANTE(1)模型并運用GLIMMIX對土柱試驗數(shù)據(jù)進行F檢驗(SAS程序II), 結(jié)果見表4。可以看出, 區(qū)組間無顯著差異。降雨強度(Rain)、施肥方式(N)和淋洗次數(shù)(Times) 3級交互作用(Rain*N*Times)效果不顯著。二級交互作用效果降雨強度與施肥方式、降雨強度與淋洗次數(shù)不顯著。施肥方式與淋洗次數(shù)之間的交互作用效果極顯著(P<0.001), 主效應(yīng)降雨強度、施肥方式、淋洗次數(shù)不同水平間差異達極顯著水平(P<0.001), 因此有必須進一步分析。

表4 F檢驗結(jié)果(III型, ANTE1, 土柱試驗)Table 4 Output of F test (type III, ANTE1, soil column experiment)

田間試驗的F檢驗結(jié)果顯示, 生物炭用量、取樣時間以及他們的交互作用均達到極其顯著水準(P<0.01)。

2.4 處理組合均值比較結(jié)果

根據(jù)F檢驗結(jié)果, 土柱試驗數(shù)據(jù)僅僅需要對效應(yīng)顯著的均值進行顯著性檢驗, 但為了讓讀者全面掌握重復(fù)測量數(shù)據(jù)均值比較的方法, 下面針對“1.2分析比較”所提出的5種情況進行說明。表5是SAS程序II輸出的結(jié)果匯總?？梢钥闯? 5種比較都能進行。語句(4a)和(4b)所得結(jié)果完全一致, 說明采用傳統(tǒng)的定位句法和新式的非定位句法[18]能達到同樣的效果,但后者不需要用“0”占位, 更直觀簡潔, 書寫更方便,尤其是一些比較復(fù)雜的比較[18], 因此建議使用新式的非定位句法。

從表2可以看出, 第6次淋洗時, N2比N1的總氮淋失量低8.05 g (圖1), 接近5%差異顯著水準(P=0.078), 此為特定時間(第6次淋洗)兩處理水平(N1 vs.N2)間的比較; 語句(5)輸出的結(jié)果表明, N2處理總N淋失量第6次比第8次淋洗高10.94 g, 差異達1%極顯著水平(P= 0.0025), 我們還可以檢驗第6次與第8、9、10、11、12次淋洗的差異顯著性, 也可以對其他時間進行差異顯著性比較, 確定總N淋失量隨時間的變化規(guī)律。語句(6)輸出的結(jié)果表明, N1第7次淋洗比N2處理第9次淋洗全氮淋失量多6.38 g氮(表5和圖1), 差異極顯著(P= 0.0008)?？梢岳^續(xù)進行兩處理水平在各特定時間點的差異比較, 找出最優(yōu)處理組合。

語句(7)輸出結(jié)果表明, 3種氮水平平均全氮淋失量第8次與第10次淋洗相差0.94 g, 未達到5%差異顯著水平。由于N*Times交互作用效果顯著, 這種比較并不合適, 此處僅僅用來說明方法。語句(8)輸出的結(jié)果表明, 所有淋洗平均全氮淋失量N2比N1處理低1.78 g, 差異達5%顯著水平。同樣原因, 這種比較在此例中不合適, 僅用于方法說明。田間試驗均值間比較思路類似, 不再贅述。

表5 處理均值兩兩比較示例Table 5 Examples of mean comparisons between treatment means

3 討論

3.1 混合效應(yīng)模型

固定效應(yīng)是研究幾個樣本是否來自于同一個總體的推斷。隨機效應(yīng)是由2個或多個總體分別抽出的幾個樣本, 用這些樣本去研究總體是不是相同的推斷。如果一個模型中既包括固定效應(yīng), 又包括隨機效應(yīng), 為混合效應(yīng)模型[5]。本研究中, 3種降雨強度、3種施肥處理、不同測定時間以及它們的交互作用均為固定效應(yīng)。區(qū)組為隨機效應(yīng)。重復(fù)測量中的受試對象(小區(qū))是總體中的樣本, 擬通過比較這些樣本經(jīng)過不同處理后重復(fù)測量的效應(yīng), 以便推論到其所代表的總體中去, 因此重復(fù)測量中受試對象效應(yīng)為隨機效應(yīng)。

對象變異包括對象內(nèi)變異和對象間變異, 重復(fù)測量試驗不同對象間一般是相互獨立的, 但對同一對象不同時間點的測定幾乎總存在相關(guān)性, 間隔越近相關(guān)性越強, 因此, 重復(fù)測量不滿足一般線性模型方差分析對變量獨立性的要求?；旌夏Ｐ图瓤紤]了固定效應(yīng), 也考慮了隨機效應(yīng), 對于重復(fù)測量數(shù)據(jù), 通過預(yù)先指定協(xié)方差結(jié)構(gòu)模型, 能夠很好地進行重復(fù)測量數(shù)據(jù)的統(tǒng)計分析[3-5,15], 此在本例中也得到了驗證。

除了用 GLIMMIX進行重復(fù)測量的數(shù)據(jù)分析外, 也可以用 mixed程序進行。前者是 SAS公司后來開發(fā)的程序, 比mixed功能更強大, 例如, 能直接輸出類似圖1的圖; 進行均值比較時, 編寫代碼更簡單[5]。

3.2 用于重復(fù)測量方差分析的協(xié)方差結(jié)構(gòu)模型

選用適當?shù)膮f(xié)方差結(jié)構(gòu)模型對于重復(fù)測量數(shù)據(jù)的統(tǒng)計分析至關(guān)重要。忽視重復(fù)測量數(shù)據(jù)間的相關(guān)性而采用過于簡單的協(xié)方差模型, 會導致標準誤偏低, 處理效應(yīng)會犯I型錯誤; 而模型過于復(fù)雜, 檢驗效能會受到嚴重影響[1,3,19]。早已有研究表明, 重復(fù)測量分析結(jié)果會因協(xié)方差模型選擇不當而嚴重受損[19]。盡管特定數(shù)據(jù)集的真實協(xié)方差結(jié)構(gòu)鮮為人知,但必須指定接近的協(xié)方差模型才能獲得可靠的結(jié)果。Guerin等[20]的研究表明, 只要選擇合適的協(xié)方差結(jié)構(gòu), 采用混合模型進行重復(fù)測量的統(tǒng)計分析總能得到正確的結(jié)果。

SAS軟件提供了30多種協(xié)方差結(jié)構(gòu)模型[4], 讀者可以根據(jù)需要選用。本例通過方差分量結(jié)構(gòu)(VC)等7種結(jié)構(gòu)進行模擬。VC又稱簡單方差結(jié)構(gòu), 他假定同一對象各時間點重復(fù)測量值相互獨立, 在協(xié)方差矩陣中主對角線元素均為σ2, 其他為0 (圖2), 該結(jié)構(gòu)僅有一個參數(shù)σ需要估算[5]。然而, 這種情況在重復(fù)測量試驗中極少出現(xiàn)。

復(fù)合對稱結(jié)構(gòu)(compound symmetry, CS)主對角線元素為其他為σd, 這里σd為對象間方差,σe為對象內(nèi)方差, 需要對這2個參數(shù)進行估計。對于重復(fù)測量數(shù)據(jù), 這是最簡單的方差結(jié)構(gòu)。它假定不同重復(fù)測量點數(shù)據(jù)之間具有恒定的相關(guān)性, 與重復(fù)測量點之間的距離無關(guān)。過去, 我們對重復(fù)測量試驗數(shù)據(jù)進行球型檢驗(H-F檢驗), 如果滿足就可以把重復(fù)測量因素當成裂區(qū)設(shè)計來進行方差分析[9,21-22],實際上就是檢驗協(xié)方差是否與測量距離有關(guān)。如果檢驗未通過, 為了減少犯I類錯誤的幾率, 常常采用校正的方法, 但這些方法的可靠性仍然不高, 更可靠的辦法是采用恰當?shù)膮f(xié)方差結(jié)構(gòu), 按照一般線性混合模型的方法進行分析。在重復(fù)測量農(nóng)業(yè)試驗中,能滿足復(fù)合對稱結(jié)構(gòu)條件的情況不多, 除非重復(fù)測量之間的時間距離特別長, 影響可以忽略不計。

一階前依賴結(jié)構(gòu)[ANTE(1)]允許不同重復(fù)測量點間的方差不同, 以及不同測量點對之間相關(guān)性和協(xié)方差不等。此很符合本研究的實際情況, 在所有協(xié)方差結(jié)構(gòu)模型中, 該模型最優(yōu)(表 3), 進一步證實了該模型的實用性。使用 ANTE(1)時, 測量點必需要按先后順序正確排列, 但各測量點之間時間間隔不必相等。這也很符合一般作物類試驗的實際情況。本研究中對不規(guī)則(UN)、空間冪相關(guān)[SP(POW)]、一階自回歸[AR(1)]、循環(huán)相關(guān)(TOEP)等協(xié)方差結(jié)構(gòu)模型進行了嘗試, 限于篇幅, 此處不進行深入討論。有興趣的讀者請參閱文獻[1]、[5]、[15]和[19]。

選擇協(xié)方差結(jié)構(gòu)模型時, 首先應(yīng)排除明顯無意義的協(xié)方差結(jié)構(gòu)。如田間試驗中, 同一小區(qū)不同時間重復(fù)測量數(shù)據(jù)一般具有相關(guān)性, VC模型不合適,應(yīng)排除。對于時間間隔不等的情況, AR(1)結(jié)構(gòu)肯定不合適[1]。下一步, 可以以時間間隔為橫坐標, 以對象內(nèi)協(xié)方差為縱坐標作圖, 考察協(xié)方差的變化規(guī)律,初步判定合適的協(xié)方差結(jié)構(gòu)。最后, 查看如表 3所示的擬合性信息, 必要時進行卡方檢驗, 確定最合適的協(xié)方差結(jié)構(gòu)。

3.3 與裂區(qū)設(shè)計方差分析或多變量方差分析的比較

在進行重復(fù)測量數(shù)據(jù)的統(tǒng)計分析時, 許多教材和學術(shù)論文中主要采用以下4種方法。(1)在各時間點分別進行F檢驗和均值比較, 該方法孤立地看待各時點的數(shù)據(jù), 完全忽視了時間因素, 沒有充分利用觀察對象在不同觀察時點間的內(nèi)在聯(lián)系, 無法發(fā)現(xiàn)隨時間變化的趨勢。在農(nóng)業(yè)類等已發(fā)表的論文中非常普遍[11,23-25], 是對信息資源的極大浪費[26]。(2)裂區(qū)設(shè)計分析法, 把重復(fù)測量因子作為副區(qū)因子,把不同時間點視為副處理水平進行分析[3,9,12], 該方法需要滿足球?qū)ΨQ條件[21-22], 即滿足復(fù)合對稱協(xié)方差結(jié)構(gòu)的條件。在農(nóng)業(yè)試驗中, 協(xié)方差會隨時間間隔的大小發(fā)生變化, 條件很難滿足, 采用此方法會增大犯 I 類錯誤的概率[1,4-5,12,27]。雖然可以進行校正, 但結(jié)果往往不理想, 甚至會得出錯誤的結(jié)論[13],且該方法沒有觸及協(xié)方差不等的實質(zhì), 而是對問題進行了回避, 沒有從根本上解決問題[3,15]。(3)多變量統(tǒng)計分析法, 該方法對重復(fù)測定各時間點的相關(guān)性要求不高, 僅要求每對時點間的相關(guān)性唯一, 但采用該方法分析重復(fù)測量數(shù)據(jù)會浪費大量信息, 損失統(tǒng)計功效。另外, 如果受試對象即使只有一個時間點的數(shù)據(jù)缺失(此在農(nóng)學類試驗中常常出現(xiàn)), 也必須刪除該對象的全部數(shù)據(jù), 造成信息的不完整[8,15],該方法也是回避了不同時間點數(shù)據(jù)有相關(guān)性的問題,不是一個理想的方法[13]。(4)混合模型(mixed)分析法,該方法直接處理重復(fù)測量中數(shù)據(jù)的相關(guān)性問題[5],允許不同時間點數(shù)據(jù)存在相關(guān)性, 不必對自由度進行校正, 可以處理缺值問題, 允許時間點不等距[4,11,28],完全符合重復(fù)測量試驗的特點, 是一種理想的方法。通過事先指定一個恰當?shù)膮f(xié)方差結(jié)構(gòu), 根據(jù)此結(jié)構(gòu)計算各處理水平的最小二乘均值及其差值。協(xié)方差結(jié)構(gòu)不同, 均值和標準誤計算結(jié)果也不同, 尤其是對非平衡設(shè)計[29]。

廣義線性混合模型(GLIMMIX)是近年來SAS在mixed模型的基礎(chǔ)上開發(fā)的程序模塊, 是對mixed的擴展和改善, 均值間比較的句法更簡單, 能很方便的輸出圖形等[5], 正如本例所示, 是目前進行重復(fù)測量數(shù)據(jù)統(tǒng)計分析的有力工具, 為重復(fù)測量數(shù)據(jù)分析提供了極大方便。

4 結(jié)論

重復(fù)測量試驗對同一對象進行多次觀測, 各數(shù)據(jù)點之間存在自相關(guān)性。用傳統(tǒng)的裂區(qū)設(shè)計、多變量方差分析會造成數(shù)據(jù)的信息浪費、統(tǒng)計功效降低、缺區(qū)無法處理等問題, 甚至會導致錯誤的結(jié)論。廣義線性混合模型 GLIMMIX能很好地處理相關(guān)性問題, 功能強大, 結(jié)構(gòu)可靠性高, 使用簡單, 允許缺區(qū),是進行重復(fù)測量試驗數(shù)據(jù)方差分析和均值比較理想的方法, 值得推廣運用。