需要JoVE订阅才能观看此内容。 请登录或开始免费试用

方法文章

基于线性混合效应方法的单木胸高断面积生长量模型构建

3.1K 次观看

⸱

DOI：

10.3791/60827

⸱

2020年7月3日

 ,  ,  , 
1Research Center of Forest Management Engineering of National Forestry and Grassland Administration, Beijing Forestry University, 2Research Institute of Forestry, Chinese Academy of Forestry

本文内容

摘要

混合效应模型是分析林业中具有层次化随机结构数据的灵活且实用的工具，也可用于显著提升森林生长模型的性能。本文提供了一项综合线性混合效应模型相关信息的实验方案。

摘要

本研究基于中国西北地区新疆省779个样地内21898株粗枝云杉（Picea asperata）的数据集，构建了单木5年期胸径断面积增量模型。为避免同一抽样单元内观测值之间的高度相关性，采用具有随机样地效应的线性混合效应方法建模，以反映随机变异。模型将树木大小、竞争状况和立地条件等个体树和林分水平的多种变量作为固定效应，用于解释残差变异。此外，通过引入方差函数和自相关结构来描述异方差性和自相关性。最优线性混合效应模型的确定依据多个拟合统计指标：赤池信息准则（AIC）、贝叶斯信息准则（BIC）、对数似然值以及似然比检验。结果表明，影响单木胸径断面积增量的显著变量包括胸径的倒数变换、大于目标树的树木断面积、每公顷株数以及海拔。此外，指数函数能最有效地模拟方差结构中的误差，而一阶自回归结构（AR(1)）可显著校正自相关性。与普通最小二乘回归模型相比，线性混合效应模型的拟合性能显著提升。

引言

与同龄单一树种经营相比，近年来，具有多重目标的异龄混交林经营受到了越来越多的关注1,2,3。为了制定稳健的森林经营策略，尤其是针对复杂的异龄混交林，有必要预测不同经营方案的效果4。森林生长与收获模型已被广泛用于预测在不同经营措施下树木或林分的生长发育及采伐结果5,6,7。森林生长与收获模型可分为单木模型、径级模型和全林分生长模型6,7,8。然而，径级模型和全林分模型并不适用于异龄混交林，因为后者需要更详细的描述以支持森林经营决策过程。因此，在过去几十年中，单木生长与收获模型因其能够对具有多种物种组成、结构和经营策略的林分进行预测而受到越来越多的关注9,10,11。

普通最小二乘法（OLS）回归是构建单木生长模型最常用的方法12,13,14,15。单木生长模型的数据集通常在固定时间间隔内对同一采样单元（即样地或单株树木）进行重复观测，具有层次化的随机结构，导致观测值之间缺乏独立性，并表现出较高的空间和时间相关性10,16。这种层次化随机结构违背了OLS回归的基本假设，即残差独立、数据正态分布且方差齐性。因此，对这类数据使用OLS回归不可避免地会导致参数估计标准误的估计出现偏倚13,14。

混合效应模型为分析具有复杂结构的数据（如重复测量数据、纵向数据和多层次数据）提供了强有力的工具。混合效应模型包含固定效应部分和随机效应部分，其中固定效应适用于整个总体，而随机效应则针对每个抽样层次特有。此外，混合效应模型通过定义非对角线的方差-协方差结构矩阵，能够考虑空间和时间上的异方差性及自相关性17,18,19。因此，混合效应模型已在林业研究中得到广泛应用，例如用于胸径-树高模型20,21、树冠模型22,23、自然稀疏模型24,25以及生长模型26,27。

本研究的主要目标是采用线性混合效应方法建立单木胸高断面积生长量模型。我们希望混合效应方法能够得到广泛应用。

访问受限。请登录或开始试用以查看此内容。

方案

1. 数据准备

  1. 准备建模数据，包括单株树木信息（树种和胸径，胸径测量高度为1.3米）以及样地信息（坡度、坡向和海拔）。本研究的数据来源于中国西北部新疆维吾尔自治区的第八次（2009年）和第九次（2014年）全国森林资源清查，共包含779个样地的21,898个观测值。这些样地为正方形，面积为1亩（中国面积单位，相当于0.067公顷），在4 km × 8 km的网格上系统布设。
    注意：用于建模（胸高断面积）增量分析的数据至少需要一个生长周期（即两次观测）。
  2. 将数据随机分为两个数据集，其中80%的样地数据用于模型拟合（模型构建数据集），包含623个样地的17,145个观测值；20%用于模型验证（模型验证数据集），包含156个样地的4,753个观测值。关键变量的描述性统计信息见表1。
    注意：建模过程中的此步骤可以省略，全部数据可用于模型构建。
变量拟合数据验证数据
最小值最大值平均值标准差最小值最大值平均值标准差
DBH1 (cm)5124.819.913.25101.519.513.4
QMD (cm)6.782.322.58.59.273.321.89.2
ID (cm)0.114.41.110.116.911.1
BAL (m3)05.21.70.905.41.71
NT (trees/ha)14.936421072673.714.934181205829.3
BA (m2/ha)0.177.534.213.90.180.634.515.3
EL (m)233022189340.3144133802256308.3

表1. 拟合与验证数据的统计摘要。DBH1：初始胸径（距地面1.3 m处的胸径），DBH2：生长5年后的胸径，QMD：平方平均直径，ID：5年期直径生长量（DBH2 – DBH1），BAL：大于目标树木的邻体树木的断面积（目标树木：用于计算竞争指数的树木），NT：每公顷树木株数，BA：每公顷断面积，EL：海拔，S.D.：标准差。

2. 基础模型构建

  1. 参考文献以确定影响单株树木胸径生长量的变量。
  2. 根据数据选择并计算变量。通常，单株树木的胸高断面积生长量受三类变量影响：树木大小、竞争程度和立地条件27,28,29,30.
    1. 考虑树木尺寸效应，例如胸径（DBH）1，胸径的平方1 (DBH 平方方程；用于森林资源计算与数据分析的数学符号。)，DBH 的逆变换1 (1/DBH1），以及DBH的常用对数1 (logDBH1）或其组合。
    2. 考虑竞争效应时，应同时采用单侧和双侧竞争指数，以更全面地量化树木所经历的竞争程度及其在林分中的社会地位。单侧竞争包括林木竞争指数（BAL）和相对密度指数（RD=DBH1/QMD）；双向竞争包括 NT 和 BA。
      注：若有数据可用，应考虑距离依赖性竞争指数。
    3. 考虑位点效应，如坡向（ASP）、坡度（SL）和海拔（EL）。应使用Stage的转换方法纳入SL和ASP31.
  3. 选择对数( DBH符号的平方，公式表示，数学符号。 - DBH 平方方程；用于林业计算和数据分析的数学符号。 +1) ( DBH 符号的平方，公式表示，数学符号。 表示胸径的平方2作为因变量。
  4. 使用逐步回归法构建基本模型。确保模型具有生物学合理性，并且自变量之间存在显著差异。利用方差膨胀因子（VIF）检验多重共线性。
  5. 保持自变量不变 p < 0.05 和 VIF < 5 在基础模型中。
  6. 输出基础模型的结果和残差图。此处建立的基础模型可作为进一步构建混合效应模型的依据。

3. 使用 R 软件中的“nlme”包构建线性混合效应模型

  1. 读取模型开发数据集并加载“nlme”包。
    >model.development.dataset=read.csv("E:/DATA/JoVE/modelingdata.csv", 
    header=TRUE)
    >library(nlme)
  2. 选择样地作为随机效应，以构建混合效应模型。
  3. 使用最大似然法（ML）拟合所有可能的随机效应组合，并输出结果。
    >Model<-lme(Y~1/DBH1+BAL+NT+EL,data=model.development.dataset, 
    method="ML", random =~1|PLOT)
    >summary(Model)
    1. 将 random =~1 设为随机参数的截距项。修改随机效应语句，直至所有组合均被拟合。例如，若将 1/DBH1 和 BAL 设为随机参数，则代码如下：random =~1/DBH1+BAL-1。此外，在拟合过程中，由于模型未能收敛，代码可能会报错。
  4. 根据赤池信息准则（AIC）、贝叶斯信息准则（BIC）、对数似然值（Loglik）以及似然比检验（LRT）选择最优模型。
    >anova(Model.1, Model.6)
    >anova(Model.6, Model.23)
    >anova(Model.23, Model.30)
  5. 确定 Ri 的结构，处理 Ri32 的异方差性和自相关性。Ri 的表达式如下：
    静态平衡公式 \(R_i = \sigma^2 G_t^{0.5} \Gamma_i G_i^{0.5}\)，数学方程。   (1)
    其中 σ2 是一个未知的缩放因子，等于模型残差方差，Gi 是描述异方差性的对角矩阵，Γi 是描述自相关性的矩阵。
    1. 通过残差图判断残差是否存在异方差性。若存在异方差性（残差呈现明显模式或趋势），则引入三种常用的方差函数——常数加幂函数、幂函数和指数函数——以建模误差的方差结构。
      >Model.30.1<-lme(Y~1/DBH1+BAL+NT+EL,data=model.development.dataset, method="ML",random=~1/DBH1+BAL+NT|PLOT,
      weights=varConstPower(form=~ fitted(.)))
      >summary(Model.30.1)
      >Model.30.2<-lme(Y~1/DBH1+BAL+NT+EL,data=model.development.dataset, method="ML",random=~1/DBH1+BAL+NT|PLOT,
      weights=varPower(form=~ fitted(.)))
      >summary(Model.30.2)
      >Model.30.3<-lme(Y~1/DBH1+BAL+NT+EL,data=model.development.dataset, method="ML",random=~1/DBH1+BAL+NT|PLOT,
      weights=varExp(form=~ fitted(.)))
      >summary(Model.30.3)
    2. 根据 AIC、BIC、Loglik 和 LRT 确定模型的最佳方差函数。
      >anova(Model.30, Model.30.1)
      >anova(Model.30, Model.30.2)
      >anova(Model.30, Model.30.3)
    3. 引入三种常用的自相关结构——复合对称结构（CS）、一阶自回归结构 [AR(1)] 以及一阶自回归与移动平均结构的组合 [ARMA(1,1)]——以处理自相关性。
      >Model.30.3.1<-lme(Y~1/DBH1+BAL+NT+EL,data=model.development.dataset, method="ML",
      random=~1/DBH1+BAL+NT|PLOT, weights=varExp(form=~fitted(.)), corr= corCompSymm())
      >summary(Model.30.3.1)
      >Model.30.3.2<-lme(Y~1/DBH1+BAL+NT+EL,data=model.development.dataset,  method="ML",
      random=~1/DBH1+BAL+NT|PLOT,weights=varExp(form=~ fitted(.)), corr=corAR1())
      >summary(Model.30.3.2)
      >Model.30.3.3<-lme(Y~1/DBH1+BAL+NT+EL,data=model.development.dataset, method="ML",
      random=~1/DBH1+BAL+NT|PLOT,weights=varExp(form=~ fitted(.)), corr=corARMA(q=1,p=1))
      >summary(Model.30.3.3)
    4. 根据 AIC、BIC、Loglik 和 LRT 确定最优的自相关结构。
      >anova(Model.30.3, Model.30.3.2)
      注：若不存在异方差性和自相关性，则无法定义 Gi 和 Γi。
    5. 使用限制性最大似然法（REML）输出混合效应模型的最终结果。
      >Mixed.model<-lme(Y~1/DBH1+BAL+NT+EL,data=model.development.dataset, method="REML",random=~1/DBH1+BAL+NT|PLOT,
      weights=varExp(form=~ fitted(.)), corr=corAR1())
      >summary(Mixed.model)

4. 偏倚校正

  1. 将最终模型在对数尺度上预测的胸径增量值转换回原始尺度。然而，将对数转换模型的预测值进行线性反向转换会产生相应的对数转换偏差。为校正该对数偏差，推导出一个校正因子并将其整合到预测方程中，以估计特定树木的实际预测胸径增量 [公式 (2)]：
    静态平衡；公式：BA = exp(BA1) + (σ²plot + σ²) / 2 - 1；数学概念。   (2)
    其中 静态平衡；带角度的向量表示法，结构分析示意图。 为模型预测的胸径增量对数值，而 带重音符号的BAI角度表示法，几何图解中的数学符号表示。 为经对数转换偏差校正后，预测的5年胸径增量反向转换值。用于统计分析图的σ²plot符号。 为样地随机效应的方差，σ2 为残差方差。
  2. 将胸径增量（带重音符号的BAI角度表示法，几何图解中的数学符号表示。）转换为直径增量。

5. 模型预测与评估

  1. 准备第1.2节中生成的模型验证数据集以进行预测。
  2. 使用线性混合效应模型预测单株树木的胸高断面积增量。随机成分通过以下最佳线性无偏预测器计算：
    涉及线性回归估计的统计分析数学公式。   (3)
    其中 静力平衡方程 ΣFx=0, ΣFy=0 示意图；用于物理受力平衡分析的工具。 是随机成分的向量；静力平衡；方程 ΣFx=0；示意图；物理概念；力的平衡 是不同样地间变异的方差-协方差矩阵；静力平衡，ΣFx=0，示意图，用向量分量说明受力平衡。 是作用于互补观测值上的随机成分的设计矩阵；平衡方程，ΣFx=0，ΣFy=0，力与力矩示意图，静态分析 是残差向量，其分量由胸高断面积增量与使用固定效应模型预测的增量之间的差值给出。
  3. 使用以下三个统计指标评估并比较基础模型和线性混合效应模型的预测能力23,33。
    R² 公式，统计分析，方程示意图，模型拟合评估，数据相关性评价。   (4)
    偏差公式，\(\text{Bias} = \frac{\Sigma_{i=1}^n |(ob_i-est_i)|}{N}\)，统计分析。   (5)
    RMSE 方程公式；统计误差分析；数据预测准确性度量。   (6)
    其中 obji 是观测的胸高断面积增量，esti 是预测的胸高断面积增量，静力平衡方程 Σobj；用于平衡分析的示意图 是观测值的均值，N 是观测值的数量。

访问受限。请登录或开始试用以查看此内容。

结果

粗略林分面积增量模型针对 P. asperata 表示为公式（7）。参数估计值、相应的标准误以及拟合不足统计量见表2。残差图见图1。残差表现出明显的异方差性。
森林生长方程，公式；测树参数分析。   (7)

估计值标准误差t检验P值VIF

访问受限。请登录或开始试用以查看此内容。

讨论

构建混合效应模型时，一个关键问题是如何确定哪些参数可作为随机效应，哪些应作为固定效应34,35。目前已有两种方法被提出。最常用的方法是将所有参数均视为随机效应，然后通过AIC、BIC、对数似然（Loglik）和似然比检验（LRT）选择最优模型。本研究采用的正是这一方法35。另一种方法是针对每个样本样地，使用普通最小二乘法（OLS）回归拟合胸径增量模型；在这些模型中，若某些参数在不同样地间的变异性较高，且其置信区间重叠较少，则可将这些参数视为随机效应17。

为了校正异方差性和自相关性，引入了三种方差函数和三种自相关结构。与 Calama 和 Montero17 以及 Uzoh 和 Oliver27 的研究结果一致，指数函数和一阶自回归模型 AR(1) 分别被确定为最优的方差函数和...

访问受限。请登录或开始试用以查看此内容。

披露

作者无任何利益冲突需要披露。

致谢

本研究由中央高校基本科研业务费专项资金资助，项目编号 2019GJZL04。感谢中国国家林业和草原局森林资源监测与规划研究院曾伟生教授提供数据支持。

访问受限。请登录或开始试用以查看此内容。

材料

本文使用的材料清单
姓名公司目录编号评论
计算机acer
Microsoft Office 2013
R x64 3.5.1

参考文献

  1. Meng, J., Lu, Y., Ji, Z. Transformation of a Degraded Pinus massoniana Plantation into a Mixed-Species Irregular Forest: Impacts on Stand Structure and Growth in Southern China. Forests. 5 (12), 3199-3221 (2014).
  2. Sharma, A., Bohn, K., Jose, S., Cropper, W. P. Converting even-aged plantations to uneven-aged stand conditions: A simulation analysis of silvicultural regimes with slash pine (Pinus elliottii Engelm). Forest Science. 60 (5), 893-906 (2014).
  3. Zhu, J., et al. Feasibility of implementing thinning in even-aged Larix olgensis plantations to develop uneven-aged larch–broadleaved mixed forests. Journal of Forest Research. 15 (1), 71-80 (2010).
  4. Leites, L. P., Robinson, A. P., Crookston, N. L. Accuracy and equivalence testing of crown ratio models and assessment of their impact on diameter growth and basal area increment predictions of two variants of the Forest Vegetation Simulator. Canadian Journal of Forest Research. 39 (3), 655-665 (2009).
  5. Pretzsch, H. Forest Dynamics, Growth and Yield. , (2009).
  6. Weiskittel, A. R., et al. Forest growth and yield modeling. Forest Growth & Yield Modeling. 7 (2), 223-233 (2002).
  7. Burkhart, H. E., Tomé, M. Modeling Forest Trees and Stands. , Springer. Netherlands. (2012).
  8. Zhang, X. Chinese Academy Of Forestry. A linkage among whole-stand model, individual-tree model and diameter-distribution model. Journal of Forest Science. 56 (56), 600-608 (2010).
  9. Peng, C. Growth and yield models for uneven-aged stands: past, present and future. Forest Ecology & Management. 132 (2), 259-279 (2000).
  10. Lhotka, J. M., Loewenstein, E. F. An individual-tree diameter growth model for managed uneven-aged oak-shortleaf pine stands in the Ozark Highlands of Missouri, USA. Forest Ecology & Management. 261 (3), 770-778 (2011).
  11. Porté, A., Bartelink, H. H. Modelling mixed forest growth: a review of models for forest management. Ecological Modelling. 150 (1), 141-188 (2002).
  12. Moses, L. E., Gale, L. C., Altmann, J. Methods for analysis of unbalanced, longitudinal, growth data. American Journal of Primatology. 28 (1), 49-59 (2010).
  13. Biging, G. S. Improved Estimates of Site Index Curves Using a Varying-Parameter Model. Forest Science. 31 (31), 248-259 (1985).
  14. Kowalchuk, R. K., Keselman, H. J. Mixed-model pairwise multiple comparisons of repeated measures means. Psychological Methods. 6 (3), 282-296 (2001).
  15. Hayes, A. F., Cai, L. Using heteroskedasticity-consistent standard error estimators in OLS regression: An introduction and software implementation. Behavior Research Methods. 39 (4), 709-722 (2007).
  16. Gutzwiller, K. J., Riffell, S. K. Using Statistical Models to Study Temporal Dynamics of Animal-Landscape Relations. , Springer. Boston, MA. (2007).
  17. Calama, R., Montero, G. Multilevel linear mixed model for tree diameter increment in stone pine (Pinus pinea): a calibrating approach. 39, (2005).
  18. Vonesh, E. F., Chinchilli, V. M. Linear and nonlinear models for the analysis of repeated measurements. Journal of Biopharmaceutical Statistics. 18 (4), 595-610 (1996).
  19. Zobel, J. M., Ek, A. R., Burk, T. E. Comparison of Forest Inventory and Analysis surveys, basal area models, and fitting methods for the aspen forest type in Minnesota. Forest Ecology & Management. 262 (2), 188-194 (2011).
  20. Sharma, M., Parton, J. Height-diameter equations for boreal tree species in Ontario using a mixed-effects modeling approach. Forest Ecology & Management. 249 (3), 187-198 (2007).
  21. Crecente-Campo, F., Tomé, M., Soares, P., Diéguez-Aranda, U. A generalized nonlinear mixed-effects height–diameter model for Eucalyptus globulus L. in northwestern Spain. Forest Ecology & Management. 259 (5), 943-952 (2010).
  22. Fu, L., Sharma, R. P., Hao, K., Tang, S. A generalized interregional nonlinear mixed-effects crown width model for Prince Rupprecht larch in northern China. Forest Ecology & Management. 389 (2017), 364-373 (2017).
  23. Hao, X., Yujun, S., Xinjie, W., Jin, W., Yao, F. Linear mixed-effects models to describe individual tree crown width for China-fir in Fujian Province, southeast China. Plos One. 10 (4), 0122257(2015).
  24. Vanderschaaf, C. L., Burkhart, H. E. Comparing methods to estimate Reineke's Maximum Size-Density Relationship species boundary line slope. Forest Science. 53 (3), 435-442 (2007).
  25. Zhang, L., Bi, H., Gove, J. H., Heath, L. S. A comparison of alternative methods for estimating the self-thinning boundary line. Canadian Journal of Forest Research. 35 (6), 1507-1514 (2005).
  26. Hart, D. R., Chute, A. S. Estimating von Bertalanffy growth parameters from growth increment data using a linear mixed-effects model, with an application to the sea scallop Placopecten magellanicus. Ices Journal of Marine Science. 66 (9), 2165-2175 (2009).
  27. Uzoh, F. C. C., Oliver, W. W. Individual tree diameter increment model for managed even-aged stands of ponderosa pine throughout the western United States using a multilevel linear mixed effects model. Forest Ecology & Management. 256 (3), 438-445 (2008).
  28. Condés, S., Sterba, H. Comparing an individual tree growth model for Pinus halepensis Mill. in the Spanish region of Murcia with yield tables gained from the same area. European Journal of Forest Research. 127 (3), 253-261 (2008).
  29. Pokharel, B., Dech, J. P. Mixed-effects basal area increment models for tree species in the boreal forest of Ontario, Canada using an ecological land classification approach to incorporate site effects. Forestry. 85 (2), 255-270 (2012).
  30. Wykoff, W. R. A basal area increment model for individual conifers in the northern Rocky Mountains. Forest Science. 36 (4), 1077-1104 (1990).
  31. Stage, A. R. Notes: An Expression for the Effect of Aspect, Slope, and Habitat Type on Tree Growth. Forest Science. 22 (4), 457-460 (1976).
  32. Gregorie, T. G. Generalized Error Structure for Forestry Yield Models. Forest Science. 33 (2), 423-444 (1987).
  33. Zhao, L., Li, C., Tang, S. Individual-tree diameter growth model for fir plantations based on multi-level linear mixed effects models across southeast China. Journal of Forest Research. 18 (4), 305-315 (2013).
  34. Hall, D. B., Bailey, R. L. Modeling and Prediction of Forest Growth Variables Based on Multilevel Nonlinear Mixed Models. Forest Science. 47 (3), 311-321 (2001).
  35. Yang, Y., Huang, S., Meng, S. X., Trincado, G., Vanderschaaf, C. L. A multilevel individual tree basal area increment model for aspen in boreal mixedwood stands : Journal canadien de la recherche forestière. Revue Canadienne De Recherche Forestière. 39 (39), 2203-2214 (2009).
  36. Pinheiro, J. C., Bates, D. M. Mixed-effects models in S and S-Plus. Publications of the American Statistical Association. 96 (455), 1135-1136 (2000).

访问受限。请登录或开始试用以查看此内容。

重印与许可

标签

线性混合效应模型随机区组效应方差函数自相关结构赤池信息准则贝叶斯信息准则限制性最大似然估计异方差校正一阶自回归
基于线性混合效应方法的单木胸高断面积生长量模型构建
联系我们
向图书馆推荐

隐私政策

使用条款

政策

版权所有 © 2026 MyJoVE Corporation。保留所有权利。