混合效应模型是指模型中部分或全部参数由固定效应和随机效应2部分组成。依据2种参数效应在模型中出现的形式不同，把混合效应模型分为线性混合效应模型(linear mixed effects model，简称为LMEMs)和非线性混合效应模型(nonlinear mixedeffects model，简称为NLMEMs)。2种参数效应都以线性形式出现的混合模型称为LMEMs，以非线性形式出现的混合模型称为NLMEMs。最早利用随机效应分析重复调查数据是在20世纪70年代(Harville，1977)，但没有正式给出混合模型的具体概念。直到1981年，Laird等(1981)在《Biometrics》杂志上详细地介绍了LMEMs。之后，混合模型受到越来越多的学者青睐，同时也提出了较多参数估计方法，如极大似然法(Laird et al.，1982)、Expectation-Maximization(EM)算法(Laird et al.，1982)、Newton-Raphson算法(Lindstrom et al.，1988)等。对于NLMEMs，它比LMEMs更加复杂，因此关于它的研究相对较晚，直到20世纪90年代初期，一些学者才逐渐提出该模型的参数估计方法(Sheiner et al.，1980;Lindstrom et al.，1990;Vonesh，1986;Vonesh et al.，1992)。到目前为止，计算NLMEMs参数的方法主要有广义最小二乘估计法(generalized least squares estimation)、一阶线性化算法(first-order linearization，简称为FO)、条件一阶线性化算法(first-order conditional expectationlinearization，简称为FOCE)、高斯-埃尔米特求积法(Gaussian-Hermite quadrature)、拉普拉斯近似法(Laplacian approximation)等。这些方法中，广义最小二乘估计法精度最低，高斯-埃尔米特求积法精度高但计算量大且运算速度慢，而FO和FOCE 2种算法计算简单且精度较高，故使用较广。一些国际主流统计软件，如SAS软件的NLINMIX宏、S-Plus和R软件的nlme模块等均实现了这2种算法。在林业上，最近几年国内一些学者利用FO和FOCE对NLMEMs进行了计算(李永慈等，2004;李春明等，2010;符利勇等，2011;2012)，但在选用参数估计方法上许多学者很难确定选用哪种方法最合适，到目前为止国内也没有关于这2种算法比较的报道。本研究基于FOCE算法，提出一种改进的随机效应参数计算方法，并对FO和FOCE以及改进的FOCE算法进行比较。文中所有计算均是在SAS9.0版本上实现的。

Lindstrom等(1990)定义NLMEMs为:

式中: y_ij和υ_ij分别为第i个研究对象第j次重复观测的因变量和自变量值;M为研究对象的个数;n_i为第i个研究对象重复观测次数;ε_ij为第i个研究对象第j次重复观测的误差项值;φ_ij是以非线性形式出现在函数f中的形式参数向量，写为:

式中: β为p维固定效应参数;u_i为q维随机效应参数，假定服从期望为零方差为ψ的正态分布;A_ij和B_ij分别为β和u_i的设计矩阵。假定ε_ij服从期望为零方差为σ²的正态分布，同时假定u_i之间、u_i与ε_ij之间相互独立。

FO和FOCE都是极大化似然函数求解模型参数。在NLMEMs中，因变量y的边缘密度函数由下式得到:

式中: u=(u₁，…，u_M)为随机效应参数向量;p(y |β，σ²，ψ)为y的边缘密度函数;p(y | u，β，σ²)为给定u时y的条件密度函数;p(u | ψ)为u的密度函数。由于β和u_i(i=1，…，M)是以非线性形式出现在模型(1)中，因此式(3)求解较为困难。许多学者提出了不同的计算方法，本研究中的FO是对式(3)在u=0处进行一阶泰勒展开，而FOCE算法是在u的条件期望处进行一阶泰勒展开。已知y的边缘密度函数，似然函数为:

FO和FOCE的参数估计值是由Newton-Raphson迭代方法优化式(4)(Lindstrom et al.，1990)得到的，详细介绍见符利勇等(2012)。

在FOCE中，随机效应u_i通过式(5)计算得到(符利勇等，2012):

式中: X_i和Z_i分别为β和u_i的设计矩阵;ω_i为伪向量，具体计算公式见文献(符利勇等，2012)。由式(5)得知，${\hat u_i}$与自身方差ψ、误差项方差σ² I_ni、设计矩阵Z_i 、因变量y_i以及自己本身有关。同时Z_i也与u_i有关，因此可以对式(5)进一步改进。假定FOCE中2个交替步骤运算第t步后计算收敛，对应ψ，σ²，Z_i和u_i的估计量分别为ψ^(t)，(σ^(t))²，Z_i^(t)及u_i^(t)，改进的FOCE算法详细计算过程如下。

第一步:计算随机效应u_i⁽⁰⁾的初始值。设ψ，σ²和$\hat \beta $的取值为ψ^(t)，(σ^(t))²和β^(t)，其中k=0，计算

其中:

第二步:计算

利用下式计算u_i^(k+1):

令k=k+1。

第三步:重复第二步，直到达到理想精度为止。迭代终止条件为:

δ为任意给定的正小数。

显然，改进的FOCE算法是在FOCE算法计算收敛后对随机效应参数u_i进一步迭代计算。模型中其他参数如固定效应参数和方差参数，2种算法完全相同。

本文选用中国林业科学研究院林业研究所设在江西省大冈山实验局种植的16株人工杉木(Cunninghamia lanceolata)作为试验对象，连续10次调查每株杉木树高，每次调查时杉木年龄分别为3，4，5，6，7，8，9，10，12和14年，样木树高的具体信息见表 1。

表 1 16株杉木各自树高统计量 Tab.1 Summary statistics of height for 16 Chinese fir

表 1 16株杉木各自树高统计量 Tab.1 Summary statistics of height for 16 Chinese fir 树号 Tree No. 最小值 Min./m 最大值 Max./m 平均值 Mean/m 方差 Variance 树号 Tree No. 最小值 Min./m 最大值 Max./m 平均值 Mean /m 方差 Variance 1 4.09 15.92 10.44 14.12 9 3.93 12.85 8.99 8.35 2 3.82 12.83 9.01 8.22 10 3.98 11.15 7.96 5.42 3 4.08 15.43 10.18 12.87 11 3.62 11.07 7.78 5.32 4 4.12 13.93 9.60 9.59 12 3.78 11.40 8.10 5.82 5 4.02 12.07 8.71 6.48 13 4.18 11.62 8.46 6.18 6 3.84 14.55 9.64 11.90 14 4.03 11.10 7.72 5.32 7 4.23 13.08 9.36 8.00 15 3.69 12.53 8.42 7.74 8 3.87 12.07 8.41 6.75 16 4.03 12.63 8.72 7.06

利用NLMEMs对杉木树高生长量进行分析需首先选定基础模型，基础模型的选择方法见文献(Hall et al.，2001;Sharma et al.，2007)。本文选取常用的Richards模型作为构建混合模型的基础模型，以模型中考虑2个混合效应参数为例，模型写为:

式中: H为杉木树高，m;A为杉木年龄;A₀杉木标准年龄，本文取A₀=20年;β₁，β₂和β₃为固定效应参数;u₁和u₂为随机效应参数。

选用平均偏差($\bar e$)以及结合平均偏差($\bar e$)和偏差的标准差(SD)的综合指标(total error of meanbias and standard deviation of the biases，简称TE)对3种算法预测精度进行比较，指标的计算公式为:

式中: $\bar e = \sum\limits_{i = 1}^M {\sum\limits_{j = 1}^{{n_i}} {{e_{ij}}/N = } } \sum\limits_{i = 1}^M {\sum\limits_{j = 1}^{{n_i}} {\left({{H_{ij}} - {{\hat H}_{ij}}} \right)/N;} } {\rm{SD = }}\sqrt {\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^{{n_i}} {{{\left({{e_{ij}} - \bar e} \right)}^2}/\left({N - 1} \right)} } } ;$H_ij为第i株杉木第j次观测的树高值;${{{\hat H}_{ij}}}$为第i株杉木第j次调查时树高的预测值。

本文通过拟合研究进一步比较3种算法各自的参数估计精度。根据模型(11)和实例分析中由FOCE算法计算的参数值作为真实值产生500组模拟数据，利用3种算法对每组模拟数据进行计算得到各自的参数估计值。结合参数的真实值，利用偏差、方差和均方误差对3种算法估计精度进行比较。3种算法都选用限制极大似然参数估计方法(唐守正，2002)。

根据建模数据，利用FO，FOCE和改进的FOCE算法分别对模型(11)进行计算，得到随机效应参数u₁和u₂的估计量见表 2，固定效应参数β₁，β₂和β₃以及方差参数的估计量见表 3。

表 2 3种算法各自的随机参数估计量^① Tab.2 The random parameters estimates for the three algorithms

表 2 3种算法各自的随机参数估计量① Tab.2 The random parameters estimates for the three algorithms 树号Tree No. ${\hat u}$Z1 ${\hat u}$Z2 ${\hat u}$E1 ${\hat u}$E2 ${\hat u}$I1 ${\hat u}$I2 1 4.049 8 -0.037 2 4.126 9 -0.031 8 4.108 4 -0.032 6 2 0.232 0 0.000 8 0.189 5 -0.000 5 0.190 8 0.000 1 3 3.373 3 -0.030 6 3.409 1 -0.027 2 3.398 2 -0.028 9 4 1.412 3 -0.005 6 1.361 8 -0.006 4 1.376 1 -0.006 0 5 -0.936 5 0.017 2 -0.857 8 0.014 3 -0.877 2 0.014 9 6 2.490 2 -0.029 3 2.511 9 -0.026 9 2.497 1 -0.027 1 7 0.366 3 0.010 5 0.390 2 0.007 3 0.387 3 0.007 9 8 -0.935 8 0.006 3 -0.972 6 0.005 6 -0.965 1 0.005 8 9 0.288 4 -0.001 4 0.236 2 -0.002 5 0.241 7 -0.002 2 10 -2.034 9 0.015 6 -2.003 9 0.015 5 -2.018 3 0.015 5 11 -2.263 4 0.014 5 -2.255 8 0.015 4 -2.258 6 0.015 0 12 -1.756 4 0.014 3 -1.730 3 0.013 7 -1.742 4 0.013 9 13 -1.281 1 0.060 6 -1.217 8 0.013 9 -1.219 5 0.017 3 14 -2.207 4 0.010 5 -2.256 5 0.012 2 -2.248 0 0.011 9 15 -0.339 9 -0.007 5 -0.444 3 -0.007 1 -0.440 2 -0.007 3 16 -0.457 1 0.005 8 -0.486 8 0.004 6 -0.472 0 0.004 9 ① ${\hat u}$Z1和${\hat u}$Z2是由FO算法计算出的u1和u2的估计值;${\hat u}$E1和${\hat u}$E2是由FOCE计算出的u1和u2的估计值;${\hat u}$I1和${\hat u}$I2是由改进的FOCE算法计算出的u1和u2的估计值。${\hat u}$Z1 and ${\hat u}$Z2 calculated by FO algorithm are estimators of u1 and u2,respectively.${\hat u}$E1 and ${\hat u}$E2 calculated by FOCE algorithm are estimators of u1 and u2,respectively.${\hat u}$I1 and ${\hat u}$I2 calculated by improved FOCE algorithm are estimators of u1 and u2,respectively.

表 3 FO及FOCE算法各自参数及指标^① Tab.3 The parameters and indexes of FO and FOCE algorithm

表 3 FO及FOCE算法各自参数及指标① Tab.3 The parameters and indexes of FO and FOCE algorithm 算法Algorithms β1 β2 β3 su1 su2 σu1u2 s AIC AICC BIC FO 14.093 4 0.173 3 1.480 3 3.625 4 0.000 3 -0.031 5 0.048 5 86.2 86.9 91.6 FOCE 14.054 9 0.179 9 1.513 2 3.662 3 0.000 3 -0.029 0 0.049 2 85.6 86.3 91.0 ① β1,β2和β3为固定效应参数;su1和su2分别为随机效应u1和u2的标准差;σu1u2为随机效应u1和u2的协方差;s为随机误差的标准差;AIC,AICC和BIC分别为评价指标。β1,β2 and β3 are fixed effects parameters.su1 and su2 are standard deviation of random effects u1 and u2,respectively.σu1u2 is covariance of the random effects u1 and u2.AIC,AICC and BIC are evaluation indexes.

从表 2得知，3种算法对应的随机效应差异不大，其中FOCE算法和改进的FOCE算法对应的随机效应差异最小，而FO算法与FOCE算法对应的随机效应差异最大。从表 3中可以看出FOCE算法的3个指标AIC=85.6，AICC=86.3和BIC=91比FO算法中对应的3个指标AIC=86.2，AICC=86.9和BIC=91.6略小些，根据“越小越好”原则判断FOCE算法比FO算法要好。在利用NLMEMs分析林分特征因子生长量时，由于FO算法只考虑随机效应以非线性形式出现在模型中，而FOCE算法同时考虑了固定效应和随机效应，因此理论上讲选用FOCE算法计算参数估计值比选用FO算法效果要好。但在实际应用中，这2种算法拟合效果还与具体数据和模型有关(符利勇等，2012)。

利用3种算法对实际数据进行拟合，得到各自的平均偏差和指标TE。对于FO算法2个指标值分别为: ${\bar e}$=0.081 3，TE=0.24;FOCE算法分别为:${\bar e}$=0.002 8，TE=0.203 7;改进的FOCE算法分别为: ${\bar e}$=0.002 1，TE=0.196 2。很明显改进的FOCE算法拟合效果最好，FOCE次之，而FO最差。从而说明本文给出的改进FOCE算法所计算出的随机效应参数更能真实反映研究对象间的随机差异。图 1为3种算法对应的2株样本优势高拟合曲线。

	图 1 3种算法对应的2株样木树高拟合曲线 Fig. 1 Two sample trees height fitted curves for three algorithms 试验为实际调查数据。The full lines indicate the actual survey data.

表 4和表 5分别为3种算法对应的方差参数和固定效应参数拟合结果。从表 4中可知，FO算法对应的偏差比FOCE算法要大，其中，参数s_u1最为显著，2种算法相差0.001 7，而对于参数s_u22种算法相差最小，对应的偏差为0，因此说明FOCE算法计算模型(11)方差参数时精确度要比FO算法高。对于均方误差，FO算法和FOCE算法之间非常接近，但并不由于FOCE算法精度高而对应的均方误差比FO算法小，如参数s_u1和参数δ_u1u2对应的均方误差FOCE算法比FO算法要大。对于FOCE算法和改进的FOCE算法，由于改进的FOCE算法只对随机效应参数的计算进行改进，而模型中其他参数并没有变，与此同时随机效应参数不影响拟合计算，故2种算法计算出的结果是完全一样的。从表 4中还可看出误差项方差s对应的偏差比另外几个随机效应方差参数(s_u2除外)要小，这可能是由于误差项方差s与总观测次数有关，而随机效应方差参数与观测对象个数有关。表 5中3种算法对模型(11)中固定效应参数拟合结果与表 4中方差参数拟合结果完全一致。

表 4 模型(11)中方差参数拟合结果 Tab.4 Simulation results for the variance-covariance components in model(11)

表 4 模型(11)中方差参数拟合结果 Tab.4 Simulation results for the variance-covariance components in model(11) 算法 Algorithms Su1 Su2 平均值Mean 偏差Bias 均方误差MSE 平均值Mean 偏差Bias 均方误差MSE FO 3.659 4 0.002 9 0.653 2 0.000 3 0 0.0062 FOCE 3.661 1 0.0012 0.658 1 0.000 3 0 0.005 4 改进的FOCE Improved FOCE 3.661 1 0.0012 0.658 1 0.000 3 0 0.005 4 算法Algorithms δu1u2 s 平均值Mean 偏差Bias 均方误差MSE 平均值Mean 偏差Bias 均方误差MSE FO -0.031 2 0.0022 0.0161 0.048 9 0.000 3 0.015 3 FOCE -0.029 9 0.000 9 0.016 7 0.049 0 0.0002 0.014 9 改进FOCE Improved FOCE - 0.029 9 0.000 9 0.016 7 0.049 0 0.0002 0.014 9

表 5 模型(11)中固定效应参数拟合结果 Tab.5 Simulation results for the fixed effects parameters in model(11)

表 5 模型(11)中固定效应参数拟合结果 Tab.5 Simulation results for the fixed effects parameters in model(11) 算法 Algorithms β1 β2 β3 平均值Mean 偏差Bias 均方误差MSE 平均值Mean 偏差Bias 均方误差MSE 平均值Mean 偏差Bias 均方误差MSE FO 14.051 3 0.0036 1.472 3 0.179 1 0.0008 0.0648 1.511 7 0.001 5 0.265 1 FOCE 14.054 79 0.00011 1.475 9 0.179 5 0.000 4 0.062 5 1.512 3 0.000 9 0.268 4 改进的FOCE 14.054 79 0.00011 1.475 9 0.179 5 0.000 4 0.062 5 1.512 3 0.000 9 0.268 4

由于FOCE算法和改进的FOCE算法参数估计结果完全一致，因此选用FOCE算法与FO算法进行比较分析，图 2为2种算法的散点分布图。从图 2中可以看出，每个方差参数的散点分布图几乎都在曲线y=x以下，这表明FOCE算法算出来的参数值比FO算法要大。参数s_u1具有较大的均方误差(表 4)，在图 2中可明显反映出来。鉴于文章篇幅所限，固定效应参数的散点图不再列出，但结果与方差参数完全一致。

	图 2 FO和FOCE对模型(11)方差参数估计的散点分布 Fig. 2 Scatter plots of variance-covariance components etimates for the two algorithms of FO and FOCE in the Eq. 11

线性化算法是把非线性混合模型通过一阶泰勒展开转化为线性混合模型进行求解，它避免了求解y_i边缘分布难的问题。FO算法中只考虑随机效应u_i在模型中以非线性形式出现，并且假定u_i=0时对参数进行估计;而正常情况下固定效应β是以非线性形式出现在模型中，并且随机效应不为零。许多国外学者喜欢选用FO算法来求解NLMEMs的未知参数(Fang et al.，2001)，但是Davidian等(2003)证明了E{f(β，u_i，υ_i)} ≠f(β，υ_i)，因此利用FO算法对因变量进行估计时误差将会很大，本文实例分析也证明了这点。FOCE算法是基于FO算法的基础上考虑固定效应以非线性形式出现在模型中，并且取随机效应ui的最小二乘估计量作为初始值进行迭代计算。本文基于FOCE提出一种改进的随机效应参数计算方法，利用树高生长模型实例和拟合研究对3种算法进行比较和分析。通过树高生长模型实例可知，改进的FOCE算法对应的指标e-和TE最小，拟合效果最好。FOCE算法和改进的FOCE算法唯一的不同点是二者具有不同的随机效应参数值，其中改进的FOCE算法计算出的随机效应参数更能真实反映研究对象间的随机差异，对于其他参数估计值，2种算法完全相同。通过拟合研究可得，FOCE算法计算的参数值比FO算法精度更高，同时由FOCE算法计算的参数值几乎都比FO算法计算出的要大。除随机效应参数外，对于其他参数，FOCE算法与改进的FOCE算法有相同的计算精度。总而言之，在模型应用时，利用改进的FOCE算法计算随机效应参数能进一步提高模型拟合效果。

除线性化算法来求解NLMEM外，还有一些其他算法，如高斯积分法和拉普拉斯法同样可以对随机效应参数的计算进行改进，这需进一步研究。