用SAS的mixed過程擬合林分的線性差分生長模型
2014-01-13 07:10:16王慶豐倪成才駱國民
王慶豐 倪成才 駱國民
【摘 要】本研究的目的在于研究如何用SAS的proc mixed過程擬合線性代數(shù)差分模型。所用數(shù)據(jù)來源于148個(gè)集約經(jīng)營火炬松人工林。直接擬合了一個(gè)胸高斷面積的收獲模型，而非代數(shù)差分生長模型。模型擬合過程如下：i）.同時(shí)確定隨林分變化的參數(shù)和最優(yōu)擬合的方差結(jié)構(gòu)模型；ii）.依據(jù)AIC、BIC和極大似然比檢驗(yàn)化簡期望模型；iii）.用代數(shù)差分法將擬合的收獲模型轉(zhuǎn)化為代數(shù)差分生長模型。
【關(guān)鍵詞】線性代數(shù)差分模型；mixed過程；模型篩選；林分生長與收獲預(yù)估
0.前言
在林分生長與收獲預(yù)估的模型中，差分生長模型得到了廣泛的應(yīng)用。線性差分模型基本上為Schumacher模型的變型，廣泛應(yīng)用于林分蓄積、胸高斷面積的建模，以及單位面積株數(shù)和優(yōu)勢木樹高生長模型。
差分生長模型擬合方法有“直接最小二乘估計(jì)法”和“分類變量回歸法”[1]。一般認(rèn)為后者可以獲得近似無偏的估計(jì)，而前者則導(dǎo)致檢驗(yàn)統(tǒng)計(jì)量如RMSE的失真[2]。傳統(tǒng)上差分生長模型的擬合主要是直接擬合差分生長模型，然后根據(jù)擬合統(tǒng)計(jì)量如RMSE、R2等確定最優(yōu)擬合模型。與傳統(tǒng)方法不同，本文直接擬合生長模型，在獲得參數(shù)估計(jì)值后，再用代數(shù)差分法導(dǎo)出相應(yīng)的代數(shù)差分生長模型。這樣做的優(yōu)越之處在于非常便于對期望模型和方差結(jié)構(gòu)模型進(jìn)行篩選。更為重要的是，可以通過模型擬合識別最適合的隨林分變化參數(shù)。
本文詳細(xì)討論了如何用“分類變量回歸法”和SAS的mixed過程擬合代數(shù)差分生長模型，可簡述如下：i）.直接以生長收獲模型為對象，同時(shí)確定一個(gè)隨林分變化的參數(shù)和最優(yōu)擬合方差結(jié)構(gòu)模型；ii）. 保持方差結(jié)構(gòu)模型不變，根據(jù)擬合統(tǒng)計(jì)量逐步化簡期望模型；iii）.在確定最優(yōu)擬合的期望模型后，運(yùn)用代數(shù)差分法導(dǎo)出相對應(yīng)的代數(shù)差分生長模型。所有擬合與篩選均用SAS的mixed過程完成，并給出了詳細(xì)的SAS代碼和代碼解釋。
1.方法與材料
1.1數(shù)據(jù)
數(shù)據(jù)來源于148個(gè)集約經(jīng)營的火炬松實(shí)驗(yàn)人工林逐年觀測的固定樣地?cái)?shù)據(jù)（樣地約0.152公頃）。SAS的數(shù)據(jù)集basal的內(nèi)容如表（1）。
表1 模型擬合的基本數(shù)據(jù)結(jié)構(gòu)
Table 1 data structure for mode fitting
age=林分年齡；fert=經(jīng)營措施，分別取值為H=施加除草劑以控制競爭植物、F=施肥以增加土壤肥力、HF=除草劑和施肥并用、C=對照；code=樣地代碼，每個(gè)樣地有一個(gè)唯一代碼；logba=樣地胸高斷面積的自然對數(shù)值；logdh=樣地優(yōu)勢木樹高的自然對數(shù)值；iage=林分年齡的倒數(shù)；logtpa=樣地株數(shù)的自然對數(shù)值。共有148個(gè)樣地?cái)?shù)據(jù)，1491個(gè)記錄。
1.2數(shù)學(xué)模型
所考慮的胸高斷面積收獲的數(shù)學(xué)模型為
E（lny）=α+αlnN+αlnH+α+α （1）
這里lny、lnH和lnN分別對應(yīng)logba、logdh和logtpa；1/t對應(yīng)iage，其余為模型參數(shù)。通過不同假設(shè)，以上模型可以導(dǎo)出很多廣泛應(yīng)用的差分生長模型。例如假設(shè)α0為隨林分變化的參數(shù)，運(yùn)用代數(shù)差分法可以導(dǎo)出Pienaar和Shiver （1986）的胸高斷面積模型[3]；假設(shè)α4=0且α0為隨林分變化的參數(shù)則可導(dǎo)出Forss等人（1996）提出的差分生長模型[4]；假設(shè)α4=0和α2=0且α3為隨林分變化的參數(shù)，則可以導(dǎo)出Souter （1986）的模型[5]。
1.3參數(shù)估計(jì)方法
從以上代數(shù)差分法的討論中，總是假設(shè)模型（1）中有一個(gè)參數(shù)隨林分變化而變化，因而使用“分類變量回歸法”以便考慮這一特點(diǎn)。如果用SAS的reg過程，必須構(gòu)造特殊的數(shù)據(jù)結(jié)構(gòu)。用SAS的reg過程擬合線性差分生長模型有兩點(diǎn)不足之處：1）.需要構(gòu)造復(fù)雜的數(shù)據(jù)文件。以本研究的為例，數(shù)據(jù)文件中需要增加額外的148個(gè)變量；2）.reg過程無法考慮重復(fù)觀測數(shù)據(jù)的自相關(guān)性和異質(zhì)方差結(jié)構(gòu)。由于所考慮的胸高斷面積模型均為線性模型，而且數(shù)據(jù)為典型的重復(fù)觀測數(shù)據(jù)，因而模型的擬合與篩選均用SAS的proc mixed過程。
1.4模型篩選和篩選指標(biāo)選用
模型篩選包括期望模型和最優(yōu)擬合的自相關(guān)與異質(zhì)方差結(jié)構(gòu)模型的篩選。Ngo和Brand（1997）推薦了兩種模型選擇方法[6]。一種方法就是首先列出所有可用的數(shù)學(xué)期望模型和所有可能的方差結(jié)構(gòu)模型，擬合二者間的所有組合，然后根據(jù)模型擬合指標(biāo)選擇最優(yōu)擬合模型。另一種方法為Wolfinger和Diggle提出的方法，即首先考慮最復(fù)雜的期望模型，并保持期望模型不變，然后選擇最優(yōu)擬合的方差結(jié)構(gòu)模型。一旦選出最優(yōu)擬合的方差模型，保持方差模型不變，再逐步簡化期望模型。模型篩選指標(biāo)較多，但常用的指標(biāo)為極大似然比檢驗(yàn)（LRT）、AIC、BIC（或稱為SBC）。LRT用于嵌套的模型，而AIC和BIC用于非嵌套模型。
2.結(jié)果與分析
2.1殘差的方差結(jié)構(gòu)模型
為考慮經(jīng)營措施對胸高斷面積生長過程的影響，將模型（1）改寫如下：
這里參數(shù)?為不同經(jīng)營措施對截距和各項(xiàng)回歸系數(shù)的影響，下標(biāo)k分別取值為C、F、H和HF，對應(yīng)四種經(jīng)營措施；由于差分生長模型需要確定一個(gè)隨林分變化而變化的參數(shù)，因而模型擬合的首要任務(wù)就是識別該參數(shù)。假設(shè)“1/t”的回歸系數(shù)隨林分變化而變化，那么模型（2）則可以表示如下：
這里βi為隨林分變化的參數(shù)。
根據(jù)以往的研究，所考慮的方差模型為：CS模型（同一林分不同年齡的觀測值方差相同，而且相關(guān)系數(shù)為常數(shù)）、AR（1）模型（不同年齡觀測值的方差相同，但相關(guān)系數(shù)為一階自相關(guān)模型）、ARH（1）模型（相關(guān)系數(shù)為一階自相關(guān)模型， 但不同年齡觀測值的方差相異）和ARMA（1，1） 模型（不同年齡觀測值的方差相同，但相關(guān)系數(shù)為一階自相關(guān)移動(dòng)平均模型）。
由于模型（2）中共有5個(gè)參數(shù)可以考慮為隨林分而變化的參數(shù)（b0～b4），因而相對應(yīng)地共有5個(gè)期望模型，加上所考慮的4個(gè)方差結(jié)構(gòu)模型，所有組合共計(jì)為20個(gè)模型。首先擬合這20個(gè)模型，根據(jù)擬合統(tǒng)計(jì)量AIC、BIC確定隨林分變化的參數(shù)及最優(yōu)擬合的方差結(jié)構(gòu)模型，這部分與Ngo和Brand的方法一致；然后保持方差結(jié)構(gòu)模型不變，根據(jù)LRT、AIC、BIC剔除回歸效果不顯著的因子，包括經(jīng)營措施效果等，即化簡期望模型。而這一部分則與Wolfinger和Diggle的方法一致。方差結(jié)構(gòu)篩選的基本SAS代碼如下：
以上代碼及以后代碼中的粗體字為SAS系統(tǒng)關(guān)鍵詞。class語句聲明了3個(gè)分類變量，分別為fert、code及age。model語句中既包含離散型的分類變量，又包含連續(xù)型變量（如logba、logtp等，即未在class語句中聲明的變量）。以上代碼的model語句意義如下：
i）.代碼中model語句中的分類變量表示該變量對回歸模型截距的影響。例如model語句中fert表示不同經(jīng)營措施下的logba的回歸模型的截距不同，隨fert變量的取值變化而變化。
ii）.如果一個(gè)連續(xù)型變量與一個(gè)分類變量的乘積出現(xiàn)model語句中，且僅以該形式出現(xiàn)，如code*iage，則表示該連續(xù)型變量的回歸系數(shù)隨分類變量的取值變化而變化。code是一個(gè)分類變量，每個(gè)取值對應(yīng)一個(gè)樣地，而iage為一個(gè)連續(xù)型變量，那么code*iage則表示每個(gè)樣地的iage項(xiàng)的回歸系數(shù)均不相同，各有其取值。
iii）.如果一個(gè)連續(xù)型變量與一個(gè)分類變量的乘積出現(xiàn)，而且該連續(xù)型變量同時(shí)獨(dú)立出現(xiàn)在模型中，同樣表示該連續(xù)型變量的回歸系數(shù)隨分類變量的取值變化而變化，但與ii）的意義有所不同。例如logdh和fert*logdh同時(shí)出現(xiàn)在model語句中，表示logdh的回歸系數(shù)隨fert取值的變化而變化，而且每個(gè)fert取值的回歸系數(shù)均表示為兩個(gè)分量之和。一個(gè)分量為參照回歸系數(shù)，另一個(gè)分量為與參照回歸系數(shù)相比的增量（可負(fù)可正）。mixed過程按照英文字母的排列順序，以fert的最后一個(gè)取值（即HF）為參照，其回歸系數(shù)即為參照回歸系數(shù)。而fert其它取值的回歸系數(shù)則表示為參照回歸系數(shù)（即HF）和與之相比增量之和。
以上model語句的對應(yīng)的數(shù)學(xué)模型如下：
以上四個(gè)模型分別對應(yīng)經(jīng)營措施C、F、H和HF的回歸模型。參數(shù)α0、α1、α2、α4為在HF經(jīng)營措施下，林分胸高斷面積的回歸系數(shù)（參照回歸系數(shù)）；參數(shù)φi，k為k經(jīng)營措施的第i個(gè)回歸系數(shù)與HF相對應(yīng)的回歸系數(shù)相比之增量（k取值為C、H和F）。以對照C的lnH的系數(shù)為例，其回歸系數(shù)為α2+φ2，C。這里φ2，C就是C的lnH的回歸系數(shù)與HF的lnH的回歸系數(shù)相比的增量。以上代碼中，fert*logtpa、fert*logdh及fert*logdh*iage分別對應(yīng)各項(xiàng)回歸系數(shù)與HF的回歸系數(shù)相比的增量，即φ；logtpa、logdh及l(fā)ogdh*iage則對應(yīng)參照回歸系數(shù)，即HF的回歸系數(shù)α0、α1、α2和α4等。
依據(jù)模型擬合統(tǒng)計(jì)量AIC和BIC值（其值越小，擬合效果越好）選擇模型，指定1/t的回歸系數(shù)為隨林分變化的參數(shù)和arh（1）為方差結(jié)構(gòu)模型時(shí)，AIC和BIC均取得了最小值，分別為-3629.9和-3585.0。因而arh（1）為最優(yōu)擬合的方差結(jié)構(gòu)模型，并且可以確認(rèn)1/t的回歸系數(shù)為隨林分變化的參數(shù)。
2.2最優(yōu)擬合的期望模型
保持arh（1）結(jié)構(gòu)不變，根據(jù)AIC、BIC及極大似然比檢驗(yàn)的p值進(jìn)行期望模型的篩選。模型（3）的擬合結(jié)果表明，經(jīng)營措施H和HF各個(gè)參數(shù)差異并不顯著，φ0，H、φ1，H、φ2，H 和φ4，H是否為0的t檢驗(yàn)p值均大于0.05，說明經(jīng)營措施H的各項(xiàng)回歸系數(shù)與HF并無顯著差別。將二者的回歸系數(shù)合并為一，僅保留截距不同，可擬合如下模型：
這里bi，k=αi，k+φi，k（i=0，1，2，4；k=C，F(xiàn)，H和HF）；b1，p、b2，p、b4，p為經(jīng)營措施H和HF的回歸系數(shù)，b0，H和b0，HF則分別為二者的截距。模型（3）便于比較不同經(jīng)營措施的效果異同，而模型（4）則便于化簡期望模型，即判斷一個(gè)回歸因子對胸高斷面積的影響是否顯著。注意模型（3）和（4）的唯一區(qū)別僅在于將模型（3）中第3個(gè)和第4個(gè)模型合并為一個(gè)而已（僅截距不同）。為擬合以上模型，新引進(jìn)了一個(gè)分類變量“f”，其取值分別為C、F和H，其中經(jīng)營措施H和HF的樣地“f”變量值均為H。擬合模型（4）的SAS代碼如下：
模型（4）的AIC、BIC的值均小于模型（3）的AIC和BIC值，LRT檢驗(yàn)p值為0.1573，表明模型（4）與模型（3）擬合效果基本相當(dāng)，但模型（4）更為簡潔。由于參數(shù)的t檢驗(yàn)結(jié)果表明僅對照C的lnH/t參數(shù)與0相比差異顯著，因而剔除F、H和HF的lnH/t項(xiàng)（相應(yīng)的模型稱為模型5）。相對應(yīng)地，在SAS的數(shù)據(jù)文件中引入了一個(gè)新變量“check”，當(dāng)且僅當(dāng)林分的經(jīng)營措施為C取值為1，其它情形均取值為0。對應(yīng)的SAS代碼如下：
與模型（4）的SAS代碼相比，僅用check*logdh*iage 取代了f*logdh*iage。模型（5）的參數(shù)t檢驗(yàn)結(jié)果表明H和HF的回歸模型中l(wèi)nH的回歸系數(shù)b2，p與0比較差異不顯著（p值為0.1181），但剔除該項(xiàng)后的模型擬合的LRT檢驗(yàn)平值、AIC、AICC、BIC表現(xiàn)出一種相互矛盾的結(jié)果。AIC和AICC表明兩個(gè)模型擬合效果基本一致，而BIC和LRT檢驗(yàn)結(jié)果則傾向于選擇化簡后的模型。本研究最終以模型（5）為最優(yōu)擬合模型，部分SAS輸出的參數(shù)估計(jì)值見表（2）。
表2 模型（5）的參數(shù)估計(jì)值
Table 2 parameters estimates for model （5）（下轉(zhuǎn)第164頁）
（上接第109頁）*僅列出參數(shù)估計(jì)值部分，與code值相對應(yīng)148參數(shù)估計(jì)值未列于此表。
2.3代數(shù)差分生長模型的導(dǎo)出
根據(jù)表（2）的模型結(jié)構(gòu)，HF對應(yīng)的回歸模型為：
ln=+lnN+lnH+β （6）
假設(shè)t=t1時(shí)ln1已知，求解關(guān)于βi的表達(dá)式可有βi=t1（ln--lnN-lnH）；代入模型（6）中可得當(dāng)t=t2時(shí)的lny2的估計(jì)式，即：
ln=lny+++
其中l(wèi)n1用lny1代替。其它經(jīng)營措施的代數(shù)差分模型的導(dǎo)出基本相似，在此略過。
3.討論與結(jié)論
本研究的模型篩選方法與以往研究不同。以往對差分生長模型篩選方法為：擬合選取的差分生長模型，讓后根據(jù)擬合統(tǒng)計(jì)量如RMSE、R2選擇最優(yōu)擬合模型。本文從生長收獲模型入手，而不是從差分生長模型入手。首先確定最優(yōu)擬合的方差結(jié)構(gòu)模型，同時(shí)確定最適合作為隨林分變化的參數(shù)。然后保持方差結(jié)構(gòu)模型不變的前提下，通過AIC、BIC及極大似然比檢驗(yàn)的p值逐步化簡期望模型。所有模型擬合和篩選均用SAS的proc mixed過程。該過程不僅非常便于擬合有分類變量的回歸模型（例如經(jīng)營措施），而且可以非常方便地?cái)M合方差結(jié)構(gòu)模型。本文所討論的差分生長模型的擬合方法、篩選方法和擬合過程，均可以應(yīng)用于其它林分因子如蓄積、優(yōu)勢木樹高的線性差分生長模型。對于非線性模型，除了需要構(gòu)造特殊的數(shù)據(jù)文件外，其它內(nèi)容也基本相同。 [科]
【參考文獻(xiàn)】
[1]倪成才，劉春梅，丁俊峰，潘曉茹.差分生長模型預(yù)測誤差的分析[J].北京林業(yè)大學(xué)學(xué)報(bào)，2009，31（4）：1-6.
[2]NI C，ZHANG L.An Analysis and comparison of estimation methods for self-referencing equations[J].Canadian Journal of Forest Research，2007，37：1472-1484.
[3]PIENAAR L V， SHIVER B D.Basal area prediction and projection equations for pine plantations[J].Forest Science，1986，32（3）：626-633.
[4]FORSS，E.，GADOW，K.V.and Saborowski，J.1996.Growth models for unthinned Acacia mangium plantations in South Kalimantan，Indonesia[J]. Journal of Tropical Forest Science，8（4）：449-462.
[5]SOUTER，R.A.1986.Dynamic stand structure in thinned stands of naturally regenerated loblolly pine in the Georgia Piedmont.Ph.D.Thesis，University of Georgia，Athens，GA.
[6]Ngo L.，Brand R.，1997.Model Selection in Linear Mixed Effects Models Using SAS PROC MIXED.SAS Conference Proceedings：SAS Users Group International 22，March16-19，1997，San Diego，California.
