王慶豐　倪成才　駱國民

【摘 要】本研究的目的在于研究如何用SAS的proc mixed過程擬合線性代數(shù)差分模型。所用數(shù)據(jù)來源于148個集約經(jīng)營火炬松人工林。直接擬合了一個胸高斷面積的收獲模型，而非代數(shù)差分生長模型。模型擬合過程如下：i）.同時確定隨林分變化的參數(shù)和最優(yōu)擬合的方差結構模型；ii）.依據(jù)AIC、BIC和極大似然比檢驗化簡期望模型；iii）.用代數(shù)差分法將擬合的收獲模型轉(zhuǎn)化為代數(shù)差分生長模型。

【關鍵詞】線性代數(shù)差分模型；mixed過程；模型篩選；林分生長與收獲預估

0.前言

在林分生長與收獲預估的模型中，差分生長模型得到了廣泛的應用。線性差分模型基本上為Schumacher模型的變型，廣泛應用于林分蓄積、胸高斷面積的建模，以及單位面積株數(shù)和優(yōu)勢木樹高生長模型。

差分生長模型擬合方法有“直接最小二乘估計法”和“分類變量回歸法”[1]。一般認為后者可以獲得近似無偏的估計，而前者則導致檢驗統(tǒng)計量如RMSE的失真[2]。傳統(tǒng)上差分生長模型的擬合主要是直接擬合差分生長模型，然后根據(jù)擬合統(tǒng)計量如RMSE、R2等確定最優(yōu)擬合模型。與傳統(tǒng)方法不同，本文直接擬合生長模型，在獲得參數(shù)估計值后，再用代數(shù)差分法導出相應的代數(shù)差分生長模型。這樣做的優(yōu)越之處在于非常便于對期望模型和方差結構模型進行篩選。更為重要的是，可以通過模型擬合識別最適合的隨林分變化參數(shù)。

本文詳細討論了如何用“分類變量回歸法”和SAS的mixed過程擬合代數(shù)差分生長模型，可簡述如下：i）.直接以生長收獲模型為對象，同時確定一個隨林分變化的參數(shù)和最優(yōu)擬合方差結構模型；ii）. 保持方差結構模型不變，根據(jù)擬合統(tǒng)計量逐步化簡期望模型；iii）.在確定最優(yōu)擬合的期望模型后，運用代數(shù)差分法導出相對應的代數(shù)差分生長模型。所有擬合與篩選均用SAS的mixed過程完成，并給出了詳細的SAS代碼和代碼解釋。

1.方法與材料

1.1數(shù)據(jù)

數(shù)據(jù)來源于148個集約經(jīng)營的火炬松實驗人工林逐年觀測的固定樣地數(shù)據(jù)（樣地約0.152公頃）。SAS的數(shù)據(jù)集basal的內(nèi)容如表（1）。

表1 模型擬合的基本數(shù)據(jù)結構

Table 1 data structure for mode fitting

age=林分年齡；fert=經(jīng)營措施，分別取值為H=施加除草劑以控制競爭植物、F=施肥以增加土壤肥力、HF=除草劑和施肥并用、C=對照；code=樣地代碼，每個樣地有一個唯一代碼；logba=樣地胸高斷面積的自然對數(shù)值；logdh=樣地優(yōu)勢木樹高的自然對數(shù)值；iage=林分年齡的倒數(shù)；logtpa=樣地株數(shù)的自然對數(shù)值。共有148個樣地數(shù)據(jù)，1491個記錄。

1.2數(shù)學模型

所考慮的胸高斷面積收獲的數(shù)學模型為

E（lny）=α+αlnN+αlnH+α+α （1）

這里lny、lnH和lnN分別對應logba、logdh和logtpa；1/t對應iage，其余為模型參數(shù)。通過不同假設，以上模型可以導出很多廣泛應用的差分生長模型。例如假設α0為隨林分變化的參數(shù)，運用代數(shù)差分法可以導出Pienaar和Shiver （1986）的胸高斷面積模型[3]；假設α4=0且α0為隨林分變化的參數(shù)則可導出Forss等人（1996）提出的差分生長模型[4]；假設α4=0和α2=0且α3為隨林分變化的參數(shù)，則可以導出Souter （1986）的模型[5]。

1.3參數(shù)估計方法

從以上代數(shù)差分法的討論中，總是假設模型（1）中有一個參數(shù)隨林分變化而變化，因而使用“分類變量回歸法”以便考慮這一特點。如果用SAS的reg過程，必須構造特殊的數(shù)據(jù)結構。用SAS的reg過程擬合線性差分生長模型有兩點不足之處：1）.需要構造復雜的數(shù)據(jù)文件。以本研究的為例，數(shù)據(jù)文件中需要增加額外的148個變量；2）.reg過程無法考慮重復觀測數(shù)據(jù)的自相關性和異質(zhì)方差結構。由于所考慮的胸高斷面積模型均為線性模型，而且數(shù)據(jù)為典型的重復觀測數(shù)據(jù)，因而模型的擬合與篩選均用SAS的proc mixed過程。

1.4模型篩選和篩選指標選用

模型篩選包括期望模型和最優(yōu)擬合的自相關與異質(zhì)方差結構模型的篩選。Ngo和Brand（1997）推薦了兩種模型選擇方法[6]。一種方法就是首先列出所有可用的數(shù)學期望模型和所有可能的方差結構模型，擬合二者間的所有組合，然后根據(jù)模型擬合指標選擇最優(yōu)擬合模型。另一種方法為Wolfinger和Diggle提出的方法，即首先考慮最復雜的期望模型，并保持期望模型不變，然后選擇最優(yōu)擬合的方差結構模型。一旦選出最優(yōu)擬合的方差模型，保持方差模型不變，再逐步簡化期望模型。模型篩選指標較多，但常用的指標為極大似然比檢驗（LRT）、AIC、BIC（或稱為SBC）。LRT用于嵌套的模型，而AIC和BIC用于非嵌套模型。

2.結果與分析

2.1殘差的方差結構模型

為考慮經(jīng)營措施對胸高斷面積生長過程的影響，將模型（1）改寫如下：（2）

這里參數(shù)?為不同經(jīng)營措施對截距和各項回歸系數(shù)的影響，下標k分別取值為C、F、H和HF，對應四種經(jīng)營措施；由于差分生長模型需要確定一個隨林分變化而變化的參數(shù)，因而模型擬合的首要任務就是識別該參數(shù)。假設“1/t”的回歸系數(shù)隨林分變化而變化，那么模型（2）則可以表示如下：

這里βi為隨林分變化的參數(shù)。

根據(jù)以往的研究，所考慮的方差模型為：CS模型（同一林分不同年齡的觀測值方差相同，而且相關系數(shù)為常數(shù)）、AR（1）模型（不同年齡觀測值的方差相同，但相關系數(shù)為一階自相關模型）、ARH（1）模型（相關系數(shù)為一階自相關模型， 但不同年齡觀測值的方差相異）和ARMA（1，1） 模型（不同年齡觀測值的方差相同，但相關系數(shù)為一階自相關移動平均模型）。

由于模型（2）中共有5個參數(shù)可以考慮為隨林分而變化的參數(shù)（b0～b4），因而相對應地共有5個期望模型，加上所考慮的4個方差結構模型，所有組合共計為20個模型。首先擬合這20個模型，根據(jù)擬合統(tǒng)計量AIC、BIC確定隨林分變化的參數(shù)及最優(yōu)擬合的方差結構模型，這部分與Ngo和Brand的方法一致；然后保持方差結構模型不變，根據(jù)LRT、AIC、BIC剔除回歸效果不顯著的因子，包括經(jīng)營措施效果等，即化簡期望模型。而這一部分則與Wolfinger和Diggle的方法一致。方差結構篩選的基本SAS代碼如下：

以上代碼及以后代碼中的粗體字為SAS系統(tǒng)關鍵詞。class語句聲明了3個分類變量，分別為fert、code及age。model語句中既包含離散型的分類變量，又包含連續(xù)型變量（如logba、logtp等，即未在class語句中聲明的變量）。以上代碼的model語句意義如下：

i）.代碼中model語句中的分類變量表示該變量對回歸模型截距的影響。例如model語句中fert表示不同經(jīng)營措施下的logba的回歸模型的截距不同，隨fert變量的取值變化而變化。

ii）.如果一個連續(xù)型變量與一個分類變量的乘積出現(xiàn)model語句中，且僅以該形式出現(xiàn)，如code*iage，則表示該連續(xù)型變量的回歸系數(shù)隨分類變量的取值變化而變化。code是一個分類變量，每個取值對應一個樣地，而iage為一個連續(xù)型變量，那么code*iage則表示每個樣地的iage項的回歸系數(shù)均不相同，各有其取值。

iii）.如果一個連續(xù)型變量與一個分類變量的乘積出現(xiàn)，而且該連續(xù)型變量同時獨立出現(xiàn)在模型中，同樣表示該連續(xù)型變量的回歸系數(shù)隨分類變量的取值變化而變化，但與ii）的意義有所不同。例如logdh和fert*logdh同時出現(xiàn)在model語句中，表示logdh的回歸系數(shù)隨fert取值的變化而變化，而且每個fert取值的回歸系數(shù)均表示為兩個分量之和。一個分量為參照回歸系數(shù)，另一個分量為與參照回歸系數(shù)相比的增量（可負可正）。mixed過程按照英文字母的排列順序，以fert的最后一個取值（即HF）為參照，其回歸系數(shù)即為參照回歸系數(shù)。而fert其它取值的回歸系數(shù)則表示為參照回歸系數(shù)（即HF）和與之相比增量之和。

以上model語句的對應的數(shù)學模型如下：

以上四個模型分別對應經(jīng)營措施C、F、H和HF的回歸模型。參數(shù)α0、α1、α2、α4為在HF經(jīng)營措施下，林分胸高斷面積的回歸系數(shù)（參照回歸系數(shù)）；參數(shù)φi，k為k經(jīng)營措施的第i個回歸系數(shù)與HF相對應的回歸系數(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分別對應各項回歸系數(shù)與HF的回歸系數(shù)相比的增量，即φ；logtpa、logdh及l(fā)ogdh*iage則對應參照回歸系數(shù)，即HF的回歸系數(shù)α0、α1、α2和α4等。

依據(jù)模型擬合統(tǒng)計量AIC和BIC值（其值越小，擬合效果越好）選擇模型，指定1/t的回歸系數(shù)為隨林分變化的參數(shù)和arh（1）為方差結構模型時，AIC和BIC均取得了最小值，分別為-3629.9和-3585.0。因而arh（1）為最優(yōu)擬合的方差結構模型，并且可以確認1/t的回歸系數(shù)為隨林分變化的參數(shù)。

2.2最優(yōu)擬合的期望模型

保持arh（1）結構不變，根據(jù)AIC、BIC及極大似然比檢驗的p值進行期望模型的篩選。模型（3）的擬合結果表明，經(jīng)營措施H和HF各個參數(shù)差異并不顯著，φ0，H、φ1，H、φ2，H 和φ4，H是否為0的t檢驗p值均大于0.05，說明經(jīng)營措施H的各項回歸系數(shù)與HF并無顯著差別。將二者的回歸系數(shù)合并為一，僅保留截距不同，可擬合如下模型：

（4）

這里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）則便于化簡期望模型，即判斷一個回歸因子對胸高斷面積的影響是否顯著。注意模型（3）和（4）的唯一區(qū)別僅在于將模型（3）中第3個和第4個模型合并為一個而已（僅截距不同）。為擬合以上模型，新引進了一個分類變量“f”，其取值分別為C、F和H，其中經(jīng)營措施H和HF的樣地“f”變量值均為H。擬合模型（4）的SAS代碼如下：

模型（4）的AIC、BIC的值均小于模型（3）的AIC和BIC值，LRT檢驗p值為0.1573，表明模型（4）與模型（3）擬合效果基本相當，但模型（4）更為簡潔。由于參數(shù)的t檢驗結果表明僅對照C的lnH/t參數(shù)與0相比差異顯著，因而剔除F、H和HF的lnH/t項（相應的模型稱為模型5）。相對應地，在SAS的數(shù)據(jù)文件中引入了一個新變量“check”，當且僅當林分的經(jīng)營措施為C取值為1，其它情形均取值為0。對應的SAS代碼如下：

與模型（4）的SAS代碼相比，僅用check*logdh*iage 取代了f*logdh*iage。模型（5）的參數(shù)t檢驗結果表明H和HF的回歸模型中l(wèi)nH的回歸系數(shù)b2，p與0比較差異不顯著（p值為0.1181），但剔除該項后的模型擬合的LRT檢驗平值、AIC、AICC、BIC表現(xiàn)出一種相互矛盾的結果。AIC和AICC表明兩個模型擬合效果基本一致，而BIC和LRT檢驗結果則傾向于選擇化簡后的模型。本研究最終以模型（5）為最優(yōu)擬合模型，部分SAS輸出的參數(shù)估計值見表（2）。

表2 模型（5）的參數(shù)估計值

Table 2 parameters estimates for model （5）（下轉(zhuǎn)第164頁）

（上接第109頁）*僅列出參數(shù)估計值部分，與code值相對應148參數(shù)估計值未列于此表。

2.3代數(shù)差分生長模型的導出

根據(jù)表（2）的模型結構，HF對應的回歸模型為：

ln=+lnN+lnH+β （6）

假設t=t1時ln1已知，求解關于βi的表達式可有βi=t1（ln--lnN-lnH）；代入模型（6）中可得當t=t2時的lny2的估計式，即：

ln=lny+++

其中l(wèi)n1用lny1代替。其它經(jīng)營措施的代數(shù)差分模型的導出基本相似，在此略過。

3.討論與結論

本研究的模型篩選方法與以往研究不同。以往對差分生長模型篩選方法為：擬合選取的差分生長模型，讓后根據(jù)擬合統(tǒng)計量如RMSE、R2選擇最優(yōu)擬合模型。本文從生長收獲模型入手，而不是從差分生長模型入手。首先確定最優(yōu)擬合的方差結構模型，同時確定最適合作為隨林分變化的參數(shù)。然后保持方差結構模型不變的前提下，通過AIC、BIC及極大似然比檢驗的p值逐步化簡期望模型。所有模型擬合和篩選均用SAS的proc mixed過程。該過程不僅非常便于擬合有分類變量的回歸模型（例如經(jīng)營措施），而且可以非常方便地擬合方差結構模型。本文所討論的差分生長模型的擬合方法、篩選方法和擬合過程，均可以應用于其它林分因子如蓄積、優(yōu)勢木樹高的線性差分生長模型。對于非線性模型，除了需要構造特殊的數(shù)據(jù)文件外，其它內(nèi)容也基本相同。 [科]

【參考文獻】

[1]倪成才，劉春梅，丁俊峰，潘曉茹.差分生長模型預測誤差的分析[J].北京林業(yè)大學學報，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.