用SAS的mixed过程拟合林分的线性差分生长模型
时间：2020-12-06 12:05:44 来源：达达文档网 本文已影响 人 
王庆丰 倪成才 骆国民
【摘 要】本研究的目的在于研究如何用SAS的proc mixed过程拟合线性代数差分模型。所用数据来源于148个集约经营火炬松人工林。直接拟合了一个胸高断面积的收获模型，而非代数差分生长模型。模型拟合过程如下：i）.同时确定随林分变化的参数和最优拟合的方差结构模型；
ii）.依据AIC、BIC和极大似然比检验化简期望模型；
iii）.用代数差分法将拟合的收获模型转化为代数差分生长模型。
【关键词】线性代数差分模型；
mixed过程；
模型筛选；
林分生长与收获预估
0.前言
在林分生长与收获预估的模型中，差分生长模型得到了广泛的应用。线性差分模型基本上为Schumacher模型的变型，广泛应用于林分蓄积、胸高断面积的建模，以及单位面积株数和优势木树高生长模型。
差分生长模型拟合方法有“直接最小二乘估计法”和“分类变量回归法”[1]。一般认为后者可以获得近似无偏的估计，而前者则导致检验统计量如RMSE的失真[2]。传统上差分生长模型的拟合主要是直接拟合差分生长模型，然后根据拟合统计量如RMSE、R2等确定最优拟合模型。与传统方法不同，本文直接拟合生长模型，在获得参数估计值后，再用代数差分法导出相应的代数差分生长模型。这样做的优越之处在于非常便于对期望模型和方差结构模型进行筛选。更为重要的是，可以通过模型拟合识别最适合的随林分变化参数。
本文详细讨论了如何用“分类变量回归法”和SAS的mixed过程拟合代数差分生长模型，可简述如下：i）.直接以生长收获模型为对象，同时确定一个随林分变化的参数和最优拟合方差结构模型；
ii）. 保持方差结构模型不变，根据拟合统计量逐步化简期望模型；
iii）.在确定最优拟合的期望模型后，运用代数差分法导出相对应的代数差分生长模型。所有拟合与筛选均用SAS的mixed过程完成，并给出了详细的SAS代码和代码解释。
1.方法与材料
1.1数据
数据来源于148个集约经营的火炬松实验人工林逐年观测的固定样地数据（样地约0.152公顷）。SAS的数据集basal的内容如表（1）。
表1 模型拟合的基本数据结构
Table 1 data structure for mode fitting
age=林分年龄；
fert=经营措施，分别取值为H=施加除草剂以控制竞争植物、F=施肥以增加土壤肥力、HF=除草剂和施肥并用、C=对照；
code=样地代码，每个样地有一个唯一代码；
logba=样地胸高断面积的自然对数值；
logdh=样地优势木树高的自然对数值；
iage=林分年龄的倒数；
logtpa=样地株数的自然对数值。共有148个样地数据，1491个记录。
1.2数学模型
所考虑的胸高断面积收获的数学模型为
E（lny）=α+αlnN+αlnH+α+α （1）
这里lny、lnH和lnN分别对应logba、logdh和logtpa；
1/t对应iage，其余为模型参数。通过不同假设，以上模型可以导出很多广泛应用的差分生长模型。例如假设α0为随林分变化的参数，运用代数差分法可以导出Pienaar和Shiver （1986）的胸高断面积模型[3]；
假设α4=0且α0为随林分变化的参数则可导出Forss等人（1996）提出的差分生长模型[4]；
假设α4=0和α2=0且α3为随林分变化的参数，则可以导出Souter （1986）的模型[5]。
1.3参数估计方法
从以上代数差分法的讨论中，总是假设模型（1）中有一个参数随林分变化而变化，因而使用“分类变量回归法”以便考虑这一特点。如果用SAS的reg过程，必须构造特殊的数据结构。用SAS的reg过程拟合线性差分生长模型有两点不足之处：1）.需要构造复杂的数据文件。以本研究的为例，数据文件中需要增加额外的148个变量；
2）.reg过程无法考虑重复观测数据的自相关性和异质方差结构。由于所考虑的胸高断面积模型均为线性模型，而且数据为典型的重复观测数据，因而模型的拟合与筛选均用SAS的proc mixed过程。
1.4模型筛选和筛选指标选用
模型筛选包括期望模型和最优拟合的自相关与异质方差结构模型的筛选。Ngo和Brand（1997）推荐了两种模型选择方法[6]。一种方法就是首先列出所有可用的数学期望模型和所有可能的方差结构模型，拟合二者间的所有组合，然后根据模型拟合指标选择最优拟合模型。另一种方法为Wolfinger和Diggle提出的方法，即首先考虑最复杂的期望模型，并保持期望模型不变，然后选择最优拟合的方差结构模型。一旦选出最优拟合的方差模型，保持方差模型不变，再逐步简化期望模型。模型筛选指标较多，但常用的指标为极大似然比检验（LRT）、AIC、BIC（或称为SBC）。LRT用于嵌套的模型，而AIC和BIC用于非嵌套模型。
2.结果与分析
2.1残差的方差结构模型
为考虑经营措施对胸高断面积生长过程的影响，将模型（1）改写如下：
这里参数?为不同经营措施对截距和各项回归系数的影响，下标k分别取值为C、F、H和HF，对应四种经营措施；
由于差分生长模型需要确定一个随林分变化而变化的参数，因而模型拟合的首要任务就是识别该参数。假设“1/t”的回归系数随林分变化而变化，那么模型（2）则可以表示如下：
这里βi为随林分变化的参数。
根据以往的研究，所考虑的方差模型为：CS模型（同一林分不同年龄的观测值方差相同，而且相关系数为常数）、AR（1）模型（不同年龄观测值的方差相同，但相关系数为一阶自相关模型）、ARH（1）模型（相关系数为一阶自相关模型， 但不同年龄观测值的方差相异）和ARMA（1，1） 模型（不同年龄观测值的方差相同，但相关系数为一阶自相关移动平均模型）。
由于模型（2）中共有5个参数可以考虑为随林分而变化的参数（b0～b4），因而相对应地共有5个期望模型，加上所考虑的4个方差结构模型，所有组合共计为20个模型。首先拟合这20个模型，根据拟合统计量AIC、BIC确定随林分变化的参数及最优拟合的方差结构模型，这部分与Ngo和Brand的方法一致；
然后保持方差结构模型不变，根据LRT、AIC、BIC剔除回归效果不显著的因子，包括经营措施效果等，即化简期望模型。而这一部分则与Wolfinger和Diggle的方法一致。方差结构筛选的基本SAS代码如下：
以上代码及以后代码中的粗体字为SAS系统关键词。class语句声明了3个分类变量，分别为fert、code及age。model语句中既包含离散型的分类变量，又包含连续型变量（如logba、logtp等，即未在class语句中声明的变量）。以上代码的model语句意义如下：
i）.代码中model语句中的分类变量表示该变量对回归模型截距的影响。例如model语句中fert表示不同经营措施下的logba的回归模型的截距不同，随fert变量的取值变化而变化。
ii）.如果一个连续型变量与一个分类变量的乘积出现model语句中，且仅以该形式出现，如code*iage，则表示该连续型变量的回归系数随分类变量的取值变化而变化。code是一个分类变量，每个取值对应一个样地，而iage为一个连续型变量，那么code*iage则表示每个样地的iage项的回归系数均不相同，各有其取值。
iii）.如果一个连续型变量与一个分类变量的乘积出现，而且该连续型变量同时独立出现在模型中，同样表示该连续型变量的回归系数随分类变量的取值变化而变化，但与ii）的意义有所不同。例如logdh和fert*logdh同时出现在model语句中，表示logdh的回归系数随fert取值的变化而变化，而且每个fert取值的回归系数均表示为两个分量之和。一个分量为参照回归系数，另一个分量为与参照回归系数相比的增量（可负可正）。mixed过程按照英文字母的排列顺序，以fert的最后一个取值（即HF）为参照，其回归系数即为参照回归系数。而fert其它取值的回归系数则表示为参照回归系数（即HF）和与之相比增量之和。
以上model语句的对应的数学模型如下：
以上四个模型分别对应经营措施C、F、H和HF的回归模型。参数α0、α1、α2、α4为在HF经营措施下，林分胸高断面积的回归系数（参照回归系数）；
参数φi，k为k经营措施的第i个回归系数与HF相对应的回归系数相比之增量（k取值为C、H和F）。以对照C的lnH的系数为例，其回归系数为α2+φ2，C。这里φ2，C就是C的lnH的回归系数与HF的lnH的回归系数相比的增量。以上代码中，fert*logtpa、fert*logdh及fert*logdh*iage分别对应各项回归系数与HF的回归系数相比的增量，即φ；
logtpa、logdh及logdh*iage则对应参照回归系数，即HF的回归系数α0、α1、α2和α4等。
依据模型拟合统计量AIC和BIC值（其值越小，拟合效果越好）选择模型，指定1/t的回归系数为随林分变化的参数和arh（1）为方差结构模型时，AIC和BIC均取得了最小值，分别为-3629.9和-3585.0。因而arh（1）为最优拟合的方差结构模型，并且可以确认1/t的回归系数为随林分变化的参数。
2.2最优拟合的期望模型
保持arh（1）结构不变，根据AIC、BIC及极大似然比检验的p值进行期望模型的筛选。模型（3）的拟合结果表明，经营措施H和HF各个参数差异并不显著，φ0，H、φ1，H、φ2，H 和φ4，H是否为0的t检验p值均大于0.05，说明经营措施H的各项回归系数与HF并无显著差别。将二者的回归系数合并为一，仅保留截距不同，可拟合如下模型：
这里bi，k=αi，k+φi，k（i=0，1，2，4；
k=C，F，H和HF）；
b1，p、b2，p、b4，p为经营措施H和HF的回归系数，b0，H和b0，HF则分别为二者的截距。模型（3）便于比较不同经营措施的效果异同，而模型（4）则便于化简期望模型，即判断一个回归因子对胸高断面积的影响是否显著。注意模型（3）和（4）的唯一区别仅在于将模型（3）中第3个和第4个模型合并为一个而已（仅截距不同）。为拟合以上模型，新引进了一个分类变量“f”，其取值分别为C、F和H，其中经营措施H和HF的样地“f”变量值均为H。拟合模型（4）的SAS代码如下：
模型（4）的AIC、BIC的值均小于模型（3）的AIC和BIC值，LRT检验p值为0.1573，表明模型（4）与模型（3）拟合效果基本相当，但模型（4）更为简洁。由于参数的t检验结果表明仅对照C的lnH/t参数与0相比差异显著，因而剔除F、H和HF的lnH/t项（相应的模型称为模型5）。相对应地，在SAS的数据文件中引入了一个新变量“check”，当且仅当林分的经营措施为C取值为1，其它情形均取值为0。对应的SAS代码如下：
与模型（4）的SAS代码相比，仅用check*logdh*iage 取代了f*logdh*iage。模型（5）的参数t检验结果表明H和HF的回归模型中lnH的回归系数b2，p与0比较差异不显著（p值为0.1181），但剔除该项后的模型拟合的LRT检验平值、AIC、AICC、BIC表现出一种相互矛盾的结果。AIC和AICC表明两个模型拟合效果基本一致，而BIC和LRT检验结果则倾向于选择化简后的模型。本研究最终以模型（5）为最优拟合模型，部分SAS输出的参数估计值见表（2）。
表2 模型（5）的参数估计值
Table 2 parameters estimates for model （5）（下转第164页）
（上接第109页）*仅列出参数估计值部分，与code值相对应148参数估计值未列于此表。
2.3代数差分生长模型的导出
根据表（2）的模型结构，HF对应的回归模型为：
ln=+lnN+lnH+β （6）
假设t=t1时ln1已知，求解关于βi的表达式可有βi=t1（ln--lnN-lnH）；
代入模型（6）中可得当t=t2时的lny2的估计式，即：
ln=lny+++
其中ln1用lny1代替。其它经营措施的代数差分模型的导出基本相似，在此略过。
3.讨论与结论
本研究的模型筛选方法与以往研究不同。以往对差分生长模型筛选方法为：拟合选取的差分生长模型，让后根据拟合统计量如RMSE、R2选择最优拟合模型。本文从生长收获模型入手，而不是从差分生长模型入手。首先确定最优拟合的方差结构模型，同时确定最适合作为随林分变化的参数。然后保持方差结构模型不变的前提下，通过AIC、BIC及极大似然比检验的p值逐步化简期望模型。所有模型拟合和筛选均用SAS的proc mixed过程。该过程不仅非常便于拟合有分类变量的回归模型（例如经营措施），而且可以非常方便地拟合方差结构模型。本文所讨论的差分生长模型的拟合方法、筛选方法和拟合过程，均可以应用于其它林分因子如蓄积、优势木树高的线性差分生长模型。对于非线性模型，除了需要构造特殊的数据文件外，其它内容也基本相同。
[科]
【参考文献】
[1]倪成才，刘春梅，丁俊峰，潘晓茹.差分生长模型预测误差的分析[J].北京林业大学学报，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.
