林业科学  2016, Vol. 52 Issue (10): 72-79   PDF    
DOI: 10.11707/j.1001-7488.20161009
0

文章信息

祖笑锋, 李秋实, 倪成才, 覃先林, NighGorden
Zu Xiaofeng, Li Qiushi, Ni Chengcai, Qin Xianlin, Nigh Gorden
非线性混合效应生长模型的拟合、随机效应预测和应变量预测间对应关系
Analysis and Comparison of Combinations among Fitting NLME and Predictors of Random Parameters and Response Variables
林业科学, 2016, 52(10): 72-79
Scientia Silvae Sinicae, 2016, 52(10): 72-79.
DOI: 10.11707/j.1001-7488.20161009

文章历史

收稿日期:2015-07-07
修回日期:2015-11-23

作者相关文章

祖笑锋
李秋实
倪成才
覃先林
NighGorden

非线性混合效应生长模型的拟合、随机效应预测和应变量预测间对应关系
祖笑锋1,2, 李秋实2, 倪成才2 , 覃先林1, NighGorden3    
1. 中国林业科学研究院资源信息研究所 北京 100091;
2. 北华大学林学院 吉林 132013;
3. British Columbia Ministry of Forests, Lands and Natural Resources Operations, Forest Analysis and Inventory Branch, P.O.BOX9512, Stn.Prov.Govt.Victoria, B.C.V8W 9C2, Canada
摘要【目的】 在非线性混合效应模型的拟合、随机效应预测和应变量预测3个环节上,常用一阶泰勒近似法将非线性模型线性化,泰勒近似的基点一般为随机效应参数的数学期望或迭代终值。在林业中,3个环节的基点常常并不完全一致,并可能影响预测精度。本研究以树高生长过程为例,分析了不一致的基点对预测精度的影响程度。 【方法】 以加拿大哥伦比亚省美国黄松解析木数据为基础,以三参数Logistic为基本模型,随机抽取49株用于拟合,30株用于验证。采用R语言的nlme函数和SAS的nlmixed过程拟合模型。nlme函数以随机效应的数学期望为基点,而nlmixed过程则以迭代终值为基点。利用SAS的IML过程预测随机效应与应变量,并计算预测精度。以预测误差均方(MSPE)、平均相对误差(MPE)、平均相对误差绝对值(MAPE)作为评价预测精度的指标。 【结果】 预测随机效应与预测应变量的基点不同,预测精度将大幅度下降。 【结论】 预测随机效应和预测应变量的基点必须一致,且仅需二者一致,与模型拟合的基点基本无关。如果二者一致,以数学期望或迭代终值为基点对预测精度基本上无显著影响;如果预测随机效应和预测应变量的基点不同,将显著降低预测精度。
关键词: 非线性混合效应模型     生长与收获预测     Logistic方程    
Analysis and Comparison of Combinations among Fitting NLME and Predictors of Random Parameters and Response Variables
Zu Xiaofeng1,2, Li Qiushi2, Ni Chengcai2 , Qin Xianlin1, Nigh Gorden3    
1. Research Institute of Forest Resource Information Techniques, CAF Beijing 100091 ;
2. College of Forestry, Beihua University Jilin 132013 ;
3. British Columbia Ministry of Forests, Lands and Natural Resources Operations, Forest Analysis and Inventory Branch, P. O. BOX9512, Stn. Prov. Govt. Victoria, B. C. V8W 9C2, Canada
Abstract: 【Objective】 Fitting non-linear mixed effects models (NLME), predicting the random effects parameters, as well as predicting the response variable, often involve a Taylor series expansion for linearization, based upon either the expected value of the random effects or final iterative value. In forestry, however, the linearization bases are not always consistent as they should be, and probably reduced the accuracy of prediction. In this paper, we investigated the tree height growth and discussed the effects of inconsistency among the linearization bases for fitting, predicting random effects the response. 【Method】 We randomly selected 49 trees for NLME-fitting and 30 trees for validation from 79 dominant trees of ponderosa pine in British Columbia, Canada. The base model was three-parameter Logistic. We used the nlme function in R and the nlmixed procedure in SAS for model fitting, respectively corresponding linearization based upon the expected value and the final iterative value. The IML procedure in SAS was employed for predicting the random effects and the response. Mean squared prediction error (MSPE), mean percentage error (MPE), and mean absolute percentage error (MAPE) were used as evaluation criteria. 【Result】 The results showed that inconsistent linearization bases between the random effects and the response significantly decreased the accuracy of the response prediction. 【Conclusion】 The linearization bases between the random effects and the response had to be consistent, and enough for obtaining predictions as accurate as possible. The accuracy of prediction was invariant to the linearization base for model-fitting and to either the expected value or the final iterative value, which was used as a linearization base.
Key words: non-linear mixed effects(NLME)     growth and yield prediction     Logistic function    

非线性混合效应模型(non-linear mixed effects, NLME)中,包括固定效应参数和随机效应参数,适于分析重复观测数据和纵向数据。由于NLME模型参数可以随研究对象变化而随机变化,在林业中得到了越来越广泛的应用(Lappi et al., 1988Fang et al., 2001Calegario et al., 2005李永慈等,2004李春明等,2010姜立春等,2011符利勇等,20122013)。

拟合非线性混合效应模型时,一般需要采用一阶泰勒级数将模型线性化,然后再用迭代算法确定极大似然函数并求解(Davidian et al., 19931995)。在线性化时,既可以用随机效应参数的数学期望为线性化的基点,亦即0值点,也可以用随机效应参数在迭代过程中的迭代值为基点。在常用的统计软件中,R和S+以0值为基点,而SAS软件以迭代终值为基点。研究表明,以迭代终值为基点的参数估计值优于以0值为基点,但需要更多的计算时间,而且往往不稳定(Lindstrom et al., 1990Wolfinger et al., 1997)。

混合效应参数模型可用于分析不同来源的随机效应参数对应变量的影响程度;且在获得一个研究对象的若干观测值后,混合效应参数模型更可以用于预测。采用非线性混合效应参数模型进行预测时,需要2步:根据已经取得的观测值向量预测随机效应参数值;在确定了一个研究对象的随机效应值后,可以进一步预测应变量的其他观测值。随机效应参数的预测方法一般称之为经验线性无偏最优估计,简称为EBLUP(empirical best linear unbiased predictor)。与拟合类似,EBLUP同样涉及利用一阶泰勒近似将非线性模型线性化的过程,线性化基点同样可为0值或迭代终值。线性化基点不同,也导致了不同的应变量预测方法。

基于线性化基点角度,从模型拟合、随机效应参数预测到应变量预测,形成了多种组合。本文重点研究各种组合对应变量预测精度的影响。

1 材料与方法 1.1 数据

本文数据来源于加拿大哥伦比亚省(British Columbia)美国黄松(Pinus ponderosa)林分的临时样地数据。哥伦比亚省林业局于2000和2001年建立了100块临时样地,样地设置时充分考虑了现有林分的立地质量、年龄、密度以及物种丰富度等因素,本文从中选取具有代表性的林分作为样地。样地形状为圆形,面积为0.01 hm-2(半径5.64 m)。每块样地挑选1株未受损或未受压的健康优势木或亚优势木作为样木,最终得到79株优势木的解析木。样木伐倒后,对于树桩和树干部分采取纵切的方法,由髓心直接查出树高-年龄值;对于不适合做纵切的树梢部分采取轮生枝法确定每年的树高生长量,可以精确得到每年的树高值。这种获得树高年龄数据的方法,其精确性远远高于截取圆盘进行树干解析的方法。本文随机选取49株优势木进行拟合,30株优势木进行预测,数据的详细介绍见文献(Nigh,2004)。

1.2 树高生长模型的选择与建立

林分优势树高生长模型的种类很多,主要有Logistic,Mitscherlich,Gompertz,Weibull,Richards和Korf方程等,但广泛应用的有Richards-Chapman,Logistic,Weibull和Korf方程,其模型结构见表 1

表 1 优势木树高生长方程和混合效应模型 Tab.1 The equations of dominant tree and mixed-effects model
1.3 随机效应参数预测和观测值预测

非线性混合参数模型可一般性地表示如式(1):

$ {{y}_{i}}=f\left( {{A}_{i}}\beta +{{B}_{i}}{{b}_{i}},{{x}_{i}} \right)+{{e}_{i}}。 $ (1)

式中:yi为第i个研究对象的应变量向量;βp维固定效应参数向量;biq维随机效应参数向量,且有bi~N(0, D),D为随机效应间的协方差矩阵;AiBi为具有适当维数的关联矩阵(即由0或1构成的矩阵);xi为自变量矩阵;ei为与yi相关联的随机误差项向量,ei~N(0, Ri),Ri为随机误差项间的协方差矩阵。一般假设biei相互独立。模型中需要估计的参数分别为βDRi。需要注意,Ri一般可用方差函数和相关系数函数表示。

bi=b0为基点,可将式(1)通过一阶泰勒近似表示为式(2):

$ {{y}_{i}}\approx f\left( {{A}_{i}}\beta +{{B}_{i}}{{b}_{0}},{{x}_{i}} \right)+{{Z}_{i}}\left( {{b}_{i}}-{{b}_{0}} \right)+{{e}_{i}}。 $ (2)

式中:Zi为式(1)关于bi的偏导数向量,为自变量xi的函数,形式如式(3):

$ {{Z}_{i}}=\frac{\partial f\left( {{A}_{i}}\beta +{{B}_{i}}{{b}_{i}},{{x}_{i}} \right)}{\partial {{b}_{i}}}\left| _{\left( {{b}_{i}}={{b}_{0}} \right)} \right.。 $ (3)

式(3)中随机变量包括biei 2个向量以及二者的线性函数yi。由于biei相互独立,从式(3)可有:

$ \begin{align} & \text{Cov}\left( {{y}_{i}} \right)={{Z}_{i}}DZ_{i}^{\text{t}}+{{R}_{i}}; \\ & \ \ \ \ \text{Cov}\left( {{b}_{i}},{{y}_{i}} \right)=DZ_{i}^{\text{t}}。 \\ \end{align} $

上式中上标t为矩阵转置运算符。随机向量biyi的联合分布为:

$ \begin{array}{l} \left( \begin{array}{l} {b_i}\\ {y_i} \end{array} \right) \sim N\left[ {\left( {\mathop {f\left( {{A_i}\beta + {B_i}{b_0},{x_i}} \right) - {Z_i}{b_0}}\limits^0 } \right)} \right.,\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left( {\begin{array}{*{20}{c}} D&{DZ_i^{\rm{t}}}\\ {{Z_i}{D^{\rm{t}}}}&{{Z_i}DZ_i^{\rm{t}}{R_i}} \end{array}} \right)} \right]. \end{array} $ (4)

对于总体中的第i个研究对象,有2个应变量向量yi(1)yi(2),分别对应于xi(1)xi(2)。如下将以yi(1)为基础求出随机效应bi的预测值${{\hat{b}}_{i}}$,然后再以${{\hat{b}}_{i}}$xi(2)为基础,给出yi(2)的预测值${{\hat{y}}_{i\left(2 \right)}}$

根据多元正态分布理论,当yi=yi(1)时,从式(4)可以得出bi的条件分布数学期望为式(5),其中E(bi)=0。统计上,一般以数学期望预测随机变量的取值,因而随机变量biyi=yi(1)时的预测值可表示为式(6):

$ \begin{align} & \text{E}\left( {{b}_{i}}\left| _{{{y}_{i}}={{y}_{i}}\left( 1 \right)} \right. \right)=\text{E}\left( {{b}_{i}} \right)+DZ_{i\left( 1 \right)}^{\text{t}}{{\left( {{Z}_{i\left( 1 \right)}}DZ_{i\left( 1 \right)}^{\text{t}}+{{R}_{i}} \right)}^{-1}} \\ & \ \ \ \ \ \ \ \ \ \left[ {{y}_{i\left( 1 \right)}}-f\left( {{A}_{i}}\beta +{{B}_{i}}{{b}_{0}},{{x}_{i\left( 1 \right)}} \right)+{{Z}_{i\left( 1 \right)}}{{b}_{0}} \right]; \\ \end{align} $ (5)
$ \begin{align} & {{{\hat{b}}}_{i}}=DZ_{i\left( 1 \right)}^{\text{t}}{{\left( {{Z}_{i\left( 1 \right)}}DZ_{i\left( 1 \right)}^{\text{t}}+{{R}_{i}} \right)}^{-1}}\left[ {{y}_{i\left( 1 \right)}}- \right. \\ & \ \ \ \ \ \ \left. f\left( {{A}_{i}}\beta +{{B}_{i}}{{b}_{0}},{{x}_{i\left( 1 \right)}} \right)+{{Z}_{i\left( 1 \right)}}{{b}_{0}} \right]。 \\ \end{align} $ (6)

由式(6)可知,在已获得βDRi估计值的前提下,给定一个基点b0值,就可以求出相应的${{\hat{b}}_{i}}$。不同的b0值会在一定程度上导致${{\hat{b}}_{i}}$值的变化。一阶泰勒近似精度与基点选择有关,如果以式(6)求出的${{\hat{b}}_{i}}$值作为基点,则利用式(6)可进一步改进bi估计值。反复利用式(6),可形成逐步改进${{\hat{b}}_{i}}$值的迭代过程。以第k次迭代的bi预测值${{\hat{b}}_{i\left(k \right)}}$为基点(k=0, 1, 2,…),第k+1次迭代的bi预测值${{\hat{b}}_{i\left(k+1 \right)}}$则可表示为式(7)。${{\hat{Z}}_{i\left(1 \right)}}$可通过令式(3)中的xi=xi(1)β=$\hat{\beta }$求出。

$ \begin{align} & {{{\hat{b}}}_{i\left( k+1 \right)}}=\hat{D}\hat{Z}_{i\left( 1 \right)}^{\text{t}}{{\left( {{{\hat{Z}}}_{i\left( 1 \right)}}\hat{D}\hat{Z}_{i\left( 1 \right)}^{\text{t}}+{{{\hat{R}}}_{i}} \right)}^{-1}} \\ & \left[ {{y}_{i\left( 1 \right)}}-f\left( {{A}_{i}}\hat{\beta }+{{B}_{i}}{{b}_{i\left( k \right)}},{{x}_{i\left( 1 \right)}} \right)+{{{\hat{Z}}}_{i\left( 1 \right)}}{{b}_{i\left( k \right)}} \right]。 \\ \end{align} $ (7)

$\left| {{{\hat b}_{i\left( {k + 1} \right)}} - {{\hat b}_{i\left( {k + 1} \right)}}} \right| < \varepsilon $时,即当${{{\hat{b}}}_{i\left(k \right)}}$${{{\hat{b}}}_{i\left(k+1 \right)}}$的差值小于预先指定的很小正数ε时,二者近似相等,表示已无改进余地,此时式(6)中的b0=${{{\hat{b}}}_{i\left(k \right)}}$bi=${{{\hat{b}}}_{i\left(k+1 \right)}}$。令b0=${{{\hat{b}}}_{i\left(k \right)}}$bi=${{{\hat{b}}}_{i\left(k+1 \right)}}$xi=xi(2),可由式(2)得出${{{\hat{y}}}_{i\left(2 \right)}}$,表示如式(8)。忽略${{{\hat{b}}}_{i\left(k \right)}}$${{{\hat{b}}}_{i\left(k+1 \right)}}$间的差异,式(8)中的${{Z}_{i}}\left({{{\hat{b}}}_{i\left(k+1 \right)}}-{{{\hat{b}}}_{i\left(k \right)}} \right)\approx 0$,并将${{{\hat{b}}}_{i\left(k+1 \right)}}$一般性地表示为${{{\hat{b}}}_{i}}$,式(8)可表示为式(9):

$ \begin{align} & {{{\hat{y}}}_{i\left( 2 \right)}}\approx f\left( {{A}_{i}}\hat{\beta }+{{B}_{i}}{{b}_{i\left( k \right)}},{{x}_{i\left( 2 \right)}} \right)+ \\ & \ \ \ \ \ \ {{{\hat{Z}}}_{i\left( 2 \right)}}\left( {{{\hat{b}}}_{i\left( k+1 \right)}}-{{{\hat{b}}}_{i\left( k \right)}} \right); \\ \end{align} $ (8)
$ {{{\hat{y}}}_{i\left( 2 \right)}}\approx f\left( {{A}_{i}}\hat{\beta }+{{B}_{i}}{{b}_{i\left( k \right)}},{{x}_{i\left( 2 \right)}} \right)。 $ (9)

在迭代过程中,一般基点的迭代初值取0。如果不继续迭代,那么b0=0,式(2)和(6)则分别变为式(10)和(11)。与式(10)相对应地将${{{\hat{y}}}_{i\left(2 \right)}}$表达为式(12):

$ {{y}_{i}}\approx f\left( {{A}_{i}}\beta ,{{x}_{i\left( 1 \right)}} \right)+{{Z}_{i}}{{b}_{i}}+{{e}_{i}}; $ (10)
$ \begin{align} & {{{\hat{b}}}_{i}}=\hat{D}\hat{Z}_{i\left( 1 \right)}^{\text{t}}{{\left( {{{\hat{Z}}}_{i\left( 1 \right)}}\hat{D}\hat{Z}_{i\left( 1 \right)}^{\text{t}}+{{{\hat{R}}}_{i}} \right)}^{-1}} \\ & \ \ \ \ \ \left[ {{y}_{i\left( 1 \right)}}-f\left( {{A}_{i}}\hat{\beta },{{x}_{i\left( 1 \right)}} \right) \right]; \\ \end{align} $ (11)
$ {{{\hat{y}}}_{i\left( 2 \right)}}\approx f\left( {{A}_{i}}\hat{\beta },{{x}_{i\left( 2 \right)}} \right)+{{{\hat{Z}}}_{i\left( 2 \right)}}{{{\hat{b}}}_{i}}。 $ (12)

给定一个研究对象,如果已经获得了观测值向量yi(1),则可以用式(7)通过迭代求出bi的预测值,也可以用式(11)直接求出bi的预测值${{{\hat{b}}}_{i}}$,区别在于一阶泰勒近似的基点不同。获取${{{\hat{b}}}_{i}}$值后,对yi(2)的预测可用式(9)或(12)。从以上分析可知,式(7)和(9)自然地联系起来,而式(11)和(12)也同样如此。

式(9)对于${{{\hat{b}}}_{i}}$而言是非线性的,而式(12)则为线性的。为后续讨论描述简洁,将式(9)和(12)分别称为非线性预测和线性预测。在林业中,用式(11)预测随机效应的应用更多一些(Fang et al., 2001Hall et al., 2001Calama et al., 2004Sharma et al., 20072009Temesgen et al., 2008Huang et al., 2009bJiang et al., 2010祖笑锋等,2015),利用式(7)的也有(Hall et al., 2004Meng et al., 2009Huang et al., 2009aNi et al., 2012Yang et al., 2013)。

从模型拟合、随机效应参数预测和应变量yi(2)的预测,每个环节都包含了2种方法,共计形成8种组合,本文重点研究各种组合对应变量预测精度的影响。当基点b0=b^ib^i接近于bi的真值时,由于泰勒近似在基点附近的近似效果较好,式(2)对模型(1)的近似效果要优于式(10)。同样由于近似基点不同的原因,相对应地,式(7)对bi的预测效果也应该优于式(11)。直观上,均以迭代终值为基点的模型拟合(SAS的nlmixed过程)、随机效应预测[式(7)]、相关联的应变量预测[式(9)],此组合应优于其他任何组合。事实是否如此,为本文的研究目的之一。

虽然式(7)与(9)之间以及式(11)和(12)之间存在承接关系,但在林业上却多采用式(11)预测随机效应参数,采用式(9)预测yi(2)(Fang et al., 2001Sharma et al., 20072009Jiang et al., 2010)。此种组合对预测精度影响如何,也是本文的研究目的之一。

本文采用3种预测精度评价指标,其定义如表 2。3种评价指标值越小,表明预测精度越高,其最优值都为0。

表 2 预测精度的评价指标 Tab.2 Evaluation indexes for prediction accuracy
2 结果与分析 2.1 树高生长模型的比较

采用以上4个生长模型作为基本模型,利用R软件的nlme过程对49株解析木进行拟合,其统计量见表 3

表 3 拟合统计量 Tab.3 Fit statistics

使用AIC,BIC和Loglik 3个模型评价指标对拟合结果进行评价。由表 3可知,Logistic方程的AIC,BIC最小且logLik最大,表明Logistic方程拟合效果最好,因此将其作为混合效应模型的基本模型。

2.2 拟合非线性混合效应模型

在林业中,广泛应用Logistic方程分析优势木树高生长与年龄的关系。采用三参数Logistic非线性混合效应模型具有相当大的灵活性,Logistic模型如式(13):

$ {{y}_{i,j}}=\frac{{{\beta }_{1}}}{1+{{e}^{{{\beta }_{2}}+{{\beta }_{3}}\ln \left( {{x}_{i,j}} \right)}}}+{{e}_{i,j}}。 $ (13)

若将模型中的3个参数都处理为混合参数,如式(14)(Calegario et al., 2005):

$ {{y}_{i,j}}=\frac{{{\beta }_{1}}+{{b}_{1,i}}}{1+{{e}^{\left( {{\beta }_{2}}+{{b}_{2,i}} \right)+\left( {{\beta }_{3}}+{{b}_{3,i}} \right)\ln \left( {{x}_{i,j}} \right)}}}+{{e}_{i,j}}。 $ (14)

式中:yi, j表示第i个林分在年龄xi, j时的树高;β1, β2, β3为方程的固定效应参数;与之相对应的随机效应参数分别为bi, 1bi, 2bi, 3ei, j为方程误差项,且与bi, 1bi, 2bi, 3任何一个均相互独立。

随机抽取49株解析木用于模型拟合,采用R语言的nlme函数和SAS的nlmixed过程进行拟合,分别对应于以0值为基点和以迭代终值为基点的拟合方法。参数估计值如表 4所示,其中${{{\hat{\beta }}}_{1}}$, ${{{\hat{\beta }}}_{2}}$, ${{{\hat{\beta }}}_{3}}$为固定效应参数估计值;随机效应参数协方差矩阵估计值${\hat{D}}$为三维对称矩阵;${{{\hat{\sigma }}}^{2}}$为随机误差项估计值。以上参数主要用于计算新个体的随机效应参数值。

表 4 模型拟合结果 Tab.4 Model fitting results
2.3 预测精度分析

本文对30株验证数据设置了不同的起始年龄和间隔,从多角度评价预测优势木树高生长精度。在表 5中,字母组合代码的第1位代表预测随机效应参数的线性化方法,Z值表示以期望值(0值)为线性化基点,而E值则代表以迭代终值为基点;第2位代表应变量预测方法,N表示非线性预测,L表示线性预测。表中的年龄值为观测值yi(1)对应的年龄值,预测应变量yi(2)为其后面连续5年的观测值。如以第(1)行中30, 31, 32, 33, 34作为观测值,预测35, 36, 37, 38, 39年的优势木树高值。8种组合均计算了MSPE,MPE以及MAPE。

表 5 不同参数估计方法时预测随机效应bi和因变量各种组合值 Tab.5 All combinations of predictors for predictors for bi and the response

从列于表 5的结果可见,无论模型拟合以0值为基点还是以迭代终值为基点,ZL组合在2种模型拟合方法中的表现大致相同,EN组合也是如此,与模型拟合方法无关,ZL和EN组合总是明显地优于ZN和EL组合。根据表 5中的组合代码定义,EN组合对应于式(7)与(9),而ZL组合对应于式(11)和(12)。表 5预测精度表明,随机效应预测方法决定了应变量的预测方法。拟合模型(10)对应于用模型(11)计算bi值,在模型预测时自然应该以模型(10)为基础,对应的模型预测公式则为式(12);同理,拟合模型(2)对应于用模型(7)计算bi值,自然以式(2)为基础,并由式(2)进一步导向了式(9)。因而,式(7)与(9)以及式(11)与(12)之间必须配合使用,而与模型拟合方法的关系不大。

2.4 单株解析木示例 2.4.1 随机效应参数预测分析

为分析基于0值和迭代终值的2种线性化方法对随机效应参数预测的影响,对同一解析木,在xi(1)相同的前提下,分别用2种方法获取随机效应参数预测值,共计有30株解析木的30对预测值,二者对比结果如图 1所示。基于0值和迭代终值的2种线性化方法预测随机效应参数的结果基本上一致,说明式(2)和(10)对模型(1)的近似效果对大多数解析木而言非常接近,表明多数解析木的bi位于其期望值附近。

图 1 单株bi参数估计值对比分析 Fig.1 Comparisons of parameters estimates of bi
2.4.2 应变量预测结果分析

为分析随机效应和应变量预测间各种组合对预测结果的影响,本文选取1株解析木,用4种不同组合进行匹配验证,结果如图 2所示。从图 2可见,ZL和EN组合几乎完全重合,并且均能紧密跟踪数据散点;与此相反,ZN和EL组合则均明显地偏离数据散点。

图 2 4种不同匹配方法预测优势木树高 Fig.2 Predicting dominant height from the four combinations of predictors

以ZN和ZL组合为例,由于模拟拟合和随机效应参数的预测方法完全一致,二者的所有参数估计值均相同,唯有预测yi(2)的模型不同,分别为式(9)和(12)。ZN和ZL组合曲线间的差异,完全由式(9)和(12)间差异所致。

ZN组合,亦即式(13)与(11)搭配使用,在林业中误用较多。由表 5可知,ZN组合与ZL和EN组合相比,其预测精度急剧下降;且ZL组合预测精度略微优于EN组合,但尚不足以得出ZL优于EN的结论,可认为二者预测精度处于同一量级水平。但从应用方便角度看,式(9)有时无法收敛,而不能求出bi值,因而ZL组合优于EN组合。

表 5中,第(1)行观测间隔为1年,而(2)~(10)行观测间隔为2年,其目的是预测不同年龄段内MSPE,MPE以及MAPE的变化。对于正确的匹配方法,MSPE,MPE和MAPE 3种评价指标值表现得都很稳定,且明显小于错误匹配。但特别强调的是第(11)~(14)行,观测数目相同,观测间隔不同,随着观测间隔增加,3种评价指标的预测精度均有所提升。

3 结论与讨论

本文从线性化基点角度出发,分析了NLME模型拟合、随机效应参数预测以及应变量预测间多种线性化基点组合对应变量预测精度的影响,并以MSPE,MPE和MAPE 3种指标作为预测精度评价指标进行讨论。从所用数据和模型角度看,经验性分析结果表明,模型拟合与随机效应参数预测方法间线性化基点一致与否对预测精度的影响非常有限;与此相反,如果随机效应参数预测和应变量预测间线性化基点不同,将极为显著地降低应变量的预测精度。

混合效应模型拟合和随机效应参数预测可以视为2个不同的拟合阶段,分别对于总体特征和1个特定个体。混合效应参数模型拟合为第1阶段,其目的在于获取总体特征参数,如总体参数βD以及与Ri相关联的方差函数和自相关函数的参数估计值。结果表明,以迭代终值为基点优于以0值为基点的参数估计值(Lindstrom et al., 1990Wolfinger et al., 1997)。第2阶段拟合,即在已经获得总体参数估计值的前提下,对1个特定研究对象进行再次拟合,获取其对应的随机效应参数bi值,并直接应用于此特定个体的应变量预测。除总体参数β估计值外,DRi的估计值仅通过影响bi值而间接影响应变量的预测。结果表明,这种总体参数对应变量预测的间接作用,在第2阶段拟合时已极大的弱化。

经验性分析结果表明,模型拟合、随机效应参数预测和应变量预测三者间线性化基点完全一致,与仅仅后二者相一致的预测精度相比,并无任何明显优势(参考表 5);与此相反,如果模型拟合和随机效应预测的线性化基点一致,同一个体的所有参数估计值均一致,但在预测应变量时,不同基点的预测精度却完全不同(参考图 2)。这一结果表明,线性化基点不同有时导致了近似程度差异较大,二者间的差异程度取决于bi值与期望值的差异程度。这种差异存在于模拟拟合、随机效应预测和应变量预测中任何二者之间,但对应变量预测精度的影响却截然不同。随机效应参数预测和应变量预测间线性化基点不一致,更为直接和显著地影响应变量预测精度,而模型拟合和随机效应间的不一致,则基本上微乎其微。

从应变量预测角度看,随机效应参数预测和应变量预测间线性化基点一致更为重要。在预测随机效应时,虽然以迭代终值为基点可能获取更为精确的预测值,但在迭代运算时有时无法收敛,而以0点为基点总能给出随机效应的预测值,因而建议使用0值为基点[(式11)]预测随机效应,并以与之相对应的应变量线性预测[(式12)]进行应变量预测。在林业中,以0值为基点预测随机效应却配以不相应的非线性预测[(式12)]的组合应用较多,对此组合对应变量预测精度的影响,尚未引起足够的重视。

参考文献(References)
[] 符利勇, 张会儒, 唐守正. 2012. 基于非线性混合模型的杉木优势木平均高. 林业科学 , 48 (7) : 66–71.
( Fu L Y, Zhang H R, Tang S Z.2012. Dominant height for Chinese fir plantation using nonlinear mixed effects model based on linearization algorithm. Scientia Silvae Sinicae , 48 (7) : 66–71. [in Chinese] )
[] 符利勇, 孙华. 2013. 基于混合效应模型的杉木单木冠幅预测模型. 林业科学 , 49 (8) : 65–74.
( Fu L Y, Sun H.2013. Individual crown diameter prediction for Cunninghamia lanceolata forests based on mixed effects models. Scientia Silvae Sinicae , 49 (8) : 65–74. [in Chinese] )
[] 姜立春, 刘瑞龙. 2011. 基于非线性混合模型的落叶松树干削度模型. 林业科学 , 47 (4) : 101–106.
( Jiang L C, Liu R L.2011. A stem taper model with nonlinear mixed effects for dahurian larch. Scientia Silvae Sinicae , 47 (4) : 101–106. [in Chinese] )
[] 李春明, 张会儒. 2010. 利用非线性混合模型模拟杉木林优势木平均高. 林业科学 , 46 (3) : 89–95.
( Li C M, Zhang H R.2010. Modeling dominant height for Chinese Fir plantation using a nonlinear mixed-effects modeling approach. Scientia Silvae Sinicae , 46 (3) : 89–95. [in Chinese] )
[] 李永慈, 唐守正. 2004. 用Mixed和Nlmixed过程建立混合生长模型. 林业科学研究 , 17 (3) : 279–283.
( Li Y C, Tang S Z.2004. Establishment of tree height growth model based on mixed and Nlmixed of SAS. Forest Research , 17 (3) : 279–283. [in Chinese] )
[] 祖笑锋, 倪成才, NighGorden, 等. 2015. 基于混合效应模型及EBLUP预测美国黄松林分优势木树高生长过程. 林业科学 , 51 (3) : 89–95.
( Zu X F, Ni C C, Nigh G, et al.2015. Based on mixed-effects model and empirical best linear unbiased predictor to predict growth profile of dominant height. Scientia Silvae Sinicae , 51 (3) : 89–95. [in Chinese] )
[] Calama R, Montero G.2004. Interregional nonlinear height-diameter model with random coefficients for stone pine in Spain. Can J For Res , 34 (1) : 150–163. DOI:10.1139/x03-199
[] Calegario N, Daniels R F, Maestri R, et al.2005. Modeling dominant height growth based on nonlinear mixed-effects model:a clonal eucalyptus plantation case study. Forest Ecology and Management , 204 (1) : 11–21. DOI:10.1016/j.foreco.2004.07.051
[] Ni C, Nigh G.2012. An analysis and comparison of predictors of random parameters demonstrated on planted loblolly pine diameter growth prediction. Forestry , 85 (2) : 271–280. DOI:10.1093/forestry/cps001
[] Davidian M, Giltinan D M.1993. Analysis of repeated measurement data using the nonlinear mixed effects model. Chemom. Intell Lab Syst , 20 (1) : 1–24. DOI:10.1016/0169-7439(93)80017-C
[] Davidian M, Giltinan D M. 1995. Nonlinear models for repeated measurement data. Chapman & Hall, New York, NY, 359. https://www.crcpress.com/nonlinear-models-for-repeated-measurement-data/davidian-giltinan/p/book/9780412983412
[] Fang Z, Bailey R L.2001. Nonlinear mixed effects modeling for slash pine dominant height growth following intensive silvicultural treatments. Forest Science , 47 (3) : 287–300.
[] Nigh G.2004. A comparison of fitting techniques for ponderosa pine height-age models in British Columbia. Ann For Sci , 61 (7) : 609–615. DOI:10.1051/forest:2004064
[] Hall D B, Bailey R L.2001. Modeling and prediction of forest growth variables based on multilevel nonlinear mixed models. For Sci , 47 (3) : 311–321.
[] Hall D B, Clutter M.2004. Multivariate multilevel nonlinear mixed effects models for timber yield predictions. Biometrics , 60 : 16–24. DOI:10.1111/biom.2004.60.issue-1
[] Huang S, Meng S X, Yang Y.2009a. Estimating a non-linear mixed volume-age model with and without taking into account serially-correlated errors:differences and implications. Modern Appl Sci , 3 (5) : 3–20.
[] Huang S, Wiens D P, Yang Y, et al.2009b. Assessing the impacts of species composition, top height and density on individual tree height prediction of quaking aspen in boreal mixed woods. For Ecol Manage , 258 (7) : 1235–1247. DOI:10.1016/j.foreco.2009.06.017
[] Jiang L, Li Y.2010. Application of nonlinear mixed-effects modeling approach in tree height prediction. J Computers , 5 (10) : 1575–1580.
[] Lappi J, Bailey R L.1988. A height prediction model with random stand and tree parameter:an alternative to traditional site methods. Forest Science , 34 (4) : 907–927.
[] Lindstrom M L, Bates D M.1990. Nonlinear mixed effects models for repeated measures data. Biometrics , 46 (3) : 673–687. DOI:10.2307/2532087
[] Meng S X, Huang S.2009. Improved calibration of nonlinear mixed-effects models demonstrated on a height growth function. For Sci , 55 (3) : 238–248.
[] Sharma M, Parton J.2007. Height-diameter equations for boreal tree species in Ontario using a mixed-effects modeling approach. For Ecol Manage , 249 (3) : 187–198. DOI:10.1016/j.foreco.2007.05.006
[] Sharma M, Parton J.2009. Modeling stand density effects on taper for jack pine and black spruce plantations using dimensional analysis. For Sci , 55 (3) : 268–282.
[] Temesgen H, Monleon V J, Hann D W.2008. Analysis and comparison of nonlinear tree height prediction strategies for Douglas-fir forests. Can J For Res , 38 (3) : 553–565. DOI:10.1139/X07-104
[] Wolfinger R D, Lin X.1997. Two taylor-series approximation methods for nonlinear mixed models. Comput Stat Data Anal , 25 (4) : 465–490. DOI:10.1016/S0167-9473(97)00012-1
[] Yang Y, Huang S.2013. On the statistical and biological behaviors of nonlinear mixed forest models. European Journal of Forest Research , 132 (5/6) : 727–736.