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

張久權 閆慧峰 褚繼登 李彩斌

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

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

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

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

1 材料與方法

1.1 試驗設計

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

表1 數據Microsoft Excel輸入表示例Table 1 Data entry example for Microsoft Excel

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

1.2 分析比較

圖1顯示了土柱試驗3種不同施肥方式下總N淋失量隨淋次數的變化, 淋洗次數為重復測量的時間因素。我們可能需要確定如下處理組合間的差異是否顯著: (1)特定時間兩處理水平間的差異,如第6次淋洗時, 全氮淋失量 N1與N2處理的差異; (2)特定處理隨時間的變化, 如 N2處理, 全氮淋失量第6次與第8次淋洗的差異; (3)兩處理水平在各特定時間點的差異, 以便確定最優處理組合,如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個時間因素)隨機區組設計(土柱試驗), 其線性可加效應模型如下:

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

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

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

1.4 協方差結構的篩選

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

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

對于田間試驗, SAS代碼與上述類似, 現就VC結構的代碼列出如右上:

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

1.5 方差分析和均值比較

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

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

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

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

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

2 結果與分析

2.1 不同協方差結構模型引起的F檢驗結果差異

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

2.2 協方差結構模型的優選

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

田間試驗結果顯示, UN結構無法收斂, 無法計算擬合參數, 其余 6種協方差結構模型均收斂。綜合比較, VC結構模型最合適。

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

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

2.3 F檢驗結果

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

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

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

2.4 處理組合均值比較結果

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

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

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

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

3 討論

3.1 混合效應模型

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

對象變異包括對象內變異和對象間變異, 重復測量試驗不同對象間一般是相互獨立的, 但對同一對象不同時間點的測定幾乎總存在相關性, 間隔越近相關性越強, 因此, 重復測量不滿足一般線性模型方差分析對變量獨立性的要求。混合模型既考慮了固定效應, 也考慮了隨機效應, 對于重復測量數據, 通過預先指定協方差結構模型, 能夠很好地進行重復測量數據的統計分析[3-5,15], 此在本例中也得到了驗證。

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

3.2 用于重復測量方差分析的協方差結構模型

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

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

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

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

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

3.3 與裂區設計方差分析或多變量方差分析的比較

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

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

4 結論

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