Based on mixed-effects model and empirical best linear unbiased predictor predicting growth profile of height for Chinese fir

研究区位于福建省将乐县国有林场，26°26′~27°04′N，117°05′~117°40′E。将乐县位于福建省西北部，地处武夷山脉东南部，以中、低山为主，最高峰海拔为1 640.2 m。属亚热带季风气候，具有海洋性和大陆性气候特点，年平均气温为19.8 ℃，年平均降水量为2 027.0 mm。境内气温较高，夏季时间长，冬季较温暖，霜冻较少，生长期长。土壤以红壤为主，并分布有黄红壤，土层深厚，土质较好的一般为沙壤土或轻壤土，水分充足，土壤肥沃。植被以亚热带植物区系为主，植被种类非常丰富。其中人工次生林的乔木主要有杉木、马尾松Pinus massoniana等；灌木主要有粗叶榕Ficus hirta，黄毛楤木Aralia decaisneana，檵木Loropetalum chinensis等；草本主要有乌毛蕨Blechnum orientale，乌蕨Stenoloma chusanum，铁线蕨Adiantum capillus-veneris等。

2010-2012年，根据林分不同年龄和密度，以典型抽样原则设置了15块20 m × 30 m标准地。对标准地内的林木进行每木检尺，选取标准木2株·标准地^-1，共30株，伐倒并进行树干解析，解析木编号为F01~F30。根据解析木内业表内插法获得各年龄对应的树高值。随机选择20株进行拟合，10株进行预测（表 1）。

选择Richards方程、Weibull方程、Korf方程以及Logistic方程的2种变化形式等5个方程作为基础模型，对树高生长进行模拟（表 2）。运用非线性最小二乘法对5个树高生长方程进行拟合，选出拟合效果最好的模型作为基础模型，构建非线性混合效应模型^[17]。

非线性混合效应模型主要特点是参数分为固定效应和随机效应，固定效应可以反映研究对象的总体变化规律，随机效应随个体的不同而变化，反映总体中不同个体的变化规律，从而获得较好的拟合效果。笔者主要研究利用EBLUP预测树高生长过程，并探究提高预测精度的方法。因此将样木作为随机效应，建立混合效应模型后，通过选取生长过程与总体平均生长过程差异较大的样木进行预测分析。非线性混合模型的一般表达式为：

式（1）中：β为p×1维的固定效应向量；p为固定参数个数，b_i为q×1维的随机效应向量；q为随机参数个数；A_i和B_i为具有相应维度的设计矩阵，其元素通常为0，1或与固定效应和随机效应相关的协方差值；e_i为随机误差项向量；f为非线性函数；X_i为年龄；y_i为树高。其中随机向量b_i：N（0, D），e_i：N（1, R_i），并且相互独立；e_i和b_i协方差矩阵分别为R_i和D，且为n_i维和q维的对称矩阵。

构建非线性混合效应模型还要确定以下3种结构：① 混合效应参数。依据PINHEIRO等^[18]的研究，将基础模型中参数进行组合作为混合参数进行模拟，选择收敛的模型，通过比较其赤池信息量准则（AIC），贝叶斯信息准则（BIC）以及对数似然函数（Loglik）模型评价指标。AIC，BIC值越小，Loglik值越大，模型拟合效果越好。② 随机效应内的方差协方差矩阵。为了确定随机效应内的方差协方差矩阵，需要解决异方差和自相关性结构。以指数函数（Exp），幂函数（Power）和常数加幂函数（ConstPower）3种结构形式，消除数据间的异方差。表达式为：

式（2）~式（5）中：σ²为误差方差值，由模型残差方差值给出；G_i为n_i×n_i维对角矩阵来解释方差异质性，对角元素为相应误差项的标准差；由于本研究数据间不存在自相关性，因此，I_ni为n_i×n_i维的单位矩阵。β，β₁，β₂为参数，t为年龄^[19-22]。从自变量为年龄的3种函数中，根据AIC，BIC，Loglik值以及似然比检验，确定残差方差模型。③ 随机效应间的协方差矩阵。随机效应间的协方差矩阵反映了随机效应之间的变化性。本研究包含3个随机参数，因此，矩阵为3 × 3维的方差协方差矩阵，结构例如下：

式（6）中：σ_u²，σ_v²，σ_w²分别为随机参数u，v，w的方差，σ_uv=σ_vu为随机参数u和v的协方差，σ_uw=σ_wu为随机参数u和w的协方差，σ_vw=σ_wv为随机参数w和v的协方差。

推导方法基于多元正态分布理论，以b_i=0为基点，用一阶泰勒公式将非线性混合效应模型式（1）近似地表达为：

式（7）中：矩阵Z_i为q × n_i维矩阵，其元素可以通过对混合效应模型式（1）分别求关于随机效应b_i的偏导数，随后令所有b_i取0即可求出。前文提到b_i：N（0, D），e_i：N（1, R_i），且相互独立，因此对式（7）求y_i的方差矩阵为：

b_i和y_i间协方差矩阵为：

b_i和y_i的联合分布为：

根据多元正态分布理论，当 ${y_i} = {{\tilde y}_i}$ （ ${{\tilde y}_i}$ 为观测值）时，b_i的条件分布数学期望为：

对于随机变量取值最优预测值是其数学期望值，所以随机变量b_i在y_i取观测值时可以用b_i的条件分布数学期望表示：

将式中β，D，R，Z_i及f（A_i β, X_i）用相应的估计值取代后则：

式（8）称为经验线性无偏最优预测法（EBLUP）^[16]。

设定 ${y_i} = \left( {\begin{array}{*{20}{l}} {{y_{i,1}}} \\ {{y_{i,2}}} \end{array}} \right)$ ，y_{i, 1}为观测值，y_{i, 2}为预测值。则y_{i, 1}和y_{i, 2}的联合分布为：

其中：W₂₁=Cov（y_{i, 2}, y_{i, 1}）=Z_{i, 2}DZ_{i, 1}^T+Cov（e_{i, 1}^T, e_{i, 2}），W₁₁=Var（y_i）=Z_{i, 1}DZ^T_{i, 1}+R_{i, 1}，则y_{i, 2}的条件分布数学期望：

将各项参数用相应估计值带入的：

简化为：

结合公式（8）和公式（9）即可利用混合效应模型进行EBLUP预测。详细推导过程参考文献[23]。

根据EBLUP预测原理，式（9）中Z_{i, k}可通过求关于随机效应b_i的偏导数，求出偏导数后令所有b_i取0即可。本模型有3个随机效应，故k取1，2，3。Z_{i, k}的估计值为：

列向量 $f\left( {{A_i}\hat \beta ,{X_i}} \right)$ 表达式为：

式（11）~式（14）中： ${{\hat \beta }_{i,1}},{{\hat \beta }_{i,2}},{{\hat \beta }_{i,3}}$ 为β_{i, 1}，β_{i, 2}，β_{i, 3}的估计值，通过R语言软件得出。随机效应参数的预测值 ${{\hat b}_{i,1}},{{\hat b}_{i,2}},{{\hat b}_{i,3}}$ 根据式（8），利用SAS 9.3的IML过程计算获得。将各参数估计值代入式（9）进行预测。

利用R语言软件基于表 1中拟合数据，对表 2中6个树高生长方程进行拟合，拟合结果见表 3。根据决定系数R²越接近于1，均方根差R_mse以及平均绝对残差|E|越接近于0，其拟合效果越好的原则，选取方程1，方程2，方程3作为基础模型。

基于基础模型，将单木考虑作为随机效应，利用R语言软件中的nlme函数对模型进行估计，改变随机参数个数比较拟合效果。基于不同随机参数组合的混合模型的拟合效果见表 4。从表中拟合情况可以看出模型2.5的AIC，BIC最小且Loglik值最大，即当模型2以β₁，β₂和β₃同时作为混合效应参数时，模型拟合效果最佳。因此，选择模型2.5（Weibull方程）进一步构建混合效应模型。

加权残差方差模型的评价指标见表 5。考虑残差加权的3种模型中幂函数的AIC和BIC最小且Loglik值最大，故选它作为模型的残差方差模型。

利用R语言软件拟合模型的各参数列于表 6，基于非线性混合效应模型的杉木树高生长方程表达式为：

由于EBLUP算法的原因，本研究利用式（9）预测，具体原因参考NI等^[23]。从验证数据中选取F04解析木作为预测数据，选3对观测值用于预测b_i，结果统计于表 7。表中 ${{\tilde y}_i}$ 为样木在X_i时树高观测值， $f\left( {{A_i}\hat \beta ,{X_i}} \right)$ 通过式（13）计算得出， ${{\tilde Z}_i}$ 通过式（10）~式（12）计算得出， ${{\tilde b}_i}$ 是通过SAS 9.3的IML过程获得。利用表 7的 ${{\tilde b}_i}$ 预测F18样木19~22 a间的树高值，并与观测值比较列于表 8。总体估计值通过式（13）计算得到，预测值通过式（9）计算获得，校正值为预测值减总体估计值的差。可见根据3对观测值获得的对样木进行EBLUP预测树高，能较准确地反映样木生长过程。在验证数据中选择样木生长过程与总体平均生长过程有差异的F04，F09，F11，F18解析木，进行EBLUP预测，将预测生长曲线、解析木树高曲线和总体平均生长过程进行比较（图 1）。通过图 1可看出：当单木生长过程和总体估计值差距较大时，基于Weibull混合效应模型的EBLUP预测能充分地预测单木生长过程。本研究依据EBLUP预测的算法，利用样木的3组观测值，基于混合效应模型使每个样木获得对应的随机效应参数，进而充分地预测样木的生长过程。

Figure 1.  Comparison of the height growth

在进行单木树高生长预测时发现，观测值的选择方式对预测精度高低有较大的影响，因此本研究通过设置不同的观测值间隔、观测值选取次数以及预测时长，比较预测结果的精度高低。评价预测精度的指标为均方误差（E_MS）：

以解析木F04，F18，F15和F11为例（由于这几株样木生长过程与总体平均生长过程差距较大，通过不同处理方式能较明显反应预测精度的提高），通过不同的选择观测值方式来研究对E_MS的影响，其中观测间隔指每隔多少年对树高进行一次观测。本研究分别选取间隔1，3，5 a；观测次数是指观测多少个树高值，由于本研究基于具有3个随机效应参数的混合模型，预测随机效应参数至少需要3对观测值，故分别选取3，6，9次观测次数；相应地，因为每3对观测值会求得1组随机效应参数，故观测次数为6，9次时，需要分别取2，3组数据，分别按照1/2，1/3将解析木数据分段，年龄向上取整，每年龄段从最大值开始向下取观测值。以样木F15（35年生）为例，观测间隔5 a观测3次时，使用的是35，30，25年生的树高观测值；观测间隔为3 a观测6次时，分别使用35，32，29年生以及18，15，12年生这2组共6个树高观测值；观测间隔为1 a观测9次时，分别使用35，34，33年生，24，23，22年生以及12，11，10年生这3组共9个树高观测值。由于F04，F18年龄较小，故没有进行间隔5 a 9次的观测。以树高观测值为基础，根据式（8）预测b_i值，再对树高逐年预测最后求总均方误差。结果列于表 9。

观测次数相同时，随着观测间隔的增大，均方误差降低。观测次数均为3，观测间隔取1，3，5 a，各解析木的均方误差随着观测间隔的增大而明显降低。观测间隔相同时，随着观测次数的增加，均方误差降低。观测间隔均为3 a时，随着观测次数从3次增加到9次，各组的均方误差均相应降低，观测间隔为1或5 a按同样方式比较，也可以得出相同的结论。以F15为例，观测间隔均为1 a时，随着观测次数从3次增加到9次，均方误差从12.396 0降低到2.631 7。

从图 2中可以看出：“1 a 3次”在其观测值28，27，26年生附近均方误差较低，远离观测值时，均方误差开始变大。同样，“3 a 3次”在远离20年生时均方误差开始变大，而“5 a 3次”的观测值为28，23，18年生，故在图 2中均方误差较低。所以在EBLUP预测时，随着预测远离用来预测b_i的观测值，均方误差逐渐变大。

Figure 2.  The value of F18 yearly forecast E_MS

本研究选用5个生长方程，选取3个拟合效果较好的方程作为基础模型，进一步改变不同的随机效应组合，最终以式（2）为基础的Weibull方程中，3个参数都作为混合模型的模型精度最高，拟合效果最好。通过EBLUP预测表明：本研究建立的非线性混合效应模型在考虑样木随机效应时，预测效果明显优于固定模型，能更好地反映树高生长过程。将随机效应考虑到样木水平，就相当于为每株样木建立树高生长模型，可以充分预测样木的生长过程，能够为将乐林场的杉木林经营决策提供一定的依据。进行预测时，观测次数相同，增加观测间隔能够降低预测误差，提高预测精度。由于增加观测值间年龄间隔，可以更充分提供树高生长过程信息，因此可以显著降低预测误差。观测次数不同，预测精度变化较大。观测间隔相同时，增加观测次数，预测精度会提高。因此在预测杉木树高生长过程时，要根据其年龄选择适当的观测间隔和观测次数。EBLUP可以充分利用观测值所涵盖的生长过程信息降低预测误差，随着预测远离用来预测b_i的观测值，均方误差逐渐变大。

通过本次研究得到的数据不能明确观测间隔和观测次数对预测准确度影响，有关研究还需要进一步进行。