高级检索

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于Possion回归混合效应模型的长白落叶松一级枝数量模拟

王曼霖 董利虎 李凤日

王曼霖, 董利虎, 李凤日. 基于Possion回归混合效应模型的长白落叶松一级枝数量模拟[J]. 北京林业大学学报, 2017, 39(11): 45-55. doi: 10.13332/j.1000-1522.20170204
引用本文: 王曼霖, 董利虎, 李凤日. 基于Possion回归混合效应模型的长白落叶松一级枝数量模拟[J]. 北京林业大学学报, 2017, 39(11): 45-55. doi: 10.13332/j.1000-1522.20170204
WANG Man-lin, DONG Li-hu, LI Feng-ri. First-order branch number simulation for Larix olgensis plantation through Poisson regression mixed effect model[J]. Journal of Beijing Forestry University, 2017, 39(11): 45-55. doi: 10.13332/j.1000-1522.20170204
Citation: WANG Man-lin, DONG Li-hu, LI Feng-ri. First-order branch number simulation for Larix olgensis plantation through Poisson regression mixed effect model[J]. Journal of Beijing Forestry University, 2017, 39(11): 45-55. doi: 10.13332/j.1000-1522.20170204

基于Possion回归混合效应模型的长白落叶松一级枝数量模拟

doi: 10.13332/j.1000-1522.20170204
基金项目: 

“十三五”国家重点研发计划课题 2017YFD0600402

详细信息
    作者简介:

    王曼霖。主要研究方向:林分生长与收获模型。Email: 445947011@qq.com   地址: 150040 黑龙江省哈尔滨市和兴路26号东北林业大学林学院

    通讯作者:

    李凤日,教授,博士生导师。主要研究方向:林分生长与收获模型。Email: fengrili@126.com   地址:同上

  • 中图分类号: S758.5

First-order branch number simulation for Larix olgensis plantation through Poisson regression mixed effect model

  • 摘要: 利用广义线性混合模型对长白落叶松一级枝条数量进行研究,以黑龙江省佳木斯市孟家岗林场长白落叶松人工林为研究对象,基于7块标准地49株枝解析样木的596个一级枝条测定数据,利用SAS 9.3软件中的PROC GLIMMIX模块,建立了基于Poisson分布的一级枝条数量的最优基础模型。在此基础上考虑树木效应,构建每半米段一级枝条数量的广义线性混合模型,并利用AIC、BIC、-2log likelihood以及LRT检验对收敛模型的拟合优度进行比较。结果表明:任意参数组合的混合效应模型的拟合效果均好于传统模型,最终将含有DINC、LnRDINC、RDINC2这3个随机效应参数的模型作为长白落叶松每半米段一级枝条数量分布的最优混合效应模型。模型拟合结果显示,LnRDINC、CL的参数估计值为正值,DINC、RDINC2、HT/DBH、DBH的参数估计值为负值,每半米段一级枝条分布数量在树冠范围内存在峰值,模型的确定系数R2为0.669,拟合的平均绝对误差为2.250,均方根误差为3.012。从总体上看,所建立的一级枝条分布数量混合模型不但可以反映总体枝条数量的变化趋势,还可以反映树木之间的个体差异,说明广义线性混合模型确实可以提高模型的模拟精度。所得出的混合模型可以很好地预估该研究区内人工长白落叶松每半米段一级枝条数量的分布情况,为定量研究长白落叶松树冠构筑型和三维可视化模拟提供了基础。
  • 图  1  人工长白落叶松每半米段一级枝条数量分布

    高径比、冠长模拟值来自胸径为13.5 cm的人工长白落叶松解析木。Simulated values of HT/DBH,CL based on a Larix olgensis with DBH=13.5 cm.

    Figure  1.  Predicted value of number of first-order branches per 0.5 m

    图  2  人工长白落叶松一级枝条数量的实测值与预测值散点图

    Figure  2.  Observed values and predicted values of first-order branch for Larix olgensis planation

    图  3  人工长白落叶松相对冠层深度与每半米段的预测枝条数量关系

    Figure  3.  Relationship between relative distance from tree apex and the number of predicted branches per 0.5 m section for Larix olgensis planation

    表  1  长白落叶松人工林林分因子统计表

    Table  1.   Statistics of stand variables for Larix olgensis planation

    项目
    Item
    年龄/a
    Age/year
    平均胸径
    Mean DBH/cm
    平均树高
    Mean tree height/m
    平均冠幅
    Mean crown width/m
    林分密度/(株·hm-2)
    Stand density/(tree·ha-1)
    坡度
    Slope/(°)
    海拔
    Altitude/m
    最小值Minimum 9 5.40 6.20 1.10 1 016.70 5.00 141.80
    最大值Maximum 33 19.90 24.10 1.60 3 083.30 12.50 260.10
    平均值Mean 19.19 11.30 13.31 1.27 2 223.80 6.30 193.11
    标准差Std 7.53 4.19 5.20 0.19 579.72 2.90 35.35
    变异系数CV/% 39.24 37.08 39.07 14.96 26.07 46.03 18.31
    下载: 导出CSV

    表  2  长白落叶松人工林解析木和枝条分布统计表

    Table  2.   Statistics of sample trees and branch distribution for Larix olgensis planation

    项目
    Item
    胸径
    DBH/cm
    树高
    Tree height
    (HT)/cm
    冠长
    Crown length
    (CL)/m
    冠幅
    Crown width
    (CW)/m
    高径比
    HT/DBH
    每半米段的一级枝个数
    Number of first-order
    branch per 0.5 m
    建模数据(样本数量=40)
    Fitting data(sample size=40)
    最小值Minimum 2.00 3.80 2.30 0.55 0.66 1
    最大值Maximum 27.00 21.50 11.50 3.18 1.41 26
    平均值Mean 13.68 12.45 6.24 1.42 0.85 7.56
    标准差Std 5.42 4.45 2.38 0.56 0.26 4.21
    变异系数CV/% 39.61 35.74 38.14 39.44 30.59 55.69
    检验数据(样本数量=9)
    Validation data(sample size=9)
    最小值Minimum 3.80 5.30 2.52 0.51 0.70 1
    最大值Maximum 25.20 21.30 10.51 3.20 1.39 23
    平均值Mean 14.03 12.86 6.01 1.49 0.81 7.83
    标准差Std 5.88 4.13 2.49 0.59 0.22 3.86
    变异系数CV/% 37.56 32.12 41.43 39.60 27.16 49.30
    下载: 导出CSV

    表  3  长白落叶松一级枝条数量模型所需变量及相关描述

    Table  3.   Symbol and associated description for variables tested in the first-order branch model for Larix olgensis planation

    变量符号
    Variable symbol
    变量描述
    Description
    DINC 冠层深度Distance from tree apex
    RDINC 相对冠层深度Relative distance from tree apex
    LnDINC 冠层深度的对数值Logarithm of distance from tree apex
    LnRDINC 相对冠层深度的对数值Logarithm of relative distance from tree apex
    RDINC2 相对冠层深度的平方Relative distance from tree apex squared
    DBH 树木胸高处的直径Diameter at breast height/cm
    HT 树高Tree height/m
    HT/DBH 树高和胸径的比值Ratio of tree height to diameter at breast height
    CL 树冠冠长Crown length/m
    下载: 导出CSV

    表  4  基于树木效应的混合模型拟合结果

    Table  4.   Fitting results of mixed model based on individual tree effects

    参数个数
    Parameter
    number
    模型
    Model
    随机参数Random parameter AIC BIC -2log
    likelihood
    截距Intercept HT/DBH DINC LnRDINC CL RDINC2 DBH
    7 (7) 无None 2 610.85 2 640.17 2 596.85
    (8) 2 600.01 2 617.92 2 580.11
    (9) 2 599.16 2 619.67 2 582.16
    (10) 2 584.90 2 598.41 2 568.90
    8 (11) 2 600.07 2 618.98 2 581.47
    (12) 2 600.19 2 616.50 2 582.99
    (13) 2 580.35 2 593.86 2 564.35
    (14) 2 599.50 2 616.01 2 586.50
    10 (15) 2 575.01 2 591.90 2 555.01
    (16) 2 571.42 2 588.31 2 551.42
    (17) 2 574.70 2 591.98 2 554.70
    (18) 2 575.65 2 590.85 2 557.65
    (19) 2 576.51 2 593.40 2 556.61
    (20) 2 576.65 2 593.54 2 556.65
    (21) 2 572.97 2 589.86 2 552.97
    (22) 2 574.16 2 591.05 2 554.16
    (23) 2 568.91 2 585.80 2 548.91
    (24) 2 573.76 2 590.65 2 553.76
    (25) 2 573.08 2 589.97 2 553.08
    (26) 2 574.18 2 591.06 2 554.18
    (27) 2 583.32 2 600.20 2 563.32
    (28) 2 574.62 2 591.51 2 554.62
    (29) 2 580.54 2 592.43 2 561.54
    (30) 2 572.95 2 589.84 2 552.95
    13 25(31) 2 568.95 2 590.91 2 542.95
    26(32) 2 566.07 2 588.03 2 540.07
    27(33) 2 570.77 2 592.73 2 544.77
    28(34) 2 568.37 2 590.32 2 542.37
    29(35) 2 565.65 2 587.61 2 539.65
    30(36) 2 572.20 2 594.15 2 546.20
    31(37) 2 569.04 2 590.99 2 543.04
    32(38) 2 567.45 2 589.40 2 541.45
    33(39) 2 565.74 2 587.70 2 539.74
    34(40) 2 568.37 2 590.32 2 542.37
    注:参数个数包括固定效应参数个数、随机效应方差-协方差参数个数;带★变量表示包含此随机参数;表中的模型(7)即正文中的公式(7),即不含随机效应参数的最优基础模型。Notes: the number of parameters includes the fixed parameters, parameters in covariance-structure G for random effect; variables with ★ are included in the random parameters; model (7) in the table is the formula (7) model, which is the optimal basic model without random effect parameters.
    下载: 导出CSV

    表  5  不同参数个数的模型固定参数和方差组成的估计值

    Table  5.   Fixed parameters and variance component estimates of models with different numbers of parameters

    固定参数
    Fixed parameter
    模型(7)
    Model (7)
    模型(13)
    Model (13)
    模型(23)
    Model (23)
    模型(35)
    Model (35)
    截距Int 3.058 2 2.941 7 3.178 4 3.086 6
    DINC -0.002 1 -0.002 1 -0.002 2 -0.002 5
    LnRDINC 0.342 3 0.361 7 0.370 9 0.415 3
    RDINC2 -1.044 4 -1.153 0 -1.129 3 -1.046 1
    HT/DBH -0.282 2 -0.147 1 -0.280 3 -0.245 4
    DBH -0.015 3 -0.004 34 -0.007 4 -0.011 2
    CL 0.094 9 0.080 4 0.076 7 0.105 0
    随机效应方差-协方差结构
    Variance of random effect-covariance
    structure
    [0.279 1] $\left[ {\begin{array}{*{20}{c}} {0.000\;01}&{ - 0.000\;99}\\ { - 0.000\;99}&{1.446\;3} \end{array}} \right]$ $ \left[ {\begin{array}{*{20}{c}} {0.000\;001}&{0.000\;033}&{ - 0.001\;16}\\ {0.000\;033}&{0.004\;852}&{ - 0.011\;22}\\ { - 0.001\;16}&{ - 0.011\;22}&{1.565\;0} \end{array}} \right]$
    LRT 32.5 15.44 9.26
    P <0.001 <0.001 0.026
    Pearson χ2/自由度 1.73 1.50 1.38 1.20
    Pearson χ2/degree
    下载: 导出CSV

    表  6  基于随机参数不同方差-协方差结构的混合模型模拟结果比较

    Table  6.   Comparisons of mixed model based on different variance-covariance structure of random parameters

    方差-协方差结构
    Variance-covariance structure
    参数个数
    Parameter number
    AIC BIC -2log likelihood LRT P
    复合对称矩阵Composite symmetric matrix (CS) 9 2 586.71 2 601.91 2 579.71
    对角矩阵Diagonal matrix (UN(1)) 13 2 580.99 2 597.87 2 560.99 18.72 < 0.001
    无结构矩阵Non-structural matrix (UN) 13 2 565.65 2 587.61 2 539.65 21.34
    下载: 导出CSV

    表  7  检验数据的随机参数估计结果

    Table  7.   Random parameter estimation result of validation data

    检验数据
    Validation data
    随机参数
    Random parameter
    DINC LnRDINC RDINC2
    b1 -0.000 5 -0.013 2 0.383 1
    b2 0.000 6 -0.022 6 -0.971 2
    b3 0.000 9 0.014 7 -1.068 5
    b4 0.000 8 -0.073 4 -1.525 0
    b5 0.000 1 -0.109 2 -1.883 0
    b6 0.000 7 0.044 0 -0.629 0
    b7 -0.001 6 -0.114 4 1.453 0
    b8 -0.001 6 -0.173 4 1.109 1
    b9 0.001 4 0.120 6 -1.083 0
    下载: 导出CSV

    表  8  基础模型和混合效应模型模拟结果比较

    Table  8.   Comparison of based model and mixed-effect model

    评价指标
    Evaluation index
    固定模型
    Fixed model
    混合模型
    Mixed model
    确定系数
    Determination coefficient (R2)
    0.510 0.669
    均方根误差
    Root mean square error (RMSE)
    3.396 3.012
    平均绝对误差
    Average absolute error (MAE)
    2.601 2.250
    下载: 导出CSV
  • [1] BROWN P L, DOLEY D, KEENAN R J. Stem and crown dimensions as predictors of thinning responses in a crowded tropical rainforest plantation of Flindersia brayleyana F. Muell[J]. Forest Ecology & Management, 2004, 196(2/3): 379-392. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=a09870ffce90b9fdf1ec2a36d7723f87
    [2] CRECENTE-CAMPO F, MARSHALL P, LEMAY V, et al. A crown profile model for Pinus radiata D. Don. in northwestern Spain[J]. Forest Ecology & Management, 2009, 257(12): 2370-2379. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=12965a1d48a157ea8ca6be89351acf79
    [3] KADANE J B, KRISHNAN R, SHMUELI G. A data disclosure policy for count data based on the Com-Poisson distribution[J]. Management Science, 2006, 52(10): 1610-1617. doi:  10.1287/mnsc.1060.0562
    [4] ZHANG H, XU J, JIANG N, et al. PLNseq: a multivariate Poisson lognormal distribution for high-throughput matched RNA-sequencing read count data[J]. Statistics in Medicine, 2015, 34(9): 1577-1589. doi:  10.1002/sim.6449
    [5] CHOOWOSOBA H, LEVY S M, DATTA S. Marginal regression models for clustered count data based on zero-inflated Conway-Maxwell-Poisson distribution with applications[J]. Biometrics, 2015, 72(2): 606-618. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=cc718cf3ede24b6a3facab73c180678c
    [6] PANDHARD X, SAMSON A. Extension of the SAEM algorithm for nonlinear mixed models with 2 levels of random effects[J]. Biostatistics, 2009, 10(1): 121-135. http://d.old.wanfangdata.com.cn/OAPaper/oai_arXiv.org_0803.4437
    [7] MENG S X, HUANG S. Improved calibration of nonlinear mixed-effects models demonstrated on a height growth function[J]. Forest Science, 2009, 55(3): 239-248. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=39ee357cfca94197848e4b1076473535
    [8] YANG Y, HUANG S. Comparison of different methods for fitting nonlinear mixed forest models and for making predictions[J]. Canadian Journal of Forest Research, 2011, 41(8): 1671-1686. doi:  10.1139/x11-071
    [9] HEIN S, MÄKINEN H, YUE C, et al. Modelling branch characteristics of Norway spruce from wide spacings in Germany[J]. Forest Ecology & Management, 2007, 242(2-3): 155-164. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=f8db8e620a46053051843fccd9758aee
    [10] HEIN S, WEISKITTEL A R, KOHNLE U. Branch characteristics of widely spaced Douglas-fir in south-western Germany: comparisons of modelling approaches and geographic regions[J]. Forest Ecology & Management, 2008, 256(5): 1064-1079. http://cn.bing.com/academic/profile?id=8cacef1f84fe8e52dd5dab7163b01013&encoded=0&v=paper_preview&mkt=zh-cn
    [11] BEAULIEU E, SCHNEIDER R, BERNINGER F, et al. Modeling jack pine branch characteristics in eastern Canada[J]. Fuel & Energy Abstracts, 2011, 262(9): 1748-1757. http://cn.bing.com/academic/profile?id=f68060386f408f209c6342f7a3b9c252&encoded=0&v=paper_preview&mkt=zh-cn
    [12] LINNELL N A F, PARISHROBERTA, GOUDIEJAMES W. Modelling number, vertical distribution, and size of live branches on coniferous tree species in British Columbia[J]. Canadian Journal of Forest Research, 2012, 42(42): 1072-1090. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=577c6f9204a93dac7b1af00e69fb9168
    [13] KINT V, HEIN S, CAMPIOLI M, et al. Modelling self-pruning and branch attributes for young Quercus robur L. and Fagus sylvatica L. trees[J]. Fuel & Energy Abstracts, 2010, 260(11): 2023-2034. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=860f57b85f95cc5709d4708551c1bf08
    [14] NEMEC A F L, GOUDIE J W, PARISH R. A Gamma-Poisson model for vertical location and frequency of buds on lodgepole pine (Pinus contorta) leaders[J]. Canadian Journal of Forest Research, 2009, 40(40): 2049-2058. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=21370f02b2d0389d23a208619b29b07a
    [15] SATTLER D F, COMEAU P G, ACHIM A. Branch models for white spruce (Picea glauca (Moench) Voss) in naturally regenerated stands[J]. Forest Ecology & Management, 2014, 325: 74-89. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=0e350c83cbe6d90bc1980a067f5b7303
    [16] 郭孝玉.长白落叶松人工林树冠结构及生长模型研究[D].北京: 北京林业大学, 2013. http://cdmd.cnki.com.cn/Article/CDMD-10022-1013215124.htm

    GUO X Y. Crown structure and growth modle for Larix algersis plantation[D]. Beijing: Beijing Forestry University, 2013. http://cdmd.cnki.com.cn/Article/CDMD-10022-1013215124.htm
    [17] 郑杨, 董利虎, 李凤日.黑龙江省红松人工林枝条分布数量模拟[J].应用生态学报, 2016, 27(7):2172-2180. http://d.old.wanfangdata.com.cn/Periodical/yystxb201607016

    ZHENG Y, DONG L H, LI F R. Branch quantity distribution simulation for Pinus koraiensis plantation in Heilongjiang Province, China[J]. Chinese Journal of Applied Ecology, 2016, 27(7):2172-2180. http://d.old.wanfangdata.com.cn/Periodical/yystxb201607016
    [18] BURNHAM K P, ANDERSON D R, HUYVAERT K P. AIC model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons[J]. Behavioral Ecology & Sociobiology, 2011, 65(1): 23-35. doi:  10.1007-s00265-010-1029-6/
    [19] AHO K, DERRYBERRY D W, PETERSON T. Model selection for ecologists: the worldviews of AIC and BIC[J]. Ecology, 2014, 95(3): 631-636. doi:  10.1890/13-1452.1
    [20] WATSON S P. Poisson regression models for interval censored count data[D/OL].Waco: Baylor University, 2011[2017-03-10].http://gradworks.umi.com/34/56/3456978/htmlhttps://baylor-ir.tdl.org/handle/2104/8176
    [21] ZOU G Y, DONNER A. Extension of the modified Poisson regression model to prospective studies with correlated binary data[J]. Statistical Methods in Medical Research, 2011, 22(6): 661-670. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=10.1177/0962280211427759
    [22] GHITANY M E, KARLIS D, AL-MUTAIRI D K, et al. An EM algorithm for multivariate mixed Poisson regression models and its application[J]. Applied Mathematical Sciences, 2012, 6:6843-6856. http://d.old.wanfangdata.com.cn/OAPaper/oai_doaj-articles_4980ae0ae2ec78516b3df003c08fdb15
    [23] ZAMANI H, FAROUGHI P, ISMAIL N. Estimation of count data using mixed Poisson, generalized Poisson and finite Poisson mixture regression models[J]. AIP Conference Proceedings, 2014, 1602(1): 1144-1150. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=CC0214576928
    [24] LONG D L, PREISSER J S, HERRING A H, et al. A marginalized zero-inflated Poisson regression model with random effects[J]. Journal of the Royal Statistical Society, 2015, 64(5): 815-830. doi:  10.1111/rssc.12104
    [25] PINHEIRO J C, BATES D M. Mixed-effects models in S and S-PLUS[J]. Journal of the American Statistical Association, 2000, 96(455): 1135-1136. http://d.old.wanfangdata.com.cn/OAPaper/oai_doaj-articles_024659aa5aed0dcdbb679c9e020372cc
    [26] BOND B J. Age-related changes in photosynthesis of woody plants[J]. Trends in Plant Science, 2000, 5(8): 349-353. doi:  10.1016/S1360-1385(00)01691-5
    [27] WEISKITTEL A R, MAGUIRE D A, MONSERUD R A. Response of branch growth and mortality to silvicultural treatments in coastal Douglas-fir plantations: implications for predicting tree growth[J]. Forest Ecology & Management, 2007, 251(3): 182-194. http://cn.bing.com/academic/profile?id=af31a4b4e00ff12735a2449538d8c0d5&encoded=0&v=paper_preview&mkt=zh-cn
    [28] WEISKITTEL A R, SEYMOUR R S, HOFMEYER P V. Modelling primary branch frequency and size for five conifer species in Maine, USA[J]. Forest Ecology & Management, 2010, 259(10): 1912-1921. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=6525a57d6fbb2148c9c7838c5bf3a336
    [29] ISHⅡ H, MCDOWELL N. Age-related development of crown structure in coastal Douglas-fir trees[J]. Forest Ecology and Management, 2002, 169(3): 257-270. doi:  10.1016/S0378-1127(01)00751-4
    [30] THORPE H C, ASTRUP R, TROWBRIDGE A, et al. Competition and tree crowns: a neighborhood analysis of three boreal tree species[J]. Forest Ecology & Management, 2010, 259(8): 1586-1596. http://cn.bing.com/academic/profile?id=e01e55e9ca64ee12f100228cbeb0efbf&encoded=0&v=paper_preview&mkt=zh-cn
  • [1] 贺梦莹, 董利虎, 李凤日.  长白落叶松−水曲柳混交林冠幅预测模型 . 北京林业大学学报, 2020, 42(7): 23-32. doi: 10.12171/j.1000-1522.20190250
    [2] 牛亦龙, 董利虎, 李凤日.  基于广义代数差分法的长白落叶松人工林地位指数模型 . 北京林业大学学报, 2020, 42(2): 9-18. doi: 10.12171/j.1000-1522.20190036
    [3] 靳晓娟, 孙玉军, 潘磊.  基于混合效应的长白落叶松一级枝条基径预估模型 . 北京林业大学学报, 2020, 42(): 1-10. doi: 10.12171/j.1000-1522.20200133
    [4] 贾炜玮, 冯万举, 李凤日.  基于节子剖析数据的长白落叶松人工林枝条丢失年轮数研究 . 北京林业大学学报, 2020, 42(3): 87-98. doi: 10.12171/j.1000-1522.20190038
    [5] 白东雪, 刘强, 董利虎, 李凤日.  长白落叶松人工林有效冠高的确定及其影响因子 . 北京林业大学学报, 2019, 41(5): 76-87. doi: 10.13332/j.1000-1522.20190016
    [6] 徐奇刚, 雷相东, 国红, 李海奎, 李玉堂.  基于多层感知机的长白落叶松人工林林分生物量模型 . 北京林业大学学报, 2019, 41(5): 97-107. doi: 10.13332/j.1000-1522.20190035
    [7] 沈剑波, 雷相东, 雷渊才, 李玉堂.  长白落叶松人工林地位指数及立地形的比较研究 . 北京林业大学学报, 2018, 40(6): 1-8. doi: 10.13332/j.1000-1522.20170400
    [8] 王涛, 董利虎, 李凤日.  基于混合效应的杂种落叶松人工幼龄林单木枯损模型 . 北京林业大学学报, 2018, 40(10): 1-10. doi: 10.13332/j.1000-1522.20170437
    [9] 王烁, 董利虎, 李凤日.  人工长白落叶松枝条存活模型 . 北京林业大学学报, 2018, 40(1): 57-66. doi: 10.13332/j.1000-1522.20170203
    [10] 王冬至, 张冬燕, 张志东, 黄选瑞.  塞罕坝华北落叶松人工林断面积预测模型 . 北京林业大学学报, 2017, 39(7): 10-17. doi: 10.13332/j.1000-1522.20170072
    [11] 邵英男, 田松岩, 刘延坤, 陈瑶, 孙志虎.  密度调控对长白落叶松人工林土壤呼吸的影响 . 北京林业大学学报, 2017, 39(6): 51-59. doi: 10.13332/j.1000-1522.20170029
    [12] 陈东升, 孙晓梅, 李凤日, 贾炜玮.  落叶松人工林节子内部特征变化规律研究 . 北京林业大学学报, 2015, 37(2): 16-23. doi: 10.13332/j.cnki.jbfu.2015.02.014
    [13] 高慧淋, 李凤日, 董利虎.  基于分段回归的人工红松冠形预估模型 . 北京林业大学学报, 2015, 37(3): 76-83. doi: 10.13332/j.1000-1522.20140324
    [14] 李丹, 戴巍, 闫志刚, 王霓虹.  基于模糊层次分析法的落叶松人工林生境评价系统 . 北京林业大学学报, 2014, 36(4): 75-81. doi: 10.13332/j.cnki.jbfu.2014.04.015
    [15] 张建亮, 崔国发, 黄祥童, 郭子良, 周海城.  国家一级保护植物长白松种群结构与动态预测 . 北京林业大学学报, 2014, 36(3): 26-33. doi: 10.13332/j.cnki.jbfu.2014.03.004
    [16] 刘艳红, 马炜.  长白落叶松人工林可燃物碳储量分布及燃烧性 . 北京林业大学学报, 2013, 35(3): 32-38.
    [17] 李耀翔, 姜立春.  基于非线性混合模型的落叶松木材管胞长度模拟 . 北京林业大学学报, 2013, 35(3): 18-23.
    [18]
    孙志虎, 毕永娟, 牟长城, 蔡体久
    基于FORECAST模型的长白落叶松人工林经营措施对长期生产力的影响 . 北京林业大学学报, 2012, 34(6): 1-6.
    [19] 孙志虎, 牟长城, 张彦东.  用地统计学方法估算长白落叶松人工林凋落物现存量 . 北京林业大学学报, 2008, 30(4): 59-64.
    [20] 崔彬彬, 李贤军, 宗世祥, 赵俊卉, 肖化顺, 陈伟, 刘志军, 王志玲, 曹伟, 黄心渊, 张煜星, 周国模, 李国平, 江泽慧, 雷相东, 刘智, 施婷婷, 张展羽, 于寒颖, 周志强, 杜官本, 徐剑琦, 程金新, 雷霆, 程丽莉, 曹金珍, 关德新, 刘童燕, 张贵, 苏里坦, 吴家森, 骆有庆, 王正, 丁立建, 王正, 张则路, 张彩虹, 王海, 杨谦, 张璧光, 苏淑钗, 李云, 张璧光, 郭广猛, 郝雨, 黄群策, 雷洪, 李云, 张国华, 刘彤, 金晓洁], 吴家兵, 黄晓丽, 贺宏奎, 王勇, 张书香, 张慧东, 常亮, 秦岭, 方群, 秦广雍, 张佳蕊, 许志春, 张大红, 陈晓光, 宋南, 刘大鹏, 姜培坤, 李文军, 周晓燕, 李延军, 高黎, 刘海龙, 蔡学理, 陈燕, 姜静, 姜金仲, 张弥, 冯慧, 苏晓华, 于兴华, 张金桐, 刘建立, 王安志, 张冰玉, 尹伟伦, 陈绪和, 周梅, 王谦, 朱彩霞, 成小芳, 王德国, 陈建伟3, 聂立水, 亢新刚, 张连生, 张勤, 冯大领, 金昌杰, 梁树军, 崔国发, 韩士杰, 胡君艳, 姚国龙.  长白落叶松等几个树种冠幅预测模型的研究 . 北京林业大学学报, 2006, 28(6): 75-79.
  • 加载中
图(3) / 表 (8)
计量
  • 文章访问数:  1670
  • HTML全文浏览量:  189
  • PDF下载量:  43
  • 被引次数: 0
出版历程
  • 收稿日期:  2017-06-12
  • 修回日期:  2017-10-19
  • 刊出日期:  2017-11-01

基于Possion回归混合效应模型的长白落叶松一级枝数量模拟

doi: 10.13332/j.1000-1522.20170204
    基金项目:

    “十三五”国家重点研发计划课题 2017YFD0600402

    作者简介:

    王曼霖。主要研究方向:林分生长与收获模型。Email: 445947011@qq.com   地址: 150040 黑龙江省哈尔滨市和兴路26号东北林业大学林学院

    通讯作者: 李凤日,教授,博士生导师。主要研究方向:林分生长与收获模型。Email: fengrili@126.com   地址:同上
  • 中图分类号: S758.5

摘要: 利用广义线性混合模型对长白落叶松一级枝条数量进行研究,以黑龙江省佳木斯市孟家岗林场长白落叶松人工林为研究对象,基于7块标准地49株枝解析样木的596个一级枝条测定数据,利用SAS 9.3软件中的PROC GLIMMIX模块,建立了基于Poisson分布的一级枝条数量的最优基础模型。在此基础上考虑树木效应,构建每半米段一级枝条数量的广义线性混合模型,并利用AIC、BIC、-2log likelihood以及LRT检验对收敛模型的拟合优度进行比较。结果表明:任意参数组合的混合效应模型的拟合效果均好于传统模型,最终将含有DINC、LnRDINC、RDINC2这3个随机效应参数的模型作为长白落叶松每半米段一级枝条数量分布的最优混合效应模型。模型拟合结果显示,LnRDINC、CL的参数估计值为正值,DINC、RDINC2、HT/DBH、DBH的参数估计值为负值,每半米段一级枝条分布数量在树冠范围内存在峰值,模型的确定系数R2为0.669,拟合的平均绝对误差为2.250,均方根误差为3.012。从总体上看,所建立的一级枝条分布数量混合模型不但可以反映总体枝条数量的变化趋势,还可以反映树木之间的个体差异,说明广义线性混合模型确实可以提高模型的模拟精度。所得出的混合模型可以很好地预估该研究区内人工长白落叶松每半米段一级枝条数量的分布情况,为定量研究长白落叶松树冠构筑型和三维可视化模拟提供了基础。

English Abstract

王曼霖, 董利虎, 李凤日. 基于Possion回归混合效应模型的长白落叶松一级枝数量模拟[J]. 北京林业大学学报, 2017, 39(11): 45-55. doi: 10.13332/j.1000-1522.20170204
引用本文: 王曼霖, 董利虎, 李凤日. 基于Possion回归混合效应模型的长白落叶松一级枝数量模拟[J]. 北京林业大学学报, 2017, 39(11): 45-55. doi: 10.13332/j.1000-1522.20170204
WANG Man-lin, DONG Li-hu, LI Feng-ri. First-order branch number simulation for Larix olgensis plantation through Poisson regression mixed effect model[J]. Journal of Beijing Forestry University, 2017, 39(11): 45-55. doi: 10.13332/j.1000-1522.20170204
Citation: WANG Man-lin, DONG Li-hu, LI Feng-ri. First-order branch number simulation for Larix olgensis plantation through Poisson regression mixed effect model[J]. Journal of Beijing Forestry University, 2017, 39(11): 45-55. doi: 10.13332/j.1000-1522.20170204
  • 树冠结构不仅影响树木对光环境的利用和水资源的获取,而且对木材质量的形成也有直接的影响[1-2]。枝条作为叶片的支撑结构,它的数量、分布决定了树冠的形状、大小,直接影响树木的光合作用效率及木材质量。因此,深入研究枝条数量的分布规律可以更透彻地理解树冠的动态变化以及树木的生长机制。

    长白落叶松(Larix olgensis)枝条数量模型可归类为计数模型,而Poisson回归模型是专门为计数变量而发展起来的统计模型,它不仅广泛应用于林业研究,还可应用于生物学、医学等更广泛的领域[3-5]。然而,在林业实际调查过程中,收集的枝条特征因子数据来自不同的林木和样地,传统模型只能表现出枝条的平均特征,从而导致预测的结果精度较低;混合模型通过引入随机效应参数,利用方差-协方差结构来表达个体间的随机效应和重复观测导致的自相关结构,可以更好地反映个体间的相关性及异方差性,进而使模型的预估精度得到提高[6-8]。目前,国外学者利用混合模型对枝条分布数量的研究较为广泛。Hein等[9]以树高年生长量和高径比为预测变量,选择Poisson回归模型,运用广义线性混合模型来模拟挪威云杉(Picea abies)每轮的枝条数量。Weiskittel等[10]运用广义线性混合模型和非线性混合模型对花旗松(Pseudotsuga menziesii)每轮的枝条数量进行模拟。Beaulieu等[11]运用广义线性混合泊松模型,以枝长和胸径生长量为预测变量对加拿大短叶松(Pinus banksiana)每年的新生枝条数量进行模拟。研究发现林木个体可能受到遗传性的控制,使个体间的随机效应对枝条的数量影响较大,而样地间的随机效应对枝条的数量影响不大。Amanda等[12]运用非齐次复合Poisson模型对4种针叶树枝条的数量、垂直分布以及大小进行模拟。Kint等[13]以Poisson分布为基础模型,对夏栎(Quercus robur)和欧洲山毛榉(Fagus sylvatica)的枝条数量建立林木水平和枝条水平的广义线性混合模型。Nemec等[14]建立了基于Gamma-Poisson分布的马尾松(Pinus massoniana)枝条数量模型,研究发现每轮内的枝条数量与树高年生长量和第一活枝高这2个因子密切相关。Sattler等[15]选择伴随log链接函数的Poisson回归模型来研究白云杉(Picea glauca)整米段内的枝条数量,建立树木水平和样地水平的混合模型。由于轮生枝条的位置不明显,因而采用整米段切割的方式来预估枝条的特征。以上研究均表明了混合效应模型在计数数据方面有很好的拟合效果。国内学者郭孝玉[16]采用冠层分析法研究了长白落叶松全树的枝条数量和每米段内的枝条数量变化规律,发现总枝条数量与胸径密切相关,每米段内的枝条数量与枝条冠层深度、高径比密切相关。郑杨等[17]建立了基于Poisson回归和负二项回归的红松(Pinus koraiensis)二级枝条数量模型,研究表明随着相对着枝深度的增加,冠内的二级枝条数量先增大后减小。因此研究枝条数量的变化对树木构筑型的预测也有一定的参考价值。

    本文以孟家岗林场49棵长白落叶松为研究对象,以Poisson回归模型为基础建立广义线性混合模型来模拟一级枝条的数量分布情况。采用SAS9.3统计软件中的PROC GLIMMIX模块,建立包含树木效应的单水平枝条数量模型。通过比较收敛模型的赤池信息量准则(AIC)、贝叶斯信息量准则(BIC)、-2倍对数的估计模型的似然函数最大值(-2log likelihood),并对模型进行LRT检验[18-19],从而选出可以模拟长白落叶松一级枝条数量的最优模型。

    • 研究地区位于黑龙江省佳木斯市孟家岗林场,地理坐标130°32′ ~130°52′ E,46°20′~46°30′ N。该地区平均海拔为250 m,变化范围为170~575 m。该区属东亚大陆性季风气候,极端最高气温为35.6 ℃,最低气温为-34.7 ℃,年平均降水量550 mm,无霜期约120 d,全年日照时数1 955 h。土壤以典型暗棕壤为主,其次是白浆化暗棕壤,还有少量的草甸土、沼泽土等。该林场地处完达山西麓余脉,以低山丘陵为主,坡度较缓。地势东北高,西南低。水系隶属于为松花江支流,最大河流为三级支流柳树河。该林场以人工针叶林为主,主要树种为长白落叶松、红松、樟子松(P. sylvestris var. mongolica)等。

    • 2015年在孟家岗林场不同立地条件、不同年龄的长白落叶松人工林中设置7块面积为0.06 hm2(20 m×30 m)的固定标准地。对每块标准地内的所有林木进行逐一编号、坐标定位、每木检尺。根据每木检尺的结果,按照林木径阶从大到小的顺序,将林木分为断面积基本相等的5个径级,也就是常说的等断面积径级标准木法;然后计算各径级林木的平均胸径,以此为标准在各径级林木中选择一棵作为解析木;再选取标准地内胸径最粗和最细的林木各6株,分别取平均值,以此为标准选出每块标准地内的优势木和劣势木,即每块标准地均选7株标准木,所有标准木均在标准地周围选取。首先测定解析木的树高、胸径、冠幅、第一活枝高等林木因子,以及相邻木的树种、树高、胸径、距离等数据,并在树干上标明北向,为测定枝条的方位角提供标准;然后将解析木伐倒,测量梢头到树冠基部的距离,即为冠长,树冠基部的位置即为第一活枝高的位置。由于长白落叶松属于轮生枝不太明显的树种(含有水枝),很难鉴别枝条的实际生长年份。为了更好地分析冠层结构,可采用树干主轴切割分层法来代替传统的以轮生枝为依据的分析方法。由于本研究只讨论活枝的数量,因此将树干的树冠部分严格以半米区分段为单位切割,从树冠基部往上把着生在每个半米区分段的枝条划为一层,不足半米区分段的部分不作研究,因而可以将树冠划分为多个冠层,并只对每层的轮生枝(除去水枝)进行测量进而分析枝条特征因子随冠层深度的变化情况。把每个半米段的顶部到梢头的距离定义为该米段的冠层深度(DINC),冠层深度与冠长的比值为该米段的相对冠层深度(RDINC),并测定每个半米段内所有一级枝条的枝长(BL)、基径(BD)、着枝角度(θ)等枝条属性因子,并标注存活状态,统计每个半米段内一级活枝条的数量,进而获得49棵长白落叶松解析木的一级活枝条数量数据。将49株解析木按照4:1分为建模数据和检验数据。长白落叶松人工林的主要林分因子详见表 1;49棵解析木的林木测定因子以及枝条数量的统计信息详见表 2。由表 2得知每半米段的一级枝条数量由1到26不等,数据呈现出较离散的状态。长白落叶松一级枝条数量模型所需的变量及相关描述见表 3

      表 1  长白落叶松人工林林分因子统计表

      Table 1.  Statistics of stand variables for Larix olgensis planation

      项目
      Item
      年龄/a
      Age/year
      平均胸径
      Mean DBH/cm
      平均树高
      Mean tree height/m
      平均冠幅
      Mean crown width/m
      林分密度/(株·hm-2)
      Stand density/(tree·ha-1)
      坡度
      Slope/(°)
      海拔
      Altitude/m
      最小值Minimum 9 5.40 6.20 1.10 1 016.70 5.00 141.80
      最大值Maximum 33 19.90 24.10 1.60 3 083.30 12.50 260.10
      平均值Mean 19.19 11.30 13.31 1.27 2 223.80 6.30 193.11
      标准差Std 7.53 4.19 5.20 0.19 579.72 2.90 35.35
      变异系数CV/% 39.24 37.08 39.07 14.96 26.07 46.03 18.31

      表 2  长白落叶松人工林解析木和枝条分布统计表

      Table 2.  Statistics of sample trees and branch distribution for Larix olgensis planation

      项目
      Item
      胸径
      DBH/cm
      树高
      Tree height
      (HT)/cm
      冠长
      Crown length
      (CL)/m
      冠幅
      Crown width
      (CW)/m
      高径比
      HT/DBH
      每半米段的一级枝个数
      Number of first-order
      branch per 0.5 m
      建模数据(样本数量=40)
      Fitting data(sample size=40)
      最小值Minimum 2.00 3.80 2.30 0.55 0.66 1
      最大值Maximum 27.00 21.50 11.50 3.18 1.41 26
      平均值Mean 13.68 12.45 6.24 1.42 0.85 7.56
      标准差Std 5.42 4.45 2.38 0.56 0.26 4.21
      变异系数CV/% 39.61 35.74 38.14 39.44 30.59 55.69
      检验数据(样本数量=9)
      Validation data(sample size=9)
      最小值Minimum 3.80 5.30 2.52 0.51 0.70 1
      最大值Maximum 25.20 21.30 10.51 3.20 1.39 23
      平均值Mean 14.03 12.86 6.01 1.49 0.81 7.83
      标准差Std 5.88 4.13 2.49 0.59 0.22 3.86
      变异系数CV/% 37.56 32.12 41.43 39.60 27.16 49.30

      表 3  长白落叶松一级枝条数量模型所需变量及相关描述

      Table 3.  Symbol and associated description for variables tested in the first-order branch model for Larix olgensis planation

      变量符号
      Variable symbol
      变量描述
      Description
      DINC 冠层深度Distance from tree apex
      RDINC 相对冠层深度Relative distance from tree apex
      LnDINC 冠层深度的对数值Logarithm of distance from tree apex
      LnRDINC 相对冠层深度的对数值Logarithm of relative distance from tree apex
      RDINC2 相对冠层深度的平方Relative distance from tree apex squared
      DBH 树木胸高处的直径Diameter at breast height/cm
      HT 树高Tree height/m
      HT/DBH 树高和胸径的比值Ratio of tree height to diameter at breast height
      CL 树冠冠长Crown length/m
    • Poisson回归模型是以Poisson分布为基础,基于事件的计数变量而发展起来的回归模型,广泛应用于计数型数据的研究[20-24]。本研究的主要目标是建立长白落叶松一级枝条数量的混合效应模型。虽然负二项模型也可以用于计数型数据的研究,但本研究若以负二项模型为基础模型,有些重要的解释变量并没有进入到模型中,并且在引入不同的随机效应参数后,模型的拟合效果并没有明显提高,因此负二项模型并不适用于本研究,所以本文采用Poisson回归模型对长白落叶松的一级枝条数量进行模拟。其概率密度方程为:

      $$ P\left( {{{ \mathit{\boldsymbol{Y}}}_{i}}} \right) = \frac{{{{\rm{e}}^{ - {{ \mathit{\boldsymbol{\mu}}}_{i}}}}{{ \mathit{\boldsymbol{\mu}}}_{i}}^{{{ \mathit{\boldsymbol{Y}}}_{i}}}}}{{{{ \mathit{\boldsymbol{Y}}}_{i}}!}} $$ (1)

      式中:P(Yi)是单位面积或单位时间内事件Yi发生次数为yi的概率,在本文中指落叶松树干的某一位置处一级枝条数量为yi的概率;μi是Poisson分布的期望,满足期望与方差相等的条件,即E(Yi)=VAR(Yi)=μi;解释变量Xiμi通过链接函数g(μi)=ln(Xiβ)连接,其中Xi是包含DINC、LnRDINC、RDINC2、HT/DBH、DBH、CL等变量的向量。

      首先在传统模型的基础上引入了树木效应,发现模型的拟合效果很好。由于枝条数量可能会受到林分密度的影响,所以也应该考虑样地效应是否对枝条的数量存在影响。当继续引入样地效应时,发现模型的拟合结果并不显著,因此本研究不考虑样地效应。最终拟建立以Poisson分布为基础的单水平广义线性混合效应模型的表达形式为:

      $$ \left\{ {\begin{array} \ln\ln {{ \mathit{\boldsymbol{Y}}}_{i}}={{ \mathit{\boldsymbol{X}}}_{i}} \mathit{\boldsymbol{\beta}}+{{\mathit{\boldsymbol{Z}}}_{i}}{{\mathit{\boldsymbol{b}}}_{i}}+{{\mathit{\boldsymbol{\varepsilon}}}_{i}} \\ {{\mathit{\pmb{b}}}_{i}}\sim N(0, \mathit{\pmb{D}}) \\ \end{array}} \right. $$ (2)

      式中:Yini×1维第i株落叶松一级枝条个数观测值,XiZi分别为ni×p维固定效应设计矩阵和ni×q维随机效应设计矩阵;βbi分别为p×1维固定效应参数向量和q×1维随机效应参数向量;εini×1维误差向量;Dq×q维随机效应协方差矩阵。

      本研究采用极大似然估计法(maximum likelihood)进行基于Poisson分布的广义线性混合效应模型的参数估计。

    • 随机效应方差-协方差结构可以反映树木间是否存在差异。在确认混合模型方差-协方差结构之前,首先要确定参数是固定效应参数还是混合效应参数。在未知协方差结构时,可先将模型内所有参数都当作混合效应参数进行不同的组合,再对收敛的模型进行精度比较,选择收敛较好、精度较高、不存在过量参数化的模型作为最优模型[25]。林业研究中通常选用的随机效应参数方差-协方差结构为无结构矩阵(UN)、复合对称矩阵(CS)、对角矩阵(UN(1))等。

    • 本研究采用AIC、BIC、-2log likelihood这3个效果评价指标对相同参数个数的模型进行比较。对参数个数不同的模型进行LRT检验,从而选出最优的混合模型。采用确定系数(R2)、均方根误差(RMSE)、平均绝对误差(MAE)这3个检验指标对模型的模拟精度进行评价。

      确定系数(R2):

      $$ {R^2} = 1 - \frac{{\sum\limits_{i = 1}^n {{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}} }}{{\sum\limits_{i = 1}^n {{{\left( {{y_i} - {{\bar y}_i}} \right)}^2}} }} $$ (3)

      均方根误差(RMSE):

      $$ {\rm{RMSE}} = \sqrt {\frac{{\sum\limits_{i = 1}^n {{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}} }}{{n - p}}} $$ (4)

      平均绝对误差(MAE):

      $$ {\rm{MAE}} = \frac{{\sum\limits_{i = 1}^n {\left| {{y_i} - {{\hat y}_i}} \right|} }}{n} $$ (5)

      式中:yi为实测值,${{\hat y}_i} $为预测值,$ {{\bar y}_i} = \frac{{\sum\limits_{i = 1}^n {{y_i}} }}{n}$,n为样本个数,p为参数个数。

    • 运用表 2中建模时未使用的独立样本数据(检验数据)对模型进行独立性检验,对所选模型的预测性能进行评价。混合效应模型中固定效应部分的检验与传统模型的检验方法一致,但是随机效应部分的检验需要重新计算随机效应参数$\hat{b}_{i} $,进而算出检验样木的预测值,具体过程运用Matlab软件进行实现。然后,再用上述3个检验指标对验证数据的模拟精度进行评价,R2越大,RMSE和MAE越小, 说明模型的模拟精度越高。

      $$ \hat{\boldsymbol{b}}_{i}=\hat{\boldsymbol{D}} \hat{\boldsymbol{Z}}_{i}^{\mathrm{T}}\left(\hat{\boldsymbol{Z}}_{i} \hat{\boldsymbol{D}} \hat{\boldsymbol{Z}}_{i}^{\mathrm{T}}+\hat{\boldsymbol{R}}_{i}\right)^{-1} \hat{\boldsymbol{e}}_{i} $$ (6)

      式中:$ \hat{\boldsymbol{D}}$为随机效应参数的方差协方差矩阵;$\hat{\boldsymbol{R}}_{i}$为样木内方差协方差结构;$\hat{\boldsymbol{Z}}_{i} $为设计矩阵;$\hat{\boldsymbol{Z}}_{i}^{\mathrm{T}}$为设计矩阵的转置;$ \hat{\boldsymbol{e}}_{i}$为样木的实际值与用固定效应参数计算的预测值的差值。

    • 利用SAS 9.3软件中的PROC GLIMMIX模块,建立基于Poisson分布的长白落叶松每半米段一级枝条数量基础模型。通过比较R2、RMSE、MAE的大小以及模型解释变量获取的难易程度,发现含有6个变量(DINC、LnRDINC、RDINC2、HT/DBH、DBH、CL)的模型的拟合效果最好,而且模型中所有参数的估计值均极显著(P<0.001),因此选择含有这6个变量的模型作为最优基础模型。模型的表达式为:

      $$ \begin{array}{l} \ln Y = {\beta _o} + {\beta _1} \cdot {\rm{DINC}} + {\beta _2}\ln {\rm{RDINC}} + {\beta _3} \cdot {\rm{RDIN}}{{\rm{C}}^2} + \\ {\beta _4} \cdot {\rm{HT}}/{\rm{DBH}} + {\beta _5} \cdot {\rm{DBH}} + {\beta _6} \cdot {\rm{CL}} \end{array} $$ (7)
    • 通常情况下,参数效应和随机效应方差-协方差结构的确定是同时进行的。首先将随机效应的方差-协方差结构设为无结构矩阵(UN),利用SAS 9.3中的PROC GLIMMIX模块对不同随机参数的组合进行拟合,通过比较AIC、BIC、-2log likelihood等统计量指标来确定模型的拟合效果。3个指标的值越小,证明模型的拟合效果越好。模拟结果见表 4,限于篇幅,只给出了收敛的随机参数组合模型。

      表 4  基于树木效应的混合模型拟合结果

      Table 4.  Fitting results of mixed model based on individual tree effects

      参数个数
      Parameter
      number
      模型
      Model
      随机参数Random parameter AIC BIC -2log
      likelihood
      截距Intercept HT/DBH DINC LnRDINC CL RDINC2 DBH
      7 (7) 无None 2 610.85 2 640.17 2 596.85
      (8) 2 600.01 2 617.92 2 580.11
      (9) 2 599.16 2 619.67 2 582.16
      (10) 2 584.90 2 598.41 2 568.90
      8 (11) 2 600.07 2 618.98 2 581.47
      (12) 2 600.19 2 616.50 2 582.99
      (13) 2 580.35 2 593.86 2 564.35
      (14) 2 599.50 2 616.01 2 586.50
      10 (15) 2 575.01 2 591.90 2 555.01
      (16) 2 571.42 2 588.31 2 551.42
      (17) 2 574.70 2 591.98 2 554.70
      (18) 2 575.65 2 590.85 2 557.65
      (19) 2 576.51 2 593.40 2 556.61
      (20) 2 576.65 2 593.54 2 556.65
      (21) 2 572.97 2 589.86 2 552.97
      (22) 2 574.16 2 591.05 2 554.16
      (23) 2 568.91 2 585.80 2 548.91
      (24) 2 573.76 2 590.65 2 553.76
      (25) 2 573.08 2 589.97 2 553.08
      (26) 2 574.18 2 591.06 2 554.18
      (27) 2 583.32 2 600.20 2 563.32
      (28) 2 574.62 2 591.51 2 554.62
      (29) 2 580.54 2 592.43 2 561.54
      (30) 2 572.95 2 589.84 2 552.95
      13 25(31) 2 568.95 2 590.91 2 542.95
      26(32) 2 566.07 2 588.03 2 540.07
      27(33) 2 570.77 2 592.73 2 544.77
      28(34) 2 568.37 2 590.32 2 542.37
      29(35) 2 565.65 2 587.61 2 539.65
      30(36) 2 572.20 2 594.15 2 546.20
      31(37) 2 569.04 2 590.99 2 543.04
      32(38) 2 567.45 2 589.40 2 541.45
      33(39) 2 565.74 2 587.70 2 539.74
      34(40) 2 568.37 2 590.32 2 542.37
      注:参数个数包括固定效应参数个数、随机效应方差-协方差参数个数;带★变量表示包含此随机参数;表中的模型(7)即正文中的公式(7),即不含随机效应参数的最优基础模型。Notes: the number of parameters includes the fixed parameters, parameters in covariance-structure G for random effect; variables with ★ are included in the random parameters; model (7) in the table is the formula (7) model, which is the optimal basic model without random effect parameters.

      表 4得知,共有34个不同随机效应参数组合的模型可以收敛。对于参数个数相同的模型,通过比较AIC、BIC、-2log likelihood来分析模型的拟合效果。这3个指标的值越小,说明模型的拟合精度越高。根据表 4可以看出基础模型(7)的拟合效果不及任何一种混合模型的模拟效果,因此混合效应模型确实可以提高模型的预测精度。当模型只包含一个随机效应参数时,模型(13)的拟合效果最佳,即RDINC2作为随机效应参数,AIC=2 580.35,BIC=2 593.86,-2log likelihood=2 564.35;当模型中包含2个随机效应参数时,模型(23)的拟合效果最佳,即DINC、RDINC2作为随机效应参数,AIC=2 568.91,BIC=2 580.80,-2log likelihood=2 548.91;当模型中包含3个随机效应参数时,模型(35)的拟合效果最佳,即DINC、LnRDINC、RDINC2作为随机效应参数,AIC=2 565.65,BIC=2 587.61,-2log likelihood=2 539.65。通过比较发现模型(35)的AIC等值最小,有较好的拟合效果。为避免过多参数化,还要对不同参数个数的模型进行似然比检验,结果详见表 5

      表 5  不同参数个数的模型固定参数和方差组成的估计值

      Table 5.  Fixed parameters and variance component estimates of models with different numbers of parameters

      固定参数
      Fixed parameter
      模型(7)
      Model (7)
      模型(13)
      Model (13)
      模型(23)
      Model (23)
      模型(35)
      Model (35)
      截距Int 3.058 2 2.941 7 3.178 4 3.086 6
      DINC -0.002 1 -0.002 1 -0.002 2 -0.002 5
      LnRDINC 0.342 3 0.361 7 0.370 9 0.415 3
      RDINC2 -1.044 4 -1.153 0 -1.129 3 -1.046 1
      HT/DBH -0.282 2 -0.147 1 -0.280 3 -0.245 4
      DBH -0.015 3 -0.004 34 -0.007 4 -0.011 2
      CL 0.094 9 0.080 4 0.076 7 0.105 0
      随机效应方差-协方差结构
      Variance of random effect-covariance
      structure
      [0.279 1] $\left[ {\begin{array}{*{20}{c}} {0.000\;01}&{ - 0.000\;99}\\ { - 0.000\;99}&{1.446\;3} \end{array}} \right]$ $ \left[ {\begin{array}{*{20}{c}} {0.000\;001}&{0.000\;033}&{ - 0.001\;16}\\ {0.000\;033}&{0.004\;852}&{ - 0.011\;22}\\ { - 0.001\;16}&{ - 0.011\;22}&{1.565\;0} \end{array}} \right]$
      LRT 32.5 15.44 9.26
      P <0.001 <0.001 0.026
      Pearson χ2/自由度 1.73 1.50 1.38 1.20
      Pearson χ2/degree

      对4个模型进行两两比较,结果显示P值均小于0.05,证明似然比检验结果均显著,也就是说模型(7)与(13)显著不同,模型(13)与(23)显著不同,模型(23)与(35)显著不同,因此选择模型(35)作为最优混合模型。模型(7)的Pearson χ2/自由度=1.73,说明该模型存在过度离散的问题。然而,随着随机效应参数的增多,Pearson χ2/自由度的值相应减小,由1.73减小到1.20,接近1,说明模型(35)几乎不存在过度离散的问题。

    • 当考虑树木效应时,随机效应参数反映了树木间的差异,随着随机参数的增多,模型的模拟效果逐步提高。选择其他2种随机效应方差-协方差结构,即复合对称矩阵(CS)和对角矩阵(UN(1))对模型进行测试,结果见表 6

      表 6  基于随机参数不同方差-协方差结构的混合模型模拟结果比较

      Table 6.  Comparisons of mixed model based on different variance-covariance structure of random parameters

      方差-协方差结构
      Variance-covariance structure
      参数个数
      Parameter number
      AIC BIC -2log likelihood LRT P
      复合对称矩阵Composite symmetric matrix (CS) 9 2 586.71 2 601.91 2 579.71
      对角矩阵Diagonal matrix (UN(1)) 13 2 580.99 2 597.87 2 560.99 18.72 < 0.001
      无结构矩阵Non-structural matrix (UN) 13 2 565.65 2 587.61 2 539.65 21.34

      表 6得知3种随机效应方差-协方差结构均可以提高混合模型的拟合效果。对3个模型进行LRT检验,结果显示P值小于0.001,说明以CS为矩阵的模型与以UN(1)为矩阵的模型显著不同,并且后者更好。以UN为矩阵的模型的AIC等值均小于以UN(1)为矩阵的模型,所以最终选择无结构矩阵(UN)作为该模型的随机效应方差-协方差结构。长白落叶松每半米段一级枝条数量最优模型(35)的固定效应部分形式如下:

      $$ \begin{array}{*{20}{c}} {\ln Y = 3.0866 - 0.0025{\rm{DINC}} + 0.4153\ln {\rm{RDINC}} - }\\ {1.0461{\rm{RDIN}}{{\rm{C}}^2} - 0.2454{\rm{HT}}/{\rm{DBH}} - 0.0112{\rm{DBH}} + }\\ {0.105{\rm{CL}}} \end{array} $$

      可以看出:CL的参数估计值为正值,说明随着冠长的增大,每半米段一级枝条数量相应增多;HT/DBH、DBH、DINC的参数估计值为负值,说明随着高径比、胸径、冠层深度的增大,每半米段一级枝条数量相应减少;LnRDINC的参数估计值为正值,RDINC2的参数估计值为负值,说明随着相对冠层深度的增大,每半米段一级枝条数量在树冠范围内存在最大值。

      本文利用最优模型的固定效应参数模拟了相对冠层深度与每半米段一级枝条数量分布的关系(图 1)。图 1a表明:随着相对冠层深度的变化,不同单木内的枝条数量变化规律一致;从梢头开始,随着相对冠层深度的增加,每半米段内的枝条数量先逐渐增大,大约在距离梢头20%的位置(即相对冠层深度为0.2处)达最大值;之后随着相对冠层深度的继续增加,每半米段内的一级枝条数量逐渐减小,并在相对冠层深度为1时与横轴交于非零处。对比不同高径比的3棵林木,在冠内的同一位置处,高径比越大的树木,每半米段的枝条数量越少(图 1a)。图 1b表明:随着相对冠层深度的变化,不同单木内的枝条数量变化规律也一致;随着相对冠层深度的增加,每半米段内的枝条数量先增大后减小,在相对冠层深度为0.2左右达到峰值。对比不同冠长的3棵林木,在树冠的中上部,冠长越大的树木,每半米段内的枝条数量越多;在树冠的中下部,冠长越大的树木,每半米段内的枝条数量越少(图 1b),即冠长较大的树木在树冠的中上部有更多的枝,而冠长较小的树木在树冠的中下部有更多的枝。

      图  1  人工长白落叶松每半米段一级枝条数量分布

      Figure 1.  Predicted value of number of first-order branches per 0.5 m

      分别绘制每半米段一级枝条数量分布基础模型即模型(7)和最优模型即模型(35)的实测值与预测值线性相关图(图 2a2b)。根据每棵长白落叶松一级枝条实测数据,绘制全树一级枝条数量分布基础模型即模型(7)和最优模型即模型(35)的实测值与预测值线性相关图(图 2c2d)。图 2a2b的结果表明每半米段一级枝条数量的实测值与预测值大致接近,大部分数据点分布在拟合直线两侧,但图 2b中的点更接近拟合直线,而且R2也由0.504 2增加到了0.660 6,说明混合模型确实可以提高模型的拟合效果。图 2d中的全树一级枝条数量的数据点较图 2c更向拟合直线靠拢,而且R2由0.823 1增加到了0.965 2,进一步证明了加入树木水平的混合效应模型可以得到更好的预测结果。

      图  2  人工长白落叶松一级枝条数量的实测值与预测值散点图

      Figure 2.  Observed values and predicted values of first-order branch for Larix olgensis planation

      绘制最优混合效应模型预测的人工长白落叶松每半米段一级枝条数量与相对冠层深度的关系图(图 3),可以看出:含有5~8个枝的米段在树冠上的分布范围最广,说明含有5~8个枝的米段最多。含有枝条数量较多的米段主要分布在相对冠层深度为0.2的位置,说明枝条数量在距离梢头20%的位置达到最大值。从图 3可以看出枝条数量主要集中在树冠的中上部,树冠下部的枝条数量较少。

      图  3  人工长白落叶松相对冠层深度与每半米段的预测枝条数量关系

      Figure 3.  Relationship between relative distance from tree apex and the number of predicted branches per 0.5 m section for Larix olgensis planation

    • 采用独立样本数据对模型进行独立性检验;运用Matlab软件计算出检验数据的随机效应参数值,结果见表 7

      表 7  检验数据的随机参数估计结果

      Table 7.  Random parameter estimation result of validation data

      检验数据
      Validation data
      随机参数
      Random parameter
      DINC LnRDINC RDINC2
      b1 -0.000 5 -0.013 2 0.383 1
      b2 0.000 6 -0.022 6 -0.971 2
      b3 0.000 9 0.014 7 -1.068 5
      b4 0.000 8 -0.073 4 -1.525 0
      b5 0.000 1 -0.109 2 -1.883 0
      b6 0.000 7 0.044 0 -0.629 0
      b7 -0.001 6 -0.114 4 1.453 0
      b8 -0.001 6 -0.173 4 1.109 1
      b9 0.001 4 0.120 6 -1.083 0

      根据计算出的随机效应参数$\hat{\boldsymbol{b}}_{i}$,计算出检验数据的预测值,从而对基础模型和混合模型的预测精度进行比较(表 8),结果表明:混合模型的确定系数R2较基础模型有所提高,均方根误差RMSE和平均绝对误差MAE均比基础模型的低,说明引入随机参数的混合模型可以提高预测精度。

      表 8  基础模型和混合效应模型模拟结果比较

      Table 8.  Comparison of based model and mixed-effect model

      评价指标
      Evaluation index
      固定模型
      Fixed model
      混合模型
      Mixed model
      确定系数
      Determination coefficient (R2)
      0.510 0.669
      均方根误差
      Root mean square error (RMSE)
      3.396 3.012
      平均绝对误差
      Average absolute error (MAE)
      2.601 2.250
    • 本研究以黑龙江省佳木斯市孟家岗林场49棵长白落叶松每半米段的一级枝条数量为研究对象,在确定每半米段一级枝条数量分布的最优基础模型后,加入树木效应,建立以Poisson分布为基础的每半米段一级枝条数量分布的混合效应模型。研究结果表明含有DINC、LnRDINC、RDINC2、HT/DBH、DBH、CL这6个变量的模型拟合效果最好,因此确定该模型为每半米段一级枝条数量分布的最优基础模型。考虑树木效应,在最优基础模型中引入不同的随机效应参数组合,发现任意随机效应参数组合获得的混合模型的拟合效果均比传统模型要好,并且随着随机效应参数的增多,模型的拟合效果逐步提高。对于参数个数相同的模型,通过比较AIC、BIC、-2log likelihood等统计量指标选取最优模型;对于参数个数不同的模型,根据LRT检验和过度离散指标Pearson χ2/自由度的结果选取最优模型。本文还选择了不同的随机效应方差-协方差结构对模型进行测试,最终将基于Poisson分布的含有DINC、LnRDINC、RDINC2这3个随机效应参数的无结构矩阵模型作为每半米段一级枝条数量分布的最优混合效应模型。

      根据获得的最优混合效应模型的固定效应参数得知,DINC的参数估计值为负值,说明冠层深度与每半米段一级枝条数量呈负相关。随着冠层深度的增大,枝条年龄逐渐增大,光合作用能力和气孔传导能力会逐渐减弱,从而造成枝条数量的损失,因此冠层深度可作为一个可以很好地解释枝条数量损失原因的变量[26]。HT/DBH、DBH的参数估计值为负值,CL的参数估计值为正值,说明高径比、胸径与每半米段一级枝条数量呈负相关,冠长与每半米段一级枝条数量呈正相关,也就是说高径比和胸径越小、冠长越大的树木,每半米段一级枝数量越多,这一结果与前人的研究结果一致[15, 27-29]。高径比与林分密度关系密切,是反映树木活力的重要指标,研究表明枝条数量会随树木活力的增强而增多[9]。相比于冠比,冠长可以更好地反映枝条数量的分布规律[15]。LnRDINC的参数估计值为正值,RDINC2的参数估计值为负值,说明在一定范围内随着相对冠层深度的增大,每半米段一级枝条数量存在峰值。从绘制的相对冠层深度与每半米段一级枝条数量分布关系图来看,起初枝条数量随相对冠层深度的增加而增多,大约在距离梢头20%的位置(即相对冠层深度为0.2处)枝条数量达到最大值;之后随着相对冠层深度的继续增加,每半米段内的一级枝条数量逐渐减小,这与他人的研究结果一致[27]。起初每半米段一级枝条数量增多是由于枝条在不断生长,后来每半米段一级枝条数量减少是由于随着相对冠层深度的增加,树冠上部枝条和相邻木的竞争逐渐激烈,导致下部枝条被遮挡因而造成了数量的损失。Thorpe等[30]同样发现枝条的数量会受到当地竞争条件的影响。

      Poisson回归模型是一种常用的计数模型,本文以Poisson回归模型为基础,运用SAS 9.3中的PROC GLIMMIX模块对每半米段一级枝条数量进行拟合,运用极大似然估计法进行参数估计,得到了长白落叶松每半米段一级枝条数量分布的预测模型,可为今后的林分生长动态研究提供理论依据。由于枝条数量受立地条件、遗传作用等很多因素的影响,随着数据的积累,今后可开展不同林分条件下枝条数量分布的研究。

参考文献 (30)

目录

    /

    返回文章
    返回