Tree biomass estimation using allometry
Tree- and Stand- Level Biomass Estimation for Scots Pine with Allometric Equations, Biomass Conversion and Expansion Factors
Viet Nguyen
Faculty of Forest and Environment Eberswalde University for Sustainable Development
Duc.Nguyen@hnee.de
June, 2021
🛈 Information: This project was hosted at GitHub.
1. Introduction
Biomass is increasingly being viewed as a future renewable energy source and is being used in biofuel production, among other things. Biomass research is a useful source of information on the total amount of carbon and nutrients in a forest ecosystem, and it is utilized in natural carbon cycle modeling (Socha J, Wężyk P, 2004). Furthermore, many ecologists are focusing their research on carbon cycling in the biosphere in the age of climate change. Forests are one of the most important carbon repositories among terrestrial ecosystems, collecting carbon globally. Thus, their contribution to mitigating the impact of climate change is critical. Adaptation strategies must be integrated into forest management since climate change will cause changes in tree species’ physiology and ranges. There is a pressing need for reliable carbon assessment tools to properly plan for enhancing forest carbon sequestration and calculating their ability to counteract climate change (A. Jagodziński et al., 2018).
Carbon reserves in forest ecosystems must be precisely estimated as a foundation for future actions. Most inventories concentrate on biomass estimation since the carbon content of plant tissues is essentially stable. Thus, biomass assessment accuracy is critical for carbon reporting. Different countries employ different methodologies to country-scale carbon assessment, with inventory levels (tree or stand level) and accuracy varying. Carbon assessment can be based on generalized coefficients, according to the Intergovernmental Panel on Climate Change (IPCC) guidelines. Therefore, carbon reservoirs in forest stands may be underestimated or overestimated. For that reason, IPCC recommends the usage of locally-fitted models of biomass estimation wherever it is possible (A. Jagodziński et al., 2018).
Various tree biomass estimation methods, both direct and indirect, are applied on forest lands. There are two main indirect approaches for biomass estimation: the use of allometric equations (AEs) or biomass conversion and expansion factors (BCEFs). Because AEs are applied at the tree level, while BCEFs are done at the tree-stand level, these two methods differ in accuracy and labor procedures. AEs are tree-level computations that are based on relationships between tree biomass and commonly used measurements, such as diameter at breast height (DBH) or total height (H) (A. M. Jagodziński et al., 2018)c. AEs have the disadvantage of requiring a large amount of labor and accompanying costs, despite their relatively high accuracy. Besides, the empirical equation is considered site-specific and thus may not accurately reflect tree biomass under other site conditions. BCEFs, on the other hand, are literally proportions of forest stand biomass and tree stand volume, allow calculations at the tree-stand level. Expansion factor (BEF) allows converting merchantable timber volume to aboveground, belowground, or total biomass, while merchantable timber volume convergence into stem biomass could be based on conversion factors (BCF). Furthermore, because forest inventories usually lack tree-level data, BCEFs are widely used for large-scale calculations (Neumann et al., 2016). However, this approach may cause double bias by errors from statistical models and errors from forest inventory. In general, two different approaches are used to convert easily measurable tree or stand attributes to total stand biomass.
When modeling biomass equation for total above/belowground biomass of trees and their components, it is worth considering the logical assumption that the sum of the biomass components of the tree estimated using equations should be equal to the estimated biomass of the whole tree (so-called additivity of the biomass models) (Bronisz et al., 2021). This assumption can be met by applying seemingly unrelated regression (SUR). The total biomass regression model is defined as the sum of the individually formulated best regression models of the biomass of its components, with cross-modal correlation used to represent the combined behavior of the separate components. Hence, SUR is increasingly being used in the development of various additive biomass models.
This study aimed to develop a comprehensive set of allometric equations as well as biomass conversion and expansion factors for Scots Pine above-ground biomass at tree and stand levels. The empirical equations and BCEFs will be developed for determining the dry mass of the aboveground part of Scots pine trees and its components (stem, branches, and leaves), using the SUR model with DBH and H as predictors for empirical equations and age, volume, the quadratic mean diameter at the breast height, mean height and basal area as predictors for BCEFs. We hypothesized that (1) Seemingly Unrelated regressions of biomass prediction offer a wide range of capability for application, (2) forest stand features: height, volume, age, diameter, and basal area will influence forest stand biomass and BCEFs.
2. Material and Method
2.1. Study site and species
Scots pine (Pinus sylvestris L.) is a pioneering coniferous tree species that covers more than 28 million hectares in Europe and accounts for more than 20% of total timber productivity. Scots pine may thrive in a variety of environments, from dry, poor arenosols to wet, rich alluvial soils, although it can only breed naturally on poor to medium fertile podzols and brunic arenosols (Ellenberg, 1988). In Poland, Scots pine occupies 5.4 million ha (58.1 percent of forest area) and has a merchantable volume of 1519 million m3 (60.8 percent of total wood resources) according to the National Forest Inventory from 2011 to 2015 (Forest Data Bank, 14-Jun-21). Therefore, Scots pine stands are important carbon reservoirs, contributing significantly to global carbon storage through the biomass of needles, branches, stems, and roots. Due to drought-mediated limitation in southern Europe and cold-mediated limitation in the north, Scots pine reaches the highest dimensions and biomass in Central Europe (Oleksyn et al., 1999).
The research area was established 18 plots in Scots pine stands ranging around 20 - 100 years old (Appendix II). The site types including dry coniferous forest, fresh coniferous forest, fresh mixed coniferous forest according to Polish site classification (Appendix I), which are facilitated condition for Scot pine growth.
#### overview of the study forest stand ####
# load library for the project
library(systemfit)
library(ggplot2)
library(pastecs)
library(gridExtra)
library(readxl)
library(lmtest)
library(dplyr)
library(tidyr)
library(grid)
library(minpack.lm)
#library(ggpmisc)
# import data
data_tree <- read.csv("data/Scots_pine_2_trees.csv", sep=";")
data_plot <- read.csv("data/Scots_pine_2_plots.csv", sep=";")
# take a look at the data
data_treedata_plot# rename tree-level data column
data_tree$D<-data_tree$DBH
data_tree$D2<-data_tree$D^2
data_tree$H<-data_tree$H_m
data_tree$H2<-data_tree$H^2
data_tree$D2H<-data_tree$D2*data_tree$H
data_tree$SDB<-data_tree$SDB_kg
data_tree$BDB<-data_tree$BDB_kg
data_tree$LDB<-data_tree$LDB_kg
# subset data
tree<-with(data_tree, data.frame(Plot,Tree,Site,Age_class,H,D,D2,H2,D2H,SDB,BDB,LDB))
attach(tree)
# data description
stat.desc(data_plot)summary(data_plot)## Plot Site Tree.species Age
## Min. : 1.00 Length:18 Length:18 Min. : 23.00
## 1st Qu.: 5.25 Class :character Class :character 1st Qu.: 35.25
## Median : 9.50 Mode :character Mode :character Median : 56.50
## Mean : 9.50 Mean : 60.72
## 3rd Qu.:13.75 3rd Qu.: 84.50
## Max. :18.00 Max. :106.00
## BA_m2_ha N_trees_ha Dg_mm Hg_m
## Min. :17.43 Min. : 535.0 Min. : 93.0 Min. : 8.49
## 1st Qu.:26.20 1st Qu.: 866.8 1st Qu.:113.0 1st Qu.:11.94
## Median :28.96 Median :1238.0 Median :167.0 Median :14.01
## Mean :29.31 Mean :1923.6 Mean :164.4 Mean :15.11
## 3rd Qu.:31.49 3rd Qu.:2729.2 3rd Qu.:205.8 3rd Qu.:17.86
## Max. :38.18 Max. :5109.0 Max. :272.0 Max. :22.99
## V_m3_ha SDB_Mg_ha BDB_Mg_ha LDB_Mg_ha
## Min. : 83.84 Min. : 45.94 Min. : 4.773 Min. : 5.638
## 1st Qu.:178.98 1st Qu.: 90.66 1st Qu.: 9.017 1st Qu.: 6.516
## Median :219.19 Median :100.73 Median : 9.917 Median : 7.139
## Mean :219.16 Mean :104.95 Mean :10.415 Mean : 7.832
## 3rd Qu.:247.41 3rd Qu.:121.13 3rd Qu.:12.175 3rd Qu.: 9.233
## Max. :370.63 Max. :177.47 Max. :17.279 Max. :13.238
## ADB_Mg_ha
## Min. : 56.57
## 1st Qu.:106.66
## Median :117.33
## Mean :123.19
## 3rd Qu.:143.91
## Max. :204.19stat.desc(data_tree)summary(data_tree)## Plot Tree Site Age_class
## Min. : 1.0 Min. :1 Length:90 Length:90
## 1st Qu.: 5.0 1st Qu.:2 Class :character Class :character
## Median : 9.5 Median :3 Mode :character Mode :character
## Mean : 9.5 Mean :3
## 3rd Qu.:14.0 3rd Qu.:4
## Max. :18.0 Max. :5
## DBH_mm H_m SDB_kg BDB_kg
## Min. : 67.0 Min. : 6.95 Min. : 8.288 Min. : 0.6685
## 1st Qu.:107.0 1st Qu.:12.66 1st Qu.: 30.671 1st Qu.: 2.8764
## Median :149.0 Median :15.03 Median : 69.549 Median : 5.5347
## Mean :157.4 Mean :15.51 Mean :100.961 Mean : 8.9351
## 3rd Qu.:194.0 3rd Qu.:18.30 3rd Qu.:134.247 3rd Qu.:11.3319
## Max. :315.0 Max. :25.65 Max. :452.091 Max. :55.4057
## LDB_kg ADB_kg D D2
## Min. : 0.5174 Min. : 9.702 Min. : 67.0 Min. : 4489
## 1st Qu.: 2.7122 1st Qu.: 37.457 1st Qu.:107.0 1st Qu.:11449
## Median : 4.5126 Median : 79.897 Median :149.0 Median :22202
## Mean : 5.8221 Mean :115.718 Mean :157.4 Mean :28372
## 3rd Qu.: 7.2788 3rd Qu.:154.302 3rd Qu.:194.0 3rd Qu.:37636
## Max. :27.7573 Max. :500.290 Max. :315.0 Max. :99225
## H H2 D2H SDB
## Min. : 6.95 Min. : 48.3 Min. : 36184 Min. : 8.288
## 1st Qu.:12.66 1st Qu.:160.3 1st Qu.: 140954 1st Qu.: 30.671
## Median :15.03 Median :225.8 Median : 337840 Median : 69.549
## Mean :15.51 Mean :255.9 Mean : 514397 Mean :100.961
## 3rd Qu.:18.30 3rd Qu.:334.9 3rd Qu.: 709631 3rd Qu.:134.247
## Max. :25.65 Max. :657.9 Max. :2413152 Max. :452.091
## BDB LDB
## Min. : 0.6685 Min. : 0.5174
## 1st Qu.: 2.8764 1st Qu.: 2.7122
## Median : 5.5347 Median : 4.5126
## Mean : 8.9351 Mean : 5.8221
## 3rd Qu.:11.3319 3rd Qu.: 7.2788
## Max. :55.4057 Max. :27.7573
Table 1. Overview of the study forest stand
2.2. Material
Tree heights and diameters were measured. Besides, the fresh weights of the biomass components were selected and collected, then they were oven-dried to constant mass and were weighed again. Using proportions of dry and fresh masses of samples and total fresh masses of biomass components obtained on the field, we calculated the total dry mass of each biomass component including stem dry biomass (SDB), branches dry biomass (BDB), leaves dry biomass (LDB), and aboveground biomass (ADB).
# data visualization
p1 <- ggplot(data = tree, aes(D)) +
geom_histogram()+theme_classic()+labs(title = "(a)", x="DBH (mm)")
p2 <- ggplot(data = tree, aes(H)) +
geom_histogram()+theme_classic()+labs(title = "(b)", x="H (m)")
p3 <-ggplot(data = tree, aes(D,H)) +
geom_point()+theme_classic()+labs(title = "(c)", x="DBH (mm)", y="H (m)")
p4 <-ggplot(data = tree, aes(D,SDB)) +
geom_point()+theme_classic()+labs(title = "(D)", x="DBH (mm)", y="SDB (kg)")
p5 <-ggplot(data = tree, aes(D,BDB)) +
geom_point()+theme_classic()+labs(title = "(e)", x="DBH (mm)", y="BDB (kg)")
p6 <-ggplot(data = tree, aes(D,LDB)) +
geom_point()+theme_classic()+labs(title = "(f)", x="DBH (mm)", y="LDB (kg)")
p7 <-ggplot(data = tree, aes(H,SDB)) +
geom_point()+theme_classic()+labs(title = "(g)", x="H (m)", y="SDB (kg)")
p8 <-ggplot(data = tree, aes(H,BDB)) +
geom_point()+theme_classic()+labs(title = "(h)", x="H (m)", y="BDB (kg)")
p9 <-ggplot(data = tree, aes(H,LDB)) +
geom_point()+theme_classic()+labs(title = "(i)", x="H (m)", y="LDB (kg)")
tree_description <- grid.arrange(p1,p2,p3,p4,p5,p6,p7,p8,p9, ncol=3,
top= textGrob("Figure 1. Distribution of DBH (a), H (b), and their relationship with tree components (c-i)\n", gp=gpar(fontsize=16)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Table 2. Overview of tree characteristics 
2.3. Method
2.3.1. Empirical equaitons
Biomass components were estimated based on equations in previous studies, some of which consider exclusively on DBH or H while others based on both parameters.
\[ (1)\quad W = a + b×DBH \\ (2)\quad W = a + b×DBH^2 \\ (3)\quad W = a + b×H \\ (4)\quad W = a + b×DBH + c×H \\ (5)\quad W = a +b×DBH^2 + c×H^2 \\ (6)\quad W = a + b×DBH^2×H \\ (7)\quad W = a + b×DBH^2 + c×H + d×DBH^2×H \\ (8)\quad W = a×DBH^b \\ (9)\quad W = a×DBH^b×H^c \\ (10)\quad W = a×(DBH^2×H)^b \\ \]
Where: W- dry mass of the considered biomass component (kg), DBH - diameter at breast height (mm), H - height (m), a,b,c - model parameters.
All models were fit to data separately for the stem, branches, and foliage biomass and compared based on the following goodness-of-fit measures: coefficient of determination (R2), residual standard error (RSE), and The Akaike information criterion (AIC).
#### regression models for stem dry tree
Ms1 <- nls(SDB~b0+b1*D,start=list(b0=1,b1=1),data=tree)
Ms2 <- nls(SDB~b0+b1*D2,start=list(b0=1,b1=1),data=tree)
Ms3 <- nls(SDB~b0+b1*H,start=list(b0=1,b1=1),data=tree)
Ms4 <- nls(SDB~b0+b1*D+b2*H,start=list(b0=1,b1=1,b2=1),data=tree)
Ms5 <- nls(SDB~b0+b1*D2+b2*H2,start=list(b0=1,b1=1,b2=1),data=tree)
Ms6 <- nls(SDB~b0+b1*D2H,start=list(b0=1,b1=1),data=tree)
Ms7 <- nls(SDB~b0+b1*D2+b2*H+b3*D2H,start=list(b0=1,b1=1,b2=1,b3=1),data=tree)
Ms8 <- nls(SDB~b0*D^b1,start=list(b0=1,b1=1),control=list(maxiter=100),data=tree)
Ms9 <- nls(SDB~b0*D^b1*H^b2,start=list(b0=1,b1=2,b2=1),data=tree,algorithm="port",control=list(maxiter=100))
Ms10 <- nls(SDB~b0*D2H^b1,start=list(b0=1,b1=1),data=tree)
## calculate goodness-of-fit measurements
# R2
Ms1r2 <- 1-(sum({summary(Ms1)$residuals}^2)/sum({SDB-mean(SDB)}^2))
Ms2r2 <- 1-(sum({summary(Ms2)$residuals}^2)/sum({SDB-mean(SDB)}^2))
Ms3r2 <- 1-(sum({summary(Ms3)$residuals}^2)/sum({SDB-mean(SDB)}^2))
Ms4r2 <- 1-(sum({summary(Ms4)$residuals}^2)/sum({SDB-mean(SDB)}^2))
Ms5r2 <- 1-(sum({summary(Ms5)$residuals}^2)/sum({SDB-mean(SDB)}^2))
Ms6r2 <- 1-(sum({summary(Ms6)$residuals}^2)/sum({SDB-mean(SDB)}^2))
Ms7r2 <- 1-(sum({summary(Ms7)$residuals}^2)/sum({SDB-mean(SDB)}^2))
Ms8r2 <- 1-(sum({summary(Ms8)$residuals}^2)/sum({SDB-mean(SDB)}^2))
Ms9r2 <- 1-(sum({summary(Ms9)$residuals}^2)/sum({SDB-mean(SDB)}^2))
Ms10r2 <- 1-(sum({summary(Ms10)$residuals}^2)/sum({SDB-mean(SDB)}^2))
# RSE
Mssrrse <- c(summary(Ms1)$sigma,summary(Ms2)$sigma,summary(Ms3)$sigma,summary(Ms4)$sigma,summary(Ms5)$sigma,summary(Ms6)$sigma,summary(Ms7)$sigma,summary(Ms8)$sigma,summary(Ms9)$sigma,summary(Ms10)$sigma)
# AIC
Mssraic <- AIC(Ms1,Ms2,Ms3,Ms4,Ms5,Ms6,Ms7,Ms8,Ms9,Ms10)$AIC
# R2
Mssrr2 <- c(Ms1r2,Ms2r2,Ms3r2,Ms4r2,Ms5r2,Ms6r2,Ms7r2,Ms8r2,Ms9r2,Ms10r2)
# show goodness-of-fit measurements
Mssrrse## [1] 31.44427 19.42744 42.86856 30.64686 17.30065 12.81482 12.54961 17.42277
## [9] 12.17437 12.42762Mssraic## [1] 880.0654 793.3900 935.8514 876.4132 773.4917 718.4948 716.6615 773.7865
## [9] 710.2377 712.9723Mssrr2## [1] 0.8932732 0.9592599 0.8016335 0.8997697 0.9680587 0.9822738 0.9833863
## [8] 0.9672339 0.9841831 0.9833288#### regression models for branches dry tree - BDB [kg]
Mb1 <- nls(BDB~b0+b1*D,start=list(b0=1,b1=1),data=tree)
Mb2 <- nls(BDB~b0+b1*D2,start=list(b0=1,b1=1),data=tree)
Mb3 <- nls(BDB~b0+b1*H,start=list(b0=1,b1=1),data=tree)
Mb4 <- nls(BDB~b0+b1*D+b2*H,start=list(b0=1,b1=1,b2=1),data=tree)
Mb5 <- nls(BDB~b0+b1*D2+b2*H2,start=list(b0=1,b1=1,b2=1),data=tree)
Mb6 <- nls(BDB~b0+b1*D2H,start=list(b0=1,b1=1),data=tree)
Mb7 <- nls(BDB~b0+b1*D2+b2*H+b3*D2H,start=list(b0=1,b1=1,b2=1,b3=1),data=tree)
#Mb8 <- nls(BDB~b0*D^b1,start=list(b0=1,b1=1),control=list(maxiter=100),data=tree)
Mb9 <- nls(BDB~b0*D^b1*H^b2,start=list(b0=1,b1=2,b2=1),data=tree,algorithm="port",control=list(maxiter=100))
Mb10 <- nls(BDB~b0*D2H^b1,start=list(b0=1,b1=0),data=tree,control=list(maxiter=100))
## calculate goodness-of-fit measurements
# R2
Mb1r2 <- 1-(sum({summary(Mb1)$residuals}^2)/sum({BDB-mean(BDB)}^2))
Mb2r2 <- 1-(sum({summary(Mb2)$residuals}^2)/sum({BDB-mean(BDB)}^2))
Mb3r2 <- 1-(sum({summary(Mb3)$residuals}^2)/sum({BDB-mean(BDB)}^2))
Mb4r2 <- 1-(sum({summary(Mb4)$residuals}^2)/sum({BDB-mean(BDB)}^2))
Mb5r2 <- 1-(sum({summary(Mb5)$residuals}^2)/sum({BDB-mean(BDB)}^2))
Mb6r2 <- 1-(sum({summary(Mb6)$residuals}^2)/sum({BDB-mean(BDB)}^2))
Mb7r2 <- 1-(sum({summary(Mb7)$residuals}^2)/sum({BDB-mean(BDB)}^2))
#Mb8r2 <- 1-(sum({summary(Mb8)$residuals}^2)/sum({BDB-mean(BDB)}^2))
Mb9r2 <- 1-(sum({summary(Mb9)$residuals}^2)/sum({BDB-mean(BDB)}^2))
Mb10r2 <- 1-(sum({summary(Mb10)$residuals}^2)/sum({BDB-mean(BDB)}^2))
# RSE
Mbsrrse <- c(summary(Mb1)$sigma,summary(Mb2)$sigma,summary(Mb3)$sigma,summary(Mb4)$sigma,summary(Mb5)$sigma,summary(Mb6)$sigma,summary(Mb7)$sigma,summary(Mb9)$sigma,summary(Mb10)$sigma)
# AIC
Mbsraic <- AIC(Mb1,Mb2,Mb3,Mb4,Mb5,Mb6,Mb7,Mb9,Mb10)$AIC
#R2
Mbsrr2 <- c(Mb1r2,Mb2r2,Mb3r2,Mb4r2,Mb5r2,Mb6r2,Mb7r2,Mb9r2,Mb10r2)
# show goodness-of-fit measurements
Mbsrrse## [1] 5.336924 4.744390 6.965418 5.124863 4.436163 4.990908 4.472874 4.402668
## [9] 4.991713Mbsraic## [1] 560.8233 539.6397 608.7588 554.4965 528.5200 548.7576 530.9629 527.1557
## [9] 548.7866Mbsrr2## [1] 0.6955300 0.7593847 0.4813710 0.7224357 0.7920236 0.7337304 0.7909975
## [8] 0.7951524 0.7336445#### regression models for leaves dry tree - LDB [kg]
Ml1 <- nls(LDB~b0+b1*D,start=list(b0=1,b1=1),data=tree)
Ml2 <- nls(LDB~b0+b1*D2,start=list(b0=1,b1=1),data=tree)
Ml3 <- nls(LDB~b0+b1*H,start=list(b0=1,b1=1),data=tree)
Ml4 <- nls(LDB~b0+b1*D+b2*H,start=list(b0=1,b1=1,b2=1),data=tree)
Ml5 <- nls(LDB~b0+b1*D2+b2*H2,start=list(b0=1,b1=1,b2=1),data=tree)
Ml6 <- nls(LDB~b0+b1*D2H,start=list(b0=1,b1=1),data=tree)
Ml7 <- nls(LDB~b0+b1*D2+b2*H+b3*D2H,start=list(b0=1,b1=1,b2=1,b3=1),data=tree)
Ml8 <- nls(LDB~b0*D^b1,start=list(b0=1,b1=0),data=tree)
Ml9 <- nls(LDB~b0*D^b1*H^b2,start=list(b0=1,b1=2,b2=1),data=tree)
Ml10 <- nls(LDB~b0*D2H^b1,start=list(b0=1,b1=0),data=tree)
## calculate goodness-of-fit measurements
# R2
Ml1r2 <- 1-(sum({summary(Ml1)$residuals}^2)/sum({LDB-mean(LDB)}^2))
Ml2r2 <- 1-(sum({summary(Ml2)$residuals}^2)/sum({LDB-mean(LDB)}^2))
Ml3r2 <- 1-(sum({summary(Ml3)$residuals}^2)/sum({LDB-mean(LDB)}^2))
Ml4r2 <- 1-(sum({summary(Ml4)$residuals}^2)/sum({LDB-mean(LDB)}^2))
Ml5r2 <- 1-(sum({summary(Ml5)$residuals}^2)/sum({LDB-mean(LDB)}^2))
Ml6r2 <- 1-(sum({summary(Ml6)$residuals}^2)/sum({LDB-mean(LDB)}^2))
Ml7r2 <- 1-(sum({summary(Ml7)$residuals}^2)/sum({LDB-mean(LDB)}^2))
Ml8r2 <- 1-(sum({summary(Ml8)$residuals}^2)/sum({LDB-mean(LDB)}^2))
Ml9r2 <- 1-(sum({summary(Ml9)$residuals}^2)/sum({LDB-mean(LDB)}^2))
Ml10r2 <- 1-(sum({summary(Ml10)$residuals}^2)/sum({LDB-mean(LDB)}^2))
# RSE
Mlsrrse <- c(summary(Ml1)$sigma,summary(Ml2)$sigma,summary(Ml3)$sigma,summary(Ml4)$sigma,summary(Ml5)$sigma,summary(Ml6)$sigma,summary(Ml7)$sigma,summary(Ml8)$sigma,summary(Ml9)$sigma,summary(Ml10)$sigma)
# AIC
Mlsraic <- AIC(Ml1,Ml2,Ml3,Ml4,Ml5,Ml6,Ml7,Ml8,Ml9,Ml10)$AIC
# R2
Mlsrr2 <- c(Ml1r2,Ml2r2,Ml3r2,Ml4r2,Ml5r2,Ml6r2,Ml7r2,Ml8r2,Ml9r2,Ml10r2)
# show goodness-of-fit measurements
Mlsrrse## [1] 2.425928 2.260627 3.201058 2.388457 2.219373 2.441832 2.234935 2.263683
## [9] 2.221036 2.393161Mlsraic## [1] 418.9050 406.2020 468.8130 417.0743 403.8583 420.0811 406.0755 406.4451
## [9] 403.9930 416.4571Mlsrr2## [1] 0.7444385 0.7780796 0.5550345 0.7550876 0.7885359 0.7410768 0.7880249
## [8] 0.7774793 0.7882190 0.7512958After assessing the goodness-of-fit of each model, the best models for each biomass component were chosen as a base for Seemingly Unrelated Regression. This approach allows fulfillment of the logical assumption that the sum of the estimated biomass components (stem, branches, and foliage) should be equal to the estimated total aboveground biomass, assuring additivity of the biomass equations (Bronisz et al., 2016). We created a set of allometric equations that took into account all of the sample trees. Although the study design may have an impact on the model, we chose a simpler approach (non-linear models) instead of mixed-effects models to maximize model applicability.
#### Seemingly Unrelated Regression - SUR
Ms.formula<- SDB~b0*D^b1*H^b2
Mb.formula<- BDB~b3*D^b4*H^b5
Ml.formula<- LDB~b6+b7*D2+b8*H2
equations <- list(Ms.formula,Mb.formula,Ml.formula)
start.values <- c(b0=0.0005357 ,b1=1.6941255 ,b2=1.2137352
,b3=5.215e-06,b4=3.609e+00,b5=-1.532e+00,
b6=1.224e+00, b7=2.484e-04,b8 = -9.576e-03)
labels <- list("SDB","BDB","LDB")
SUR<-nlsystemfit(method ="SUR",equations,start.values,data=tree,eqnlabels=labels)## Warning in nlm(knls, estols$estimate, typsize = abs(estols$estimate), iterlim =
## maxiter, : NA/Inf replaced by maximum positive valueprint(SUR)## $eq
## $eq[[1]]
## $method
## [1] "SUR"
##
## $i
## [1] 1
##
## $eqnlabel
## [1] "SDB"
##
## $formula
## SDB ~ b0 * D^b1 * H^b2
##
## $b
## b0 b1 b2
## 0.0005442026 1.6196007186 1.3436476330
##
## $se
## b0
## [1,] 0.0001111865
## [2,] 0.0840125412
## [3,] 0.1184139329
##
## $t
## b0
## [1,] 4.894503
## [2,] 19.278083
## [3,] 11.347040
##
## $p
## b0
## [1,] 4.493348e-06
## [2,] 0.000000e+00
## [3,] 0.000000e+00
##
## $n
## [1] 90
##
## $k
## [1] 3
##
## $df
## [1] 87
##
## $predicted
## [,1]
## [1,] 70.765850
## [2,] 62.651178
## [3,] 85.365667
## [4,] 95.898858
## [5,] 170.099435
## [6,] 10.582751
## [7,] 25.610589
## [8,] 27.628456
## [9,] 37.874395
## [10,] 58.473363
## [11,] 13.460839
## [12,] 15.400986
## [13,] 20.151757
## [14,] 23.469379
## [15,] 30.450151
## [16,] 59.510185
## [17,] 78.859149
## [18,] 67.071910
## [19,] 98.125818
## [20,] 131.144558
## [21,] 49.077417
## [22,] 88.080449
## [23,] 75.226504
## [24,] 102.223122
## [25,] 122.750102
## [26,] 7.897053
## [27,] 8.898923
## [28,] 14.201582
## [29,] 26.736085
## [30,] 26.672386
## [31,] 180.152694
## [32,] 238.686160
## [33,] 341.885746
## [34,] 337.533606
## [35,] 449.418923
## [36,] 14.509558
## [37,] 23.265122
## [38,] 41.483569
## [39,] 38.078980
## [40,] 45.794306
## [41,] 135.697726
## [42,] 139.023497
## [43,] 173.188127
## [44,] 186.099215
## [45,] 262.241172
## [46,] 28.478480
## [47,] 32.573656
## [48,] 51.312298
## [49,] 51.037659
## [50,] 63.407013
## [51,] 11.773959
## [52,] 21.602771
## [53,] 31.869364
## [54,] 25.312302
## [55,] 32.038732
## [56,] 19.099401
## [57,] 32.625473
## [58,] 32.506319
## [59,] 40.424627
## [60,] 51.026690
## [61,] 78.167738
## [62,] 102.908010
## [63,] 156.227479
## [64,] 148.606696
## [65,] 124.181099
## [66,] 100.697243
## [67,] 156.319434
## [68,] 139.571103
## [69,] 202.640308
## [70,] 253.747612
## [71,] 163.493203
## [72,] 205.621213
## [73,] 262.807140
## [74,] 329.102894
## [75,] 440.831764
## [76,] 24.253012
## [77,] 26.861029
## [78,] 51.294620
## [79,] 44.347524
## [80,] 53.692288
## [81,] 68.105139
## [82,] 72.847366
## [83,] 102.991212
## [84,] 114.443140
## [85,] 155.301664
## [86,] 79.799693
## [87,] 92.841124
## [88,] 122.116621
## [89,] 159.693612
## [90,] 234.321036
##
## $residuals
## [,1]
## [1,] -7.57591049
## [2,] -1.25240387
## [3,] 9.40189307
## [4,] 7.01217756
## [5,] -2.93110263
## [6,] -0.15494258
## [7,] -1.83113786
## [8,] 2.75600240
## [9,] 6.17236329
## [10,] 1.49527998
## [11,] 0.71989343
## [12,] 0.49539337
## [13,] 0.22370133
## [14,] 2.77279181
## [15,] 1.25248999
## [16,] 7.68710977
## [17,] 17.07530299
## [18,] 5.04126157
## [19,] 5.98952501
## [20,] 13.75585114
## [21,] -3.15354095
## [22,] -10.23721919
## [23,] 2.37195923
## [24,] -19.15880250
## [25,] -4.59565376
## [26,] 1.52897718
## [27,] 1.24982131
## [28,] 0.72807795
## [29,] -0.77524213
## [30,] 3.18149329
## [31,] 7.16918598
## [32,] -4.00652142
## [33,] -43.00073339
## [34,] 19.07408724
## [35,] 2.67261458
## [36,] -2.44415189
## [37,] 0.08876253
## [38,] 1.54906801
## [39,] -0.72733215
## [40,] 4.21561249
## [41,] -6.74353181
## [42,] -7.70495194
## [43,] 12.54608607
## [44,] 8.45466342
## [45,] 44.35333711
## [46,] 0.20032798
## [47,] -11.74553200
## [48,] -7.51000759
## [49,] -7.50789883
## [50,] -0.87135952
## [51,] -3.48596246
## [52,] -3.36045692
## [53,] -1.63105925
## [54,] -3.45781516
## [55,] 1.98616309
## [56,] 0.13955587
## [57,] 1.77166272
## [58,] -0.97739159
## [59,] -0.14374547
## [60,] -2.01801101
## [61,] 9.08958595
## [62,] 1.90802322
## [63,] -4.23384159
## [64,] -13.50058382
## [65,] -1.51474474
## [66,] -6.41256561
## [67,] -3.82589318
## [68,] -7.89969963
## [69,] -9.71073985
## [70,] -12.96556786
## [71,] -8.75734121
## [72,] 12.45930755
## [73,] 26.60254256
## [74,] 42.59822839
## [75,] -32.99032923
## [76,] -5.37710031
## [77,] 0.54551537
## [78,] -2.94087815
## [79,] -6.94265936
## [80,] -4.73182530
## [81,] 3.79476361
## [82,] -6.43636087
## [83,] -1.21166018
## [84,] 2.12034684
## [85,] 16.42646858
## [86,] -0.84830669
## [87,] -8.72078653
## [88,] 7.71477868
## [89,] 2.38240802
## [90,] -40.62264128
##
## $ssr
## [1] 13052.35
##
## $mse
## [1] 150.0271
##
## $s2
## [1] 150.0271
##
## $rmse
## [1] 12.24855
##
## $s
## [1] 12.24855
##
## $covb
## b0 b1 b2
## b0 1.236243e-08 -6.188508e-06 3.726210e-06
## b1 -6.188508e-06 7.058107e-03 -9.009862e-03
## b2 3.726210e-06 -9.009862e-03 1.402186e-02
##
## $r2
## [,1]
## [1,] 0.9839898
##
## $adjr2
## [,1]
## [1,] 0.9836218
##
## attr(,"class")
## [1] "nlsystemfit.equation"
##
## $eq[[2]]
## $method
## [1] "SUR"
##
## $i
## [1] 2
##
## $eqnlabel
## [1] "BDB"
##
## $formula
## BDB ~ b3 * D^b4 * H^b5
##
## $b
## b3 b4 b5
## 5.655787e-06 3.336614e+00 -1.063237e+00
##
## $se
## b3
## [1,] 4.790170e-06
## [2,] 3.296900e-01
## [3,] 4.000897e-01
##
## $t
## b3
## [1,] 1.180707
## [2,] 10.120457
## [3,] -2.657498
##
## $p
## b3
## [1,] 2.409370e-01
## [2,] 2.220446e-16
## [3,] 9.366272e-03
##
## $n
## [1] 90
##
## $k
## [1] 3
##
## $df
## [1] 87
##
## $predicted
## [,1]
## [1,] 6.8237134
## [2,] 5.8923813
## [3,] 8.0579198
## [4,] 10.6586096
## [5,] 18.5008708
## [6,] 0.7094252
## [7,] 1.6995151
## [8,] 2.0186924
## [9,] 2.8486110
## [10,] 4.9753997
## [11,] 1.0330830
## [12,] 1.2373085
## [13,] 1.6913901
## [14,] 2.0961422
## [15,] 2.9995853
## [16,] 6.5296377
## [17,] 8.1152819
## [18,] 6.9110909
## [19,] 10.7364787
## [20,] 17.2666031
## [21,] 4.5175094
## [22,] 8.0805006
## [23,] 6.6962188
## [24,] 13.9594498
## [25,] 11.7927689
## [26,] 0.8943515
## [27,] 1.6111632
## [28,] 1.5543752
## [29,] 2.5842610
## [30,] 4.3967474
## [31,] 14.3309954
## [32,] 17.2976956
## [33,] 28.2829696
## [34,] 30.0231140
## [35,] 34.9500383
## [36,] 1.0942886
## [37,] 1.5860739
## [38,] 4.0277221
## [39,] 2.6175789
## [40,] 6.4427569
## [41,] 11.1562706
## [42,] 9.9405805
## [43,] 16.7892346
## [44,] 15.5319225
## [45,] 27.9649418
## [46,] 1.3515671
## [47,] 1.6172101
## [48,] 3.1633656
## [49,] 2.5312471
## [50,] 3.8174543
## [51,] 0.5691881
## [52,] 1.0566789
## [53,] 2.2489459
## [54,] 1.5579910
## [55,] 2.2395330
## [56,] 0.7832458
## [57,] 1.6911309
## [58,] 1.7749864
## [59,] 3.1639825
## [60,] 3.9454844
## [61,] 4.7934680
## [62,] 9.0915324
## [63,] 12.8741669
## [64,] 10.6316680
## [65,] 8.9111384
## [66,] 5.3163462
## [67,] 8.8415952
## [68,] 10.6532404
## [69,] 16.4410294
## [70,] 22.3391933
## [71,] 9.6267759
## [72,] 14.9655981
## [73,] 23.3734357
## [74,] 30.1299328
## [75,] 41.1890793
## [76,] 0.9129331
## [77,] 2.0641961
## [78,] 2.6205506
## [79,] 2.4157096
## [80,] 4.2074038
## [81,] 4.2672782
## [82,] 4.4621323
## [83,] 6.3879365
## [84,] 8.1425825
## [85,] 15.4362923
## [86,] 6.2048656
## [87,] 12.0970109
## [88,] 12.1271608
## [89,] 16.1068588
## [90,] 30.0412307
##
## $residuals
## [,1]
## [1,] -3.93179912
## [2,] -2.38878922
## [3,] -2.90641972
## [4,] 0.91806971
## [5,] 1.03963576
## [6,] -0.04094521
## [7,] 0.65543318
## [8,] 0.49313933
## [9,] 3.05113825
## [10,] 1.75517555
## [11,] -0.24433992
## [12,] 0.37655643
## [13,] 0.57168453
## [14,] 0.79106527
## [15,] 0.18548701
## [16,] 0.93746464
## [17,] -0.29860302
## [18,] 2.15193252
## [19,] -1.81905807
## [20,] -5.72307724
## [21,] -0.12789087
## [22,] -1.74997131
## [23,] 2.07407358
## [24,] -3.69335340
## [25,] -5.58892645
## [26,] 0.39977797
## [27,] 0.15931049
## [28,] 1.59655083
## [29,] 2.73894352
## [30,] 6.72928838
## [31,] -1.73043896
## [32,] -1.37658623
## [33,] -5.68107446
## [34,] -0.20006655
## [35,] -2.82015187
## [36,] 0.08857148
## [37,] 1.25118378
## [38,] 0.78851809
## [39,] 0.11453972
## [40,] 4.95780927
## [41,] -0.13623192
## [42,] -1.84676366
## [43,] 2.37446331
## [44,] -2.00445225
## [45,] 27.44071462
## [46,] 0.66511710
## [47,] -0.13228736
## [48,] -0.95760577
## [49,] -1.22205958
## [50,] 1.32470643
## [51,] 0.31899675
## [52,] 1.53730948
## [53,] 1.81015196
## [54,] 1.03201206
## [55,] 1.95086216
## [56,] 0.08995252
## [57,] 0.28329583
## [58,] -0.25993330
## [59,] 0.84916211
## [60,] -0.80478499
## [61,] -0.28383349
## [62,] -1.98730629
## [63,] -5.54566836
## [64,] -1.34068536
## [65,] -3.47804891
## [66,] -0.82677054
## [67,] -3.20528016
## [68,] -6.41339000
## [69,] -4.11174278
## [70,] -5.48835886
## [71,] -2.03643009
## [72,] 2.49895022
## [73,] 3.11501732
## [74,] 5.24272446
## [75,] -4.21343249
## [76,] 0.07492994
## [77,] 1.16873169
## [78,] 0.49356427
## [79,] 0.57244893
## [80,] 1.44727509
## [81,] -0.77327555
## [82,] -1.58927557
## [83,] 1.50335197
## [84,] 4.48690888
## [85,] 16.41059882
## [86,] 0.02038738
## [87,] -2.38969998
## [88,] 5.99270926
## [89,] 1.92564332
## [90,] -10.03482353
##
## $ssr
## [1] 1708.034
##
## $mse
## [1] 19.63258
##
## $s2
## [1] 19.63258
##
## $rmse
## [1] 4.430867
##
## $s
## [1] 4.430867
##
## $covb
## b3 b4 b5
## b3 2.294573e-11 -1.309892e-06 1.048941e-06
## b4 -1.309892e-06 1.086955e-01 -1.215059e-01
## b5 1.048941e-06 -1.215059e-01 1.600717e-01
##
## $r2
## [,1]
## [1,] 0.7925199
##
## $adjr2
## [,1]
## [1,] 0.7877503
##
## attr(,"class")
## [1] "nlsystemfit.equation"
##
## $eq[[3]]
## $method
## [1] "SUR"
##
## $i
## [1] 3
##
## $eqnlabel
## [1] "LDB"
##
## $formula
## LDB ~ b6 + b7 * D2 + b8 * H2
##
## $b
## b6 b7 b8
## 0.4402074191 0.0002215587 -0.0038177178
##
## $se
## b6
## [1,] 4.347092e-01
## [2,] 2.328384e-05
## [3,] 3.689753e-03
##
## $t
## b6
## [1,] 1.012648
## [2,] 9.515559
## [3,] -1.034681
##
## $p
## b6
## [1,] 3.140364e-01
## [2,] 3.996803e-15
## [3,] 3.036854e-01
##
## $n
## [1] 90
##
## $k
## [1] 3
##
## $df
## [1] 87
##
## $predicted
## [,1]
## [1,] 5.024951
## [2,] 4.528837
## [3,] 5.713849
## [4,] 6.812000
## [5,] 10.425170
## [6,] 1.154117
## [7,] 1.965902
## [8,] 2.190962
## [9,] 2.771335
## [10,] 4.069782
## [11,] 1.425249
## [12,] 1.584754
## [13,] 1.923454
## [14,] 2.199712
## [15,] 2.768779
## [16,] 4.739770
## [17,] 5.643803
## [18,] 5.006308
## [19,] 6.873709
## [20,] 9.432787
## [21,] 3.743200
## [22,] 5.758437
## [23,] 5.030931
## [24,] 7.953822
## [25,] 7.596209
## [26,] 1.274537
## [27,] 1.673778
## [28,] 1.768245
## [29,] 2.506453
## [30,] 3.275138
## [31,] 9.185062
## [32,] 10.839754
## [33,] 15.472724
## [34,] 15.936154
## [35,] 18.538947
## [36,] 1.476412
## [37,] 1.875293
## [38,] 3.415300
## [39,] 2.636120
## [40,] 4.450040
## [41,] 7.521150
## [42,] 7.067526
## [43,] 9.940581
## [44,] 9.677517
## [45,] 14.418168
## [46,] 1.691775
## [47,] 1.912857
## [48,] 3.048793
## [49,] 2.624646
## [50,] 3.506879
## [51,] 1.005826
## [52,] 1.444451
## [53,] 2.365913
## [54,] 1.862193
## [55,] 2.361074
## [56,] 1.177250
## [57,] 1.973282
## [58,] 2.039180
## [59,] 2.970484
## [60,] 3.486603
## [61,] 4.147791
## [62,] 6.342048
## [63,] 8.393253
## [64,] 7.442917
## [65,] 6.493608
## [66,] 4.557811
## [67,] 6.701221
## [68,] 7.365161
## [69,] 10.175539
## [70,] 12.682795
## [71,] 7.119145
## [72,] 9.672687
## [73,] 13.106863
## [74,] 15.869094
## [75,] 20.166338
## [76,] 1.274371
## [77,] 2.213726
## [78,] 2.689624
## [79,] 2.530695
## [80,] 3.646542
## [81,] 3.793168
## [82,] 3.931489
## [83,] 5.148529
## [84,] 6.069018
## [85,] 9.264768
## [86,] 4.866904
## [87,] 7.221487
## [88,] 7.703844
## [89,] 9.540312
## [90,] 14.532632
##
## $residuals
## [,1]
## [1,] -2.43303618
## [2,] -1.10414076
## [3,] -1.82129820
## [4,] 0.41167713
## [5,] -1.26457465
## [6,] -0.59628017
## [7,] -0.66128351
## [8,] -0.39148921
## [9,] 1.86230858
## [10,] 0.10499639
## [11,] -0.80223582
## [12,] -0.22191243
## [13,] 0.04293155
## [14,] 0.65252079
## [15,] -0.23778750
## [16,] 0.27447404
## [17,] 0.93800058
## [18,] 0.05661833
## [19,] -1.70368321
## [20,] -0.62134827
## [21,] -1.08925744
## [22,] -2.05072238
## [23,] 3.78952563
## [24,] -3.61189339
## [25,] -1.43945244
## [26,] -0.41883753
## [27,] -0.23211010
## [28,] 0.93837052
## [29,] 1.78134713
## [30,] 3.51465323
## [31,] -3.43056496
## [32,] -1.41069760
## [33,] -2.76470559
## [34,] 2.33468154
## [35,] -2.47040031
## [36,] -0.95900467
## [37,] 1.34987690
## [38,] 2.51332171
## [39,] 1.41690089
## [40,] 3.57839126
## [41,] -1.17714785
## [42,] -1.78085750
## [43,] 0.53162058
## [44,] -1.21884679
## [45,] 13.33914160
## [46,] 0.09928473
## [47,] -0.55297846
## [48,] -0.79338429
## [49,] 0.42679005
## [50,] 2.70378829
## [51,] -0.47976136
## [52,] 0.05522547
## [53,] 0.90824992
## [54,] 0.13999323
## [55,] 1.03871622
## [56,] -0.29099487
## [57,] 0.04453121
## [58,] 0.14480655
## [59,] -0.18511051
## [60,] -0.16314586
## [61,] 1.21798324
## [62,] -1.10223094
## [63,] -1.09613641
## [64,] 0.74429666
## [65,] -1.52249964
## [66,] -1.49034898
## [67,] -1.19033778
## [68,] -3.12872415
## [69,] -1.27199025
## [70,] -0.16684863
## [71,] -1.28379818
## [72,] 3.12997556
## [73,] 0.80144362
## [74,] 1.72469472
## [75,] -3.47636928
## [76,] -0.64461902
## [77,] 0.51519905
## [78,] 1.70184539
## [79,] 1.30659251
## [80,] 0.01441570
## [81,] 0.92897463
## [82,] -0.83645112
## [83,] 0.74088517
## [84,] 3.71130053
## [85,] 5.85670113
## [86,] -0.56658934
## [87,] -1.37305939
## [88,] -0.78033737
## [89,] -2.16957268
## [90,] -0.34159378
##
## $ssr
## [1] 437.5295
##
## $mse
## [1] 5.029075
##
## $s2
## [1] 5.029075
##
## $rmse
## [1] 2.24256
##
## $s
## [1] 2.24256
##
## $covb
## b6 b7 b8
## b6 1.889721e-01 4.969270e-06 -1.195253e-03
## b7 4.969270e-06 5.421372e-10 -7.858059e-08
## b8 -1.195253e-03 -7.858059e-08 1.361428e-05
##
## $r2
## [,1]
## [1,] 0.7840943
##
## $adjr2
## [,1]
## [1,] 0.779131
##
## attr(,"class")
## [1] "nlsystemfit.equation"
##
##
## $resids
## [,1] [,2] [,3]
## [1,] -7.57591049 -3.93179912 -2.43303618
## [2,] -1.25240387 -2.38878922 -1.10414076
## [3,] 9.40189307 -2.90641972 -1.82129820
## [4,] 7.01217756 0.91806971 0.41167713
## [5,] -2.93110263 1.03963576 -1.26457465
## [6,] -0.15494258 -0.04094521 -0.59628017
## [7,] -1.83113786 0.65543318 -0.66128351
## [8,] 2.75600240 0.49313933 -0.39148921
## [9,] 6.17236329 3.05113825 1.86230858
## [10,] 1.49527998 1.75517555 0.10499639
## [11,] 0.71989343 -0.24433992 -0.80223582
## [12,] 0.49539337 0.37655643 -0.22191243
## [13,] 0.22370133 0.57168453 0.04293155
## [14,] 2.77279181 0.79106527 0.65252079
## [15,] 1.25248999 0.18548701 -0.23778750
## [16,] 7.68710977 0.93746464 0.27447404
## [17,] 17.07530299 -0.29860302 0.93800058
## [18,] 5.04126157 2.15193252 0.05661833
## [19,] 5.98952501 -1.81905807 -1.70368321
## [20,] 13.75585114 -5.72307724 -0.62134827
## [21,] -3.15354095 -0.12789087 -1.08925744
## [22,] -10.23721919 -1.74997131 -2.05072238
## [23,] 2.37195923 2.07407358 3.78952563
## [24,] -19.15880250 -3.69335340 -3.61189339
## [25,] -4.59565376 -5.58892645 -1.43945244
## [26,] 1.52897718 0.39977797 -0.41883753
## [27,] 1.24982131 0.15931049 -0.23211010
## [28,] 0.72807795 1.59655083 0.93837052
## [29,] -0.77524213 2.73894352 1.78134713
## [30,] 3.18149329 6.72928838 3.51465323
## [31,] 7.16918598 -1.73043896 -3.43056496
## [32,] -4.00652142 -1.37658623 -1.41069760
## [33,] -43.00073339 -5.68107446 -2.76470559
## [34,] 19.07408724 -0.20006655 2.33468154
## [35,] 2.67261458 -2.82015187 -2.47040031
## [36,] -2.44415189 0.08857148 -0.95900467
## [37,] 0.08876253 1.25118378 1.34987690
## [38,] 1.54906801 0.78851809 2.51332171
## [39,] -0.72733215 0.11453972 1.41690089
## [40,] 4.21561249 4.95780927 3.57839126
## [41,] -6.74353181 -0.13623192 -1.17714785
## [42,] -7.70495194 -1.84676366 -1.78085750
## [43,] 12.54608607 2.37446331 0.53162058
## [44,] 8.45466342 -2.00445225 -1.21884679
## [45,] 44.35333711 27.44071462 13.33914160
## [46,] 0.20032798 0.66511710 0.09928473
## [47,] -11.74553200 -0.13228736 -0.55297846
## [48,] -7.51000759 -0.95760577 -0.79338429
## [49,] -7.50789883 -1.22205958 0.42679005
## [50,] -0.87135952 1.32470643 2.70378829
## [51,] -3.48596246 0.31899675 -0.47976136
## [52,] -3.36045692 1.53730948 0.05522547
## [53,] -1.63105925 1.81015196 0.90824992
## [54,] -3.45781516 1.03201206 0.13999323
## [55,] 1.98616309 1.95086216 1.03871622
## [56,] 0.13955587 0.08995252 -0.29099487
## [57,] 1.77166272 0.28329583 0.04453121
## [58,] -0.97739159 -0.25993330 0.14480655
## [59,] -0.14374547 0.84916211 -0.18511051
## [60,] -2.01801101 -0.80478499 -0.16314586
## [61,] 9.08958595 -0.28383349 1.21798324
## [62,] 1.90802322 -1.98730629 -1.10223094
## [63,] -4.23384159 -5.54566836 -1.09613641
## [64,] -13.50058382 -1.34068536 0.74429666
## [65,] -1.51474474 -3.47804891 -1.52249964
## [66,] -6.41256561 -0.82677054 -1.49034898
## [67,] -3.82589318 -3.20528016 -1.19033778
## [68,] -7.89969963 -6.41339000 -3.12872415
## [69,] -9.71073985 -4.11174278 -1.27199025
## [70,] -12.96556786 -5.48835886 -0.16684863
## [71,] -8.75734121 -2.03643009 -1.28379818
## [72,] 12.45930755 2.49895022 3.12997556
## [73,] 26.60254256 3.11501732 0.80144362
## [74,] 42.59822839 5.24272446 1.72469472
## [75,] -32.99032923 -4.21343249 -3.47636928
## [76,] -5.37710031 0.07492994 -0.64461902
## [77,] 0.54551537 1.16873169 0.51519905
## [78,] -2.94087815 0.49356427 1.70184539
## [79,] -6.94265936 0.57244893 1.30659251
## [80,] -4.73182530 1.44727509 0.01441570
## [81,] 3.79476361 -0.77327555 0.92897463
## [82,] -6.43636087 -1.58927557 -0.83645112
## [83,] -1.21166018 1.50335197 0.74088517
## [84,] 2.12034684 4.48690888 3.71130053
## [85,] 16.42646858 16.41059882 5.85670113
## [86,] -0.84830669 0.02038738 -0.56658934
## [87,] -8.72078653 -2.38969998 -1.37305939
## [88,] 7.71477868 5.99270926 -0.78033737
## [89,] 2.38240802 1.92564332 -2.16957268
## [90,] -40.62264128 -10.03482353 -0.34159378
##
## $method
## [1] "SUR"
##
## $n
## [1] 270
##
## $k
## [1] 9
##
## $b
## b0 b1 b2 b3 b4
## 5.442026e-04 1.619601e+00 1.343648e+00 5.655787e-06 3.336614e+00
## b5 b6 b7 b8
## -1.063237e+00 4.402074e-01 2.215587e-04 -3.817718e-03
##
## $se
## b0 b1 b2 b3 b4 b5
## 1.111865e-04 8.401254e-02 1.184139e-01 4.790170e-06 3.296900e-01 4.000897e-01
## b6 b7 b8
## 4.347092e-01 2.328384e-05 3.689753e-03
##
## $t
## b0 b1 b2 b3 b4 b5 b6 b7
## 4.894503 19.278083 11.347040 1.180707 10.120457 -2.657498 1.012648 9.515559
## b8
## -1.034681
##
## $p
## b0 b1 b2 b3 b4 b5
## 1.727817e-06 0.000000e+00 0.000000e+00 2.387940e-01 0.000000e+00 8.357930e-03
## b6 b7 b8
## 3.121662e-01 0.000000e+00 3.017756e-01
##
## $solvtol
## [1] 2.220446e-16
##
## $covb
## b0 b1 b2 b3 b4
## b0 1.236243e-08 -6.188508e-06 3.726210e-06 2.862436e-10 -1.584492e-05
## b1 -6.188508e-06 7.058107e-03 -9.009862e-03 -1.646716e-07 1.494151e-02
## b2 3.726210e-06 -9.009862e-03 1.402186e-02 1.262198e-07 -1.753219e-02
## b3 2.862436e-10 -1.646716e-07 1.262198e-07 2.294573e-11 -1.309892e-06
## b4 -1.584492e-05 1.494151e-02 -1.753219e-02 -1.309892e-06 1.086955e-01
## b5 1.211757e-05 -1.747594e-02 2.438185e-02 1.048941e-06 -1.215059e-01
## b6 3.104501e-06 7.322139e-03 -1.518156e-02 1.468942e-07 4.080445e-02
## b7 -7.265489e-10 8.464497e-07 -1.079246e-06 -5.290996e-11 5.031038e-06
## b8 9.144138e-08 -1.282443e-04 1.767151e-04 7.200291e-09 -8.055952e-04
## b5 b6 b7 b8
## b0 1.211757e-05 3.104501e-06 -7.265489e-10 9.144138e-08
## b1 -1.747594e-02 7.322139e-03 8.464497e-07 -1.282443e-04
## b2 2.438185e-02 -1.518156e-02 -1.079246e-06 1.767151e-04
## b3 1.048941e-06 1.468942e-07 -5.290996e-11 7.200291e-09
## b4 -1.215059e-01 4.080445e-02 5.031038e-06 -8.055952e-04
## b5 1.600717e-01 -8.302141e-02 -5.973241e-06 1.042968e-03
## b6 -8.302141e-02 1.889721e-01 4.969270e-06 -1.195253e-03
## b7 -5.973241e-06 4.969270e-06 5.421372e-10 -7.858059e-08
## b8 1.042968e-03 -1.195253e-03 -7.858059e-08 1.361428e-05
##
## $rcov
## [,1] [,2] [,3]
## [1,] 150.02707 35.245652 15.550452
## [2,] 35.24565 19.632580 8.264076
## [3,] 15.55045 8.264076 5.029075
##
## $rcor
## [,1] [,2] [,3]
## [1,] 1.0000000 0.6508509 0.5670606
## [2,] 0.6508509 1.0000000 0.8314927
## [3,] 0.5670606 0.8314927 1.0000000
##
## $drcov
## [1] 2630.602
##
## $rcovest
## [,1] [,2] [,3]
## [1,] 148.21532 34.663680 15.022083
## [2,] 34.66368 19.383487 7.963019
## [3,] 15.02208 7.963019 4.925620
##
## $g
## [1] 3
##
## $nlmest
## $nlmest$minimum
## [1] 257.3351
##
## $nlmest$estimate
## [1] 5.442026e-04 1.619601e+00 1.343648e+00 5.655787e-06 3.336614e+00
## [6] -1.063237e+00 4.402074e-01 2.215587e-04 -3.817718e-03
##
## $nlmest$gradient
## [1] -1.143399e-02 -2.960711e-05 -1.561974e-05 1.399352e+00 4.365029e-05
## [6] 2.431022e-05 -2.221540e-07 -4.440257e-02 -3.143265e-04
##
## $nlmest$code
## [1] 1
##
## $nlmest$iterations
## [1] 113
##
##
## attr(,"class")
## [1] "nlsystemfit.system"summary(SUR)## Length Class Mode
## eq 3 -none- list
## resids 270 -none- numeric
## method 1 -none- character
## n 1 -none- numeric
## k 1 -none- numeric
## b 9 -none- numeric
## se 9 -none- numeric
## t 9 -none- numeric
## p 9 -none- numeric
## solvtol 1 -none- numeric
## covb 81 -none- numeric
## rcov 9 -none- numeric
## rcor 9 -none- numeric
## drcov 1 -none- numeric
## rcovest 9 -none- numeric
## g 1 -none- numeric
## nlmest 5 -none- listSUR$eq[3]## [[1]]
## $method
## [1] "SUR"
##
## $i
## [1] 3
##
## $eqnlabel
## [1] "LDB"
##
## $formula
## LDB ~ b6 + b7 * D2 + b8 * H2
##
## $b
## b6 b7 b8
## 0.4402074191 0.0002215587 -0.0038177178
##
## $se
## b6
## [1,] 4.347092e-01
## [2,] 2.328384e-05
## [3,] 3.689753e-03
##
## $t
## b6
## [1,] 1.012648
## [2,] 9.515559
## [3,] -1.034681
##
## $p
## b6
## [1,] 3.140364e-01
## [2,] 3.996803e-15
## [3,] 3.036854e-01
##
## $n
## [1] 90
##
## $k
## [1] 3
##
## $df
## [1] 87
##
## $predicted
## [,1]
## [1,] 5.024951
## [2,] 4.528837
## [3,] 5.713849
## [4,] 6.812000
## [5,] 10.425170
## [6,] 1.154117
## [7,] 1.965902
## [8,] 2.190962
## [9,] 2.771335
## [10,] 4.069782
## [11,] 1.425249
## [12,] 1.584754
## [13,] 1.923454
## [14,] 2.199712
## [15,] 2.768779
## [16,] 4.739770
## [17,] 5.643803
## [18,] 5.006308
## [19,] 6.873709
## [20,] 9.432787
## [21,] 3.743200
## [22,] 5.758437
## [23,] 5.030931
## [24,] 7.953822
## [25,] 7.596209
## [26,] 1.274537
## [27,] 1.673778
## [28,] 1.768245
## [29,] 2.506453
## [30,] 3.275138
## [31,] 9.185062
## [32,] 10.839754
## [33,] 15.472724
## [34,] 15.936154
## [35,] 18.538947
## [36,] 1.476412
## [37,] 1.875293
## [38,] 3.415300
## [39,] 2.636120
## [40,] 4.450040
## [41,] 7.521150
## [42,] 7.067526
## [43,] 9.940581
## [44,] 9.677517
## [45,] 14.418168
## [46,] 1.691775
## [47,] 1.912857
## [48,] 3.048793
## [49,] 2.624646
## [50,] 3.506879
## [51,] 1.005826
## [52,] 1.444451
## [53,] 2.365913
## [54,] 1.862193
## [55,] 2.361074
## [56,] 1.177250
## [57,] 1.973282
## [58,] 2.039180
## [59,] 2.970484
## [60,] 3.486603
## [61,] 4.147791
## [62,] 6.342048
## [63,] 8.393253
## [64,] 7.442917
## [65,] 6.493608
## [66,] 4.557811
## [67,] 6.701221
## [68,] 7.365161
## [69,] 10.175539
## [70,] 12.682795
## [71,] 7.119145
## [72,] 9.672687
## [73,] 13.106863
## [74,] 15.869094
## [75,] 20.166338
## [76,] 1.274371
## [77,] 2.213726
## [78,] 2.689624
## [79,] 2.530695
## [80,] 3.646542
## [81,] 3.793168
## [82,] 3.931489
## [83,] 5.148529
## [84,] 6.069018
## [85,] 9.264768
## [86,] 4.866904
## [87,] 7.221487
## [88,] 7.703844
## [89,] 9.540312
## [90,] 14.532632
##
## $residuals
## [,1]
## [1,] -2.43303618
## [2,] -1.10414076
## [3,] -1.82129820
## [4,] 0.41167713
## [5,] -1.26457465
## [6,] -0.59628017
## [7,] -0.66128351
## [8,] -0.39148921
## [9,] 1.86230858
## [10,] 0.10499639
## [11,] -0.80223582
## [12,] -0.22191243
## [13,] 0.04293155
## [14,] 0.65252079
## [15,] -0.23778750
## [16,] 0.27447404
## [17,] 0.93800058
## [18,] 0.05661833
## [19,] -1.70368321
## [20,] -0.62134827
## [21,] -1.08925744
## [22,] -2.05072238
## [23,] 3.78952563
## [24,] -3.61189339
## [25,] -1.43945244
## [26,] -0.41883753
## [27,] -0.23211010
## [28,] 0.93837052
## [29,] 1.78134713
## [30,] 3.51465323
## [31,] -3.43056496
## [32,] -1.41069760
## [33,] -2.76470559
## [34,] 2.33468154
## [35,] -2.47040031
## [36,] -0.95900467
## [37,] 1.34987690
## [38,] 2.51332171
## [39,] 1.41690089
## [40,] 3.57839126
## [41,] -1.17714785
## [42,] -1.78085750
## [43,] 0.53162058
## [44,] -1.21884679
## [45,] 13.33914160
## [46,] 0.09928473
## [47,] -0.55297846
## [48,] -0.79338429
## [49,] 0.42679005
## [50,] 2.70378829
## [51,] -0.47976136
## [52,] 0.05522547
## [53,] 0.90824992
## [54,] 0.13999323
## [55,] 1.03871622
## [56,] -0.29099487
## [57,] 0.04453121
## [58,] 0.14480655
## [59,] -0.18511051
## [60,] -0.16314586
## [61,] 1.21798324
## [62,] -1.10223094
## [63,] -1.09613641
## [64,] 0.74429666
## [65,] -1.52249964
## [66,] -1.49034898
## [67,] -1.19033778
## [68,] -3.12872415
## [69,] -1.27199025
## [70,] -0.16684863
## [71,] -1.28379818
## [72,] 3.12997556
## [73,] 0.80144362
## [74,] 1.72469472
## [75,] -3.47636928
## [76,] -0.64461902
## [77,] 0.51519905
## [78,] 1.70184539
## [79,] 1.30659251
## [80,] 0.01441570
## [81,] 0.92897463
## [82,] -0.83645112
## [83,] 0.74088517
## [84,] 3.71130053
## [85,] 5.85670113
## [86,] -0.56658934
## [87,] -1.37305939
## [88,] -0.78033737
## [89,] -2.16957268
## [90,] -0.34159378
##
## $ssr
## [1] 437.5295
##
## $mse
## [1] 5.029075
##
## $s2
## [1] 5.029075
##
## $rmse
## [1] 2.24256
##
## $s
## [1] 2.24256
##
## $covb
## b6 b7 b8
## b6 1.889721e-01 4.969270e-06 -1.195253e-03
## b7 4.969270e-06 5.421372e-10 -7.858059e-08
## b8 -1.195253e-03 -7.858059e-08 1.361428e-05
##
## $r2
## [,1]
## [1,] 0.7840943
##
## $adjr2
## [,1]
## [1,] 0.779131
##
## attr(,"class")
## [1] "nlsystemfit.equation"SUR$resids## [,1] [,2] [,3]
## [1,] -7.57591049 -3.93179912 -2.43303618
## [2,] -1.25240387 -2.38878922 -1.10414076
## [3,] 9.40189307 -2.90641972 -1.82129820
## [4,] 7.01217756 0.91806971 0.41167713
## [5,] -2.93110263 1.03963576 -1.26457465
## [6,] -0.15494258 -0.04094521 -0.59628017
## [7,] -1.83113786 0.65543318 -0.66128351
## [8,] 2.75600240 0.49313933 -0.39148921
## [9,] 6.17236329 3.05113825 1.86230858
## [10,] 1.49527998 1.75517555 0.10499639
## [11,] 0.71989343 -0.24433992 -0.80223582
## [12,] 0.49539337 0.37655643 -0.22191243
## [13,] 0.22370133 0.57168453 0.04293155
## [14,] 2.77279181 0.79106527 0.65252079
## [15,] 1.25248999 0.18548701 -0.23778750
## [16,] 7.68710977 0.93746464 0.27447404
## [17,] 17.07530299 -0.29860302 0.93800058
## [18,] 5.04126157 2.15193252 0.05661833
## [19,] 5.98952501 -1.81905807 -1.70368321
## [20,] 13.75585114 -5.72307724 -0.62134827
## [21,] -3.15354095 -0.12789087 -1.08925744
## [22,] -10.23721919 -1.74997131 -2.05072238
## [23,] 2.37195923 2.07407358 3.78952563
## [24,] -19.15880250 -3.69335340 -3.61189339
## [25,] -4.59565376 -5.58892645 -1.43945244
## [26,] 1.52897718 0.39977797 -0.41883753
## [27,] 1.24982131 0.15931049 -0.23211010
## [28,] 0.72807795 1.59655083 0.93837052
## [29,] -0.77524213 2.73894352 1.78134713
## [30,] 3.18149329 6.72928838 3.51465323
## [31,] 7.16918598 -1.73043896 -3.43056496
## [32,] -4.00652142 -1.37658623 -1.41069760
## [33,] -43.00073339 -5.68107446 -2.76470559
## [34,] 19.07408724 -0.20006655 2.33468154
## [35,] 2.67261458 -2.82015187 -2.47040031
## [36,] -2.44415189 0.08857148 -0.95900467
## [37,] 0.08876253 1.25118378 1.34987690
## [38,] 1.54906801 0.78851809 2.51332171
## [39,] -0.72733215 0.11453972 1.41690089
## [40,] 4.21561249 4.95780927 3.57839126
## [41,] -6.74353181 -0.13623192 -1.17714785
## [42,] -7.70495194 -1.84676366 -1.78085750
## [43,] 12.54608607 2.37446331 0.53162058
## [44,] 8.45466342 -2.00445225 -1.21884679
## [45,] 44.35333711 27.44071462 13.33914160
## [46,] 0.20032798 0.66511710 0.09928473
## [47,] -11.74553200 -0.13228736 -0.55297846
## [48,] -7.51000759 -0.95760577 -0.79338429
## [49,] -7.50789883 -1.22205958 0.42679005
## [50,] -0.87135952 1.32470643 2.70378829
## [51,] -3.48596246 0.31899675 -0.47976136
## [52,] -3.36045692 1.53730948 0.05522547
## [53,] -1.63105925 1.81015196 0.90824992
## [54,] -3.45781516 1.03201206 0.13999323
## [55,] 1.98616309 1.95086216 1.03871622
## [56,] 0.13955587 0.08995252 -0.29099487
## [57,] 1.77166272 0.28329583 0.04453121
## [58,] -0.97739159 -0.25993330 0.14480655
## [59,] -0.14374547 0.84916211 -0.18511051
## [60,] -2.01801101 -0.80478499 -0.16314586
## [61,] 9.08958595 -0.28383349 1.21798324
## [62,] 1.90802322 -1.98730629 -1.10223094
## [63,] -4.23384159 -5.54566836 -1.09613641
## [64,] -13.50058382 -1.34068536 0.74429666
## [65,] -1.51474474 -3.47804891 -1.52249964
## [66,] -6.41256561 -0.82677054 -1.49034898
## [67,] -3.82589318 -3.20528016 -1.19033778
## [68,] -7.89969963 -6.41339000 -3.12872415
## [69,] -9.71073985 -4.11174278 -1.27199025
## [70,] -12.96556786 -5.48835886 -0.16684863
## [71,] -8.75734121 -2.03643009 -1.28379818
## [72,] 12.45930755 2.49895022 3.12997556
## [73,] 26.60254256 3.11501732 0.80144362
## [74,] 42.59822839 5.24272446 1.72469472
## [75,] -32.99032923 -4.21343249 -3.47636928
## [76,] -5.37710031 0.07492994 -0.64461902
## [77,] 0.54551537 1.16873169 0.51519905
## [78,] -2.94087815 0.49356427 1.70184539
## [79,] -6.94265936 0.57244893 1.30659251
## [80,] -4.73182530 1.44727509 0.01441570
## [81,] 3.79476361 -0.77327555 0.92897463
## [82,] -6.43636087 -1.58927557 -0.83645112
## [83,] -1.21166018 1.50335197 0.74088517
## [84,] 2.12034684 4.48690888 3.71130053
## [85,] 16.42646858 16.41059882 5.85670113
## [86,] -0.84830669 0.02038738 -0.56658934
## [87,] -8.72078653 -2.38969998 -1.37305939
## [88,] 7.71477868 5.99270926 -0.78033737
## [89,] 2.38240802 1.92564332 -2.16957268
## [90,] -40.62264128 -10.03482353 -0.34159378SUR$nlmest## $minimum
## [1] 257.3351
##
## $estimate
## [1] 5.442026e-04 1.619601e+00 1.343648e+00 5.655787e-06 3.336614e+00
## [6] -1.063237e+00 4.402074e-01 2.215587e-04 -3.817718e-03
##
## $gradient
## [1] -1.143399e-02 -2.960711e-05 -1.561974e-05 1.399352e+00 4.365029e-05
## [6] 2.431022e-05 -2.221540e-07 -4.440257e-02 -3.143265e-04
##
## $code
## [1] 1
##
## $iterations
## [1] 1132.3.2. Biomass conversion and expansion factors
#### visualization of Relationship between tree volume and ADB, SDB, BDB, LDB
# subset data for forest stand level
plot <- with(data_plot, data.frame(Plot,Site, Tree.species, Age,BA_m2_ha,
N_trees_ha,Dg_mm,Hg_m,V_m3_ha,SDB_Mg_ha,
BDB_Mg_ha,LDB_Mg_ha,ADB_Mg_ha))
# rename forest stand data column
plot$age <- plot$Age
plot$BA <- plot$BA_m2_ha
plot$Dg <- plot$Dg_mm
plot$Hg <- plot$Hg_m
plot$V <- plot$V_m3_ha
plot$ADB <- plot$ADB_Mg_ha
plot$SDB <- plot$SDB_Mg_ha
plot$BDB <- plot$BDB_Mg_ha
plot$LDB <- plot$LDB_Mg_ha
# data discription
p12 <- ggplot(data = plot, aes(V,ADB))+
geom_point()+ labs(title = "(a)", x ="V (m3)" , y="ADB (Mg_ha)")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()
p13 <- ggplot(data = plot, aes(V,SDB))+
geom_point()+ labs(title = "(b)", x ="V (m3)" , y="SDB (Mg_ha)")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()
p14 <- ggplot(data = plot, aes(V,BDB))+
geom_point()+ labs(title = "(c)", x ="V (m3)" , y="BDB (Mg_ha)")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()
p15 <- ggplot(data = plot, aes(V,LDB))+
geom_point()+ labs(title = "(d)", x ="V (m3)" , y="LDB (Mg_ha)")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()
grid.arrange(p12,p13,p14,p15,ncol=2, top=textGrob('Figure 2. Relationship between tree volume (V) and \n ADB (a), SDB (b), BDB (c), and LDB (d).\n', gp=gpar(fontsize=16)))## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
For stand-level biomass analyses, we calculated BCEFs based on the dry biomass per hectare of particular components (stem, branches, foliage, and aboveground biomass) and stem volume per hectare as BCEF = W/V, where W - dry mass of the considered biomass component [Mg ha−1] and V - total stem volume of trees [m3 ha−1]. The following model types were used to analyze the correlations between BCEFs and tree stand characteristics:
\[ (11)\quad BCEF = a×x^b \\ (12)\quad BCEF = a+b×log(x) \\ (13)\quad BCEF = a+b/x \\ (14)\quad BCEF = e^{-( a+b/x)} \\ (15)\quad BCEF = a + b×e^{-(c/x)} \\ \]
where x - considered tree stand characteristic (age, basal area, volume, quadratic mean diameter, mean height), a, b and c - model coefficients, and e - base of the natural logarithm.
#### nls models for BEF_SDB
plot$BEF_SDB<-plot$SDB/plot$V
summary(plot$BEF_SDB)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4141 0.4569 0.4909 0.4860 0.5147 0.5551### age
BEF_SDB<-plot$BEF_SDB
x<-plot$age
M1<-nls(BEF_SDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_SDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_SDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_SDB~exp(a+b/x),start=list(a=1,b=1))
M5<-nlsLM(BEF_SDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.5171710 0.5107437 0.5398656 0.5393095 0.5495497# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.03047890 0.03068109 0.02975397 0.02977195 0.03040468### BA
BEF_SDB<-plot$BEF_SDB
x<-plot$BA
M1<-nls(BEF_SDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_SDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_SDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_SDB~exp(a+b/x),start=list(a=1,b=1))
# M5<-nlsLM(BEF_SDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodeness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.01485694 0.01388431 0.03319683 0.03529157 0.54954968# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.04353638 0.04355786 0.04312923 0.04308248 0.03040468### Dg
BEF_SDB<-plot$BEF_SDB
x<-plot$Dg
M1<-nls(BEF_SDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_SDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_SDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_SDB~exp(a+b/x),start=list(a=1,b=1))
# M5<-nlsLM(BEF_SDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
# M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.6472831 0.6378393 0.6753937 0.6798988 0.5495497# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.02605048 0.02639692 0.02499085 0.02481683 0.03040468### Hg
BEF_SDB<-plot$BEF_SDB
x<-plot$Hg
M1<-nls(BEF_SDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_SDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_SDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_SDB~exp(a+b/x),start=list(a=1,b=1))
# M5<-nlsLM(BEF_SDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
# M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.7114184 0.7100111 0.6992226 0.6875340 0.5495497# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.02356335 0.02362073 0.02405610 0.02451907 0.03040468### V
BEF_SDB<-plot$BEF_SDB
x<-plot$V
M1<-nls(BEF_SDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_SDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_SDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_SDB~exp(a+b/x),start=list(a=1,b=1))
#M5<-nlsLM(BEF_SDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
#M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_SDB-mean(BEF_SDB)}^2))
## show goodness-of-fit measurement
#R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.3857271 0.3876238 0.3485571 0.3376469 0.5495497# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.03437821 0.03432510 0.03540306 0.03569829 0.03040468#### nls models for BEF_BDB ####
### age
plot$BEF_BDB <- plot$BDB/plot$V
BEF_BDB<-plot$BEF_BDB
x<-plot$age
M1<-nls(BEF_BDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_BDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_BDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_BDB~exp(a+b/x),start=list(a=1,b=1))
M5<-nlsLM(BEF_BDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.6025822 0.5946508 0.6197391 0.6176302 0.6292397# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.003293628 0.003326332 0.003221749 0.003230671 0.003285579### BA
BEF_BDB<-plot$BEF_BDB
x<-plot$BA
M1<-nls(BEF_BDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_BDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_BDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_BDB~exp(a+b/x),start=list(a=1,b=1))
#M5<-nlsLM(BEF_LDB~a+b*exp(-x*c),start=list(a=0.03,b=0.7,c=0.5))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
#M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.02908348 0.02665391 0.05380396 0.05812916 0.62923971# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.005148039 0.005154476 0.005082080 0.005070451 0.003285579### Dg
BEF_BDB<-plot$BEF_BDB
x<-plot$Dg
M1<-nls(BEF_BDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_BDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_BDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_BDB~exp(a+b/x),start=list(a=1,b=1))
# M5<-nlsLM(BEF_BDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
# M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.7382197 0.7239322 0.7699823 0.7774402 0.6292397# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.002673126 0.002745104 0.002505714 0.002464758 0.003285579### Hg
BEF_BDB<-plot$BEF_BDB
x<-plot$Hg
M1<-nls(BEF_BDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_BDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_BDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_BDB~exp(a+b/x),start=list(a=1,b=1))
#M5<-nlsLM(BEF_BDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
#M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.7893533 0.7840869 0.7843557 0.7721956 0.6292397# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.002397884 0.002427674 0.002426162 0.002493630 0.003285579### V
BEF_BDB<-plot$BEF_BDB
x<-plot$V
M1<-nls(BEF_BDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_BDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_BDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_BDB~exp(a+b/x),start=list(a=1,b=1))
#M5<-nlsLM(BEF_BDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
#M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_BDB-mean(BEF_BDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.4519715 0.4534764 0.4169822 0.4023292 0.6292397# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.003867697 0.003862383 0.003989255 0.004039075 0.003285579#### nls models for BEF_LDB
### age
plot$BEF_LDB<-plot$LDB/plot$V
BEF_LDB<-plot$BEF_LDB
x<-plot$age
M1<-nls(BEF_LDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_LDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_LDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_LDB~exp(a+b/x),start=list(a=1,b=1))
M5<-nlsLM(BEF_LDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
## show goodness-of-fit measuremnt
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.7730380 0.7630210 0.7663181 0.7362044 0.7786157# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.006976330 0.007128618 0.007078855 0.007521149 0.007116037### BA
BEF_LDB<-plot$BEF_LDB
x<-plot$BA
M1<-nls(BEF_LDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_LDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_LDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_LDB~exp(a+b/x),start=list(a=1,b=1))
# M5<-nlsLM(BEF_LDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
#M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.12006334 0.08453938 0.13421377 0.17509121 0.77861571# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.013736502 0.014011037 0.013625605 0.013300054 0.007116037### Dg
BEF_LDB<-plot$BEF_LDB
x<-plot$Dg
M1<-nls(BEF_LDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_LDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_LDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_LDB~exp(a+b/x),start=list(a=1,b=1))
#M5<-nls(BEF_LDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
# M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.9366212 0.8682821 0.9346970 0.9728657 0.7786157# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.003686569 0.005314628 0.003742114 0.002412182 0.007116037### Hg
BEF_LDB<-plot$BEF_LDB
x<-plot$Hg
M1<-nls(BEF_LDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_LDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_LDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_LDB~exp(a+b/x),start=list(a=1,b=1))
#M5<-nlsLM(BEF_LDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
# M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.9137509 0.8715123 0.9110343 0.8915781 0.7786157# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.004300587 0.005249056 0.004367791 0.004821797 0.007116037### V
BEF_LDB<-plot$BEF_LDB
x<-plot$V
M1<-nls(BEF_LDB~a*x^b,start=list(a=1,b=1))
M2<-nls(BEF_LDB~a+b*log(x),start=list(a=1,b=1))
M3<-nls(BEF_LDB~a+b/x,start=list(a=1,b=1))
M4<-nls(BEF_LDB~exp(a+b/x),start=list(a=1,b=1))
#M5<-nlsLM(BEF_LDB~a+b*exp(-x*c),start=list(a=1,b=1,c=1))
## goodness-of-fit measurement
M1_r2 <- 1-(sum({summary(M1)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M2_r2 <- 1-(sum({summary(M2)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M3_r2 <- 1-(sum({summary(M3)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
M4_r2 <- 1-(sum({summary(M4)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
#M5_r2 <- 1-(sum({summary(M5)$residuals}^2)/sum({BEF_LDB-mean(BEF_LDB)}^2))
## show goodness-of-fit measurement
# R2
R2<-c(M1_r2,M2_r2,M3_r2,M4_r2,M5_r2)
R2## [1] 0.5852619 0.5891149 0.5729966 0.5242669 0.7786157# RSE
rse<- c(summary(M1)$sigma,summary(M2)$sigma,summary(M3)$sigma,
summary(M4)$sigma,summary(M5)$sigma)
rse## [1] 0.009430559 0.009386651 0.009568991 0.010100250 0.007116037The best models were chosen based on the goodness of fit of the resulting models and their parameters were assessed by examining the pattern of residuals, coefficient of determination (R2), and residual standard error (RSE). Lastly, all chosen models of every biomass component were used as a base for Seemingly Unrelated Regression (SUR). This approach allows fulfillment of the logical assumption that the sum of the estimated biomass components (stem, branches, and foliage) should be equal to the estimated total aboveground biomass, assuring additivity of the biomass equations (Bronisz et al., 2016).
#### Seemingly Unrelated Regression - SUR for BEF
Ms.formula<- BEF_SDB~b0*Hg^b1
Mb.formula<- BEF_BDB~b2*Hg^b3
Ml.formula<- BEF_LDB~exp(b4+b5/Dg)
equations <- list(Ms.formula,Mb.formula,Ml.formula)
start.values <- c(b0=0.99532 ,b1=-0.26835 ,b2=0.119390
,b3=-0.338020,b4=-4.351000,b5=150.518260)
labels <- list("BEF_SDB","BEF_BDB","BEF_LDB")
SUR<-nlsystemfit(method ="SUR",equations,start.values,data=plot,eqnlabels=labels)## Warning in nlm(knls, estols$estimate, typsize = abs(estols$estimate), iterlim =
## maxiter, : NA/Inf replaced by maximum positive valueprint(SUR)## $eq
## $eq[[1]]
## $method
## [1] "SUR"
##
## $i
## [1] 1
##
## $eqnlabel
## [1] "BEF_SDB"
##
## $formula
## BEF_SDB ~ b0 * Hg^b1
##
## $b
## b0 b1
## 1.0374630 -0.2842608
##
## $se
## b0
## [1,] 0.11559913
## [2,] 0.04198083
##
## $t
## b0
## [1,] 8.974661
## [2,] -6.771205
##
## $p
## b0
## [1,] 1.210323e-07
## [2,] 4.490560e-06
##
## $n
## [1] 18
##
## $k
## [1] 2
##
## $df
## [1] 16
##
## $predicted
## [,1]
## [1,] 0.4826666
## [2,] 0.5294458
## [3,] 0.5334670
## [4,] 0.4937191
## [5,] 0.4908735
## [6,] 0.5648337
## [7,] 0.4255401
## [8,] 0.5200902
## [9,] 0.4608357
## [10,] 0.4925923
## [11,] 0.5194407
## [12,] 0.4943387
## [13,] 0.4489660
## [14,] 0.4429839
## [15,] 0.4336085
## [16,] 0.4888855
## [17,] 0.4628156
## [18,] 0.4560482
##
## $residuals
## [,1]
## [1,] 0.007335544
## [2,] -0.010733839
## [3,] -0.007771779
## [4,] -0.001856460
## [5,] 0.001774853
## [6,] -0.016884685
## [7,] 0.053289990
## [8,] -0.004658071
## [9,] 0.021973756
## [10,] 0.013425162
## [11,] 0.035692582
## [12,] 0.018165121
## [13,] -0.019940193
## [14,] 0.006566232
## [15,] -0.019464258
## [16,] 0.001051915
## [17,] -0.037368560
## [18,] -0.033100016
##
## $ssr
## [1] 0.00896236
##
## $mse
## [1] 0.0005601475
##
## $s2
## [1] 0.0005601475
##
## $rmse
## [1] 0.02366744
##
## $s
## [1] 0.02366744
##
## $covb
## b0 b1
## b0 0.013363160 -0.004827241
## b1 -0.004827241 0.001762390
##
## $r2
## [,1]
## [1,] 0.7088632
##
## $adjr2
## [,1]
## [1,] 0.6906671
##
## attr(,"class")
## [1] "nlsystemfit.equation"
##
## $eq[[2]]
## $method
## [1] "SUR"
##
## $i
## [1] 2
##
## $eqnlabel
## [1] "BEF_BDB"
##
## $formula
## BEF_BDB ~ b2 * Hg^b3
##
## $b
## b2 b3
## 0.1248819 -0.3552512
##
## $se
## b2
## [1,] 0.01419200
## [2,] 0.04298461
##
## $t
## b2
## [1,] 8.799458
## [2,] -8.264613
##
## $p
## b2
## [1,] 1.577880e-07
## [2,] 3.626464e-07
##
## $n
## [1] 18
##
## $k
## [1] 2
##
## $df
## [1] 16
##
## $predicted
## [,1]
## [1,] 0.04799324
## [2,] 0.05387502
## [3,] 0.05438687
## [4,] 0.04937060
## [5,] 0.04901523
## [6,] 0.05841223
## [7,] 0.04100257
## [8,] 0.05268789
## [9,] 0.04529591
## [10,] 0.04922982
## [11,] 0.05260568
## [12,] 0.04944803
## [13,] 0.04384259
## [14,] 0.04311375
## [15,] 0.04197644
## [16,] 0.04876728
## [17,] 0.04553925
## [18,] 0.04470859
##
## $residuals
## [,1]
## [1,] 0.0001941516
## [2,] -0.0012926306
## [3,] -0.0005465253
## [4,] -0.0009387548
## [5,] -0.0004490905
## [6,] -0.0014889629
## [7,] 0.0056170080
## [8,] -0.0004623286
## [9,] 0.0018655551
## [10,] 0.0014192195
## [11,] 0.0040934551
## [12,] 0.0017993394
## [13,] -0.0017957228
## [14,] 0.0008388619
## [15,] -0.0016140333
## [16,] 0.0003252876
## [17,] -0.0035809939
## [18,] -0.0032986122
##
## $ssr
## [1] 9.290382e-05
##
## $mse
## [1] 5.806489e-06
##
## $s2
## [1] 5.806489e-06
##
## $rmse
## [1] 0.002409666
##
## $s
## [1] 0.002409666
##
## $covb
## b2 b3
## b2 0.0002014127 -0.0006068065
## b3 -0.0006068065 0.0018476769
##
## $r2
## [,1]
## [1,] 0.7872783
##
## $adjr2
## [,1]
## [1,] 0.7739832
##
## attr(,"class")
## [1] "nlsystemfit.equation"
##
## $eq[[3]]
## $method
## [1] "SUR"
##
## $i
## [1] 3
##
## $eqnlabel
## [1] "BEF_LDB"
##
## $formula
## BEF_LDB ~ exp(b4 + b5/Dg)
##
## $b
## b4 b5
## -4.377619 153.858544
##
## $se
## b4
## [1,] 0.05659659
## [2,] 6.44679920
##
## $t
## b4
## [1,] -77.34775
## [2,] 23.86588
##
## $p
## b4
## [1,] 0.00000e+00
## [2,] 6.17284e-14
##
## $n
## [1] 18
##
## $k
## [1] 2
##
## $df
## [1] 16
##
## $predicted
## [,1]
## [1,] 0.02980005
## [2,] 0.05021175
## [3,] 0.05848215
## [4,] 0.03009420
## [5,] 0.03120386
## [6,] 0.06566257
## [7,] 0.02210475
## [8,] 0.05150518
## [9,] 0.02443925
## [10,] 0.04477699
## [11,] 0.06235409
## [12,] 0.04431272
## [13,] 0.02752593
## [14,] 0.02526685
## [15,] 0.02210475
## [16,] 0.04574420
## [17,] 0.03190019
## [18,] 0.02621420
##
## $residuals
## [,1]
## [1,] 0.0017449830
## [2,] 0.0010778818
## [3,] -0.0002383775
## [4,] 0.0022707302
## [5,] 0.0018611220
## [6,] 0.0041992195
## [7,] 0.0033788285
## [8,] -0.0008169467
## [9,] 0.0030148899
## [10,] -0.0021042007
## [11,] -0.0026271974
## [12,] -0.0016725723
## [13,] -0.0015466847
## [14,] 0.0007510863
## [15,] 0.0004000234
## [16,] -0.0038235279
## [17,] -0.0031527518
## [18,] -0.0010726256
##
## $ssr
## [1] 9.464394e-05
##
## $mse
## [1] 5.915246e-06
##
## $s2
## [1] 5.915246e-06
##
## $rmse
## [1] 0.002432128
##
## $s
## [1] 0.002432128
##
## $covb
## b4 b5
## b4 0.003203175 -0.3536427
## b5 -0.353642712 41.5612199
##
## $r2
## [,1]
## [1,] 0.9724151
##
## $adjr2
## [,1]
## [1,] 0.970691
##
## attr(,"class")
## [1] "nlsystemfit.equation"
##
##
## $resids
## [,1] [,2] [,3]
## [1,] 0.007335544 0.0001941516 0.0017449830
## [2,] -0.010733839 -0.0012926306 0.0010778818
## [3,] -0.007771779 -0.0005465253 -0.0002383775
## [4,] -0.001856460 -0.0009387548 0.0022707302
## [5,] 0.001774853 -0.0004490905 0.0018611220
## [6,] -0.016884685 -0.0014889629 0.0041992195
## [7,] 0.053289990 0.0056170080 0.0033788285
## [8,] -0.004658071 -0.0004623286 -0.0008169467
## [9,] 0.021973756 0.0018655551 0.0030148899
## [10,] 0.013425162 0.0014192195 -0.0021042007
## [11,] 0.035692582 0.0040934551 -0.0026271974
## [12,] 0.018165121 0.0017993394 -0.0016725723
## [13,] -0.019940193 -0.0017957228 -0.0015466847
## [14,] 0.006566232 0.0008388619 0.0007510863
## [15,] -0.019464258 -0.0016140333 0.0004000234
## [16,] 0.001051915 0.0003252876 -0.0038235279
## [17,] -0.037368560 -0.0035809939 -0.0031527518
## [18,] -0.033100016 -0.0032986122 -0.0010726256
##
## $method
## [1] "SUR"
##
## $n
## [1] 54
##
## $k
## [1] 6
##
## $b
## b0 b1 b2 b3 b4 b5
## 1.0374630 -0.2842608 0.1248819 -0.3552512 -4.3776192 153.8585443
##
## $se
## b0 b1 b2 b3 b4 b5
## 0.11559913 0.04198083 0.01419200 0.04298461 0.05659659 6.44679920
##
## $t
## b0 b1 b2 b3 b4 b5
## 8.974661 -6.771205 8.799458 -8.264613 -77.347749 23.865881
##
## $p
## b0 b1 b2 b3 b4 b5
## 7.727374e-12 1.649620e-08 1.401812e-11 8.795364e-11 0.000000e+00 0.000000e+00
##
## $solvtol
## [1] 2.220446e-16
##
## $covb
## b0 b1 b2 b3 b4
## b0 0.0133631596 -0.0048272413 1.621729e-03 -0.0048761609 -9.567182e-04
## b1 -0.0048272413 0.0017623900 -5.869877e-04 0.0017837858 3.773597e-04
## b2 0.0016217285 -0.0005869877 2.014127e-04 -0.0006068065 -7.800109e-05
## b3 -0.0048761609 0.0017837858 -6.068065e-04 0.0018476769 2.567387e-04
## b4 -0.0009567182 0.0003773597 -7.800109e-05 0.0002567387 3.203175e-03
## b5 0.1266165798 -0.0480023612 1.027811e-02 -0.0324795658 -3.536427e-01
## b5
## b0 0.12661658
## b1 -0.04800236
## b2 0.01027811
## b3 -0.03247957
## b4 -0.35364271
## b5 41.56121991
##
## $rcov
## [,1] [,2] [,3]
## [1,] 5.601475e-04 5.637678e-05 1.289121e-05
## [2,] 5.637678e-05 5.806489e-06 9.020586e-07
## [3,] 1.289121e-05 9.020586e-07 5.915246e-06
##
## $rcor
## [,1] [,2] [,3]
## [1,] 1.0000000 0.9885338 0.2234253
## [2,] 0.9885338 1.0000000 0.1533946
## [3,] 0.2234253 0.1533946 1.0000000
##
## $drcov
## [1] 3.290353e-16
##
## $rcovest
## [,1] [,2] [,3]
## [1,] 5.552313e-04 5.584097e-05 1.408354e-05
## [2,] 5.584097e-05 5.749850e-06 1.076514e-06
## [3,] 1.408354e-05 1.076514e-06 5.818621e-06
##
## $g
## [1] 3
##
## $nlmest
## $nlmest$minimum
## [1] 47.5159
##
## $nlmest$estimate
## [1] 1.0374630 -0.2842608 0.1248819 -0.3552512 -4.3776192 153.8585443
##
## $nlmest$gradient
## [1] -8.766673e-04 -2.224509e-03 6.524711e-03 2.018723e-03 6.153968e-05
## [6] 6.361916e-07
##
## $nlmest$code
## [1] 2
##
## $nlmest$iterations
## [1] 46
##
##
## attr(,"class")
## [1] "nlsystemfit.system"summary(SUR)## Length Class Mode
## eq 3 -none- list
## resids 54 -none- numeric
## method 1 -none- character
## n 1 -none- numeric
## k 1 -none- numeric
## b 6 -none- numeric
## se 6 -none- numeric
## t 6 -none- numeric
## p 6 -none- numeric
## solvtol 1 -none- numeric
## covb 36 -none- numeric
## rcov 9 -none- numeric
## rcor 9 -none- numeric
## drcov 1 -none- numeric
## rcovest 9 -none- numeric
## g 1 -none- numeric
## nlmest 5 -none- listSUR$eq[3]## [[1]]
## $method
## [1] "SUR"
##
## $i
## [1] 3
##
## $eqnlabel
## [1] "BEF_LDB"
##
## $formula
## BEF_LDB ~ exp(b4 + b5/Dg)
##
## $b
## b4 b5
## -4.377619 153.858544
##
## $se
## b4
## [1,] 0.05659659
## [2,] 6.44679920
##
## $t
## b4
## [1,] -77.34775
## [2,] 23.86588
##
## $p
## b4
## [1,] 0.00000e+00
## [2,] 6.17284e-14
##
## $n
## [1] 18
##
## $k
## [1] 2
##
## $df
## [1] 16
##
## $predicted
## [,1]
## [1,] 0.02980005
## [2,] 0.05021175
## [3,] 0.05848215
## [4,] 0.03009420
## [5,] 0.03120386
## [6,] 0.06566257
## [7,] 0.02210475
## [8,] 0.05150518
## [9,] 0.02443925
## [10,] 0.04477699
## [11,] 0.06235409
## [12,] 0.04431272
## [13,] 0.02752593
## [14,] 0.02526685
## [15,] 0.02210475
## [16,] 0.04574420
## [17,] 0.03190019
## [18,] 0.02621420
##
## $residuals
## [,1]
## [1,] 0.0017449830
## [2,] 0.0010778818
## [3,] -0.0002383775
## [4,] 0.0022707302
## [5,] 0.0018611220
## [6,] 0.0041992195
## [7,] 0.0033788285
## [8,] -0.0008169467
## [9,] 0.0030148899
## [10,] -0.0021042007
## [11,] -0.0026271974
## [12,] -0.0016725723
## [13,] -0.0015466847
## [14,] 0.0007510863
## [15,] 0.0004000234
## [16,] -0.0038235279
## [17,] -0.0031527518
## [18,] -0.0010726256
##
## $ssr
## [1] 9.464394e-05
##
## $mse
## [1] 5.915246e-06
##
## $s2
## [1] 5.915246e-06
##
## $rmse
## [1] 0.002432128
##
## $s
## [1] 0.002432128
##
## $covb
## b4 b5
## b4 0.003203175 -0.3536427
## b5 -0.353642712 41.5612199
##
## $r2
## [,1]
## [1,] 0.9724151
##
## $adjr2
## [,1]
## [1,] 0.970691
##
## attr(,"class")
## [1] "nlsystemfit.equation"SUR$resids## [,1] [,2] [,3]
## [1,] 0.007335544 0.0001941516 0.0017449830
## [2,] -0.010733839 -0.0012926306 0.0010778818
## [3,] -0.007771779 -0.0005465253 -0.0002383775
## [4,] -0.001856460 -0.0009387548 0.0022707302
## [5,] 0.001774853 -0.0004490905 0.0018611220
## [6,] -0.016884685 -0.0014889629 0.0041992195
## [7,] 0.053289990 0.0056170080 0.0033788285
## [8,] -0.004658071 -0.0004623286 -0.0008169467
## [9,] 0.021973756 0.0018655551 0.0030148899
## [10,] 0.013425162 0.0014192195 -0.0021042007
## [11,] 0.035692582 0.0040934551 -0.0026271974
## [12,] 0.018165121 0.0017993394 -0.0016725723
## [13,] -0.019940193 -0.0017957228 -0.0015466847
## [14,] 0.006566232 0.0008388619 0.0007510863
## [15,] -0.019464258 -0.0016140333 0.0004000234
## [16,] 0.001051915 0.0003252876 -0.0038235279
## [17,] -0.037368560 -0.0035809939 -0.0031527518
## [18,] -0.033100016 -0.0032986122 -0.0010726256SUR$nlmest## $minimum
## [1] 47.5159
##
## $estimate
## [1] 1.0374630 -0.2842608 0.1248819 -0.3552512 -4.3776192 153.8585443
##
## $gradient
## [1] -8.766673e-04 -2.224509e-03 6.524711e-03 2.018723e-03 6.153968e-05
## [6] 6.361916e-07
##
## $code
## [1] 2
##
## $iterations
## [1] 46
3. Results and Discussion
3.1. Allometric equations
All the allometric equations developed for each biomass component explained from 69.5 to 98.4% of the variation in biomass. For SDB, BDB and LDB, R2 varied from 0.801 to 0.984, from 0.695 to 0.795, and from 0.741 to 0.789, respectively
# R2 for SDB estimation
Mssrr2## [1] 0.8932732 0.9592599 0.8016335 0.8997697 0.9680587 0.9822738 0.9833863
## [8] 0.9672339 0.9841831 0.9833288# R2 for BDB estimation
Mbsrr2## [1] 0.6955300 0.7593847 0.4813710 0.7224357 0.7920236 0.7337304 0.7909975
## [8] 0.7951524 0.7336445# R2 for LDB estimation
Mlsrr2## [1] 0.7444385 0.7780796 0.5550345 0.7550876 0.7885359 0.7410768 0.7880249
## [8] 0.7774793 0.7882190 0.75129583.1.1. Model fitting
Based on the goodness-of-fit measures and model parameters’ accuracy, the most appropriate equations to estimate stem, branches, and foliage biomass were chosen separately (Table 3). Models number 9 and 5 showed the high goodness-of-fit in most cases. Notably, in terms of coefficient of determination, SDB was predicted much more accurately compared to BDB and LDB, with the R2 value is 0.984, 0.795, 0.789, respectively. In contrast, predicted models for BDB and LDB show better RSE and AIC values than the SDB. In general, biomass components were predicted by both DBH and H present better accuracy than models with exclusively DBH or H as independent variables.
# best model (9) for SDB estimation
summary(Ms9)##
## Formula: SDB ~ b0 * D^b1 * H^b2
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b0 0.0005357 0.0001209 4.431 2.72e-05 ***
## b1 1.6941254 0.0913716 18.541 < 2e-16 ***
## b2 1.2137352 0.1303440 9.312 1.03e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.17 on 87 degrees of freedom
##
## Algorithm "port", convergence message: relative convergence (4)# best model (9) for BDB estimation
summary(Mb9)##
## Formula: BDB ~ b0 * D^b1 * H^b2
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b0 5.215e-06 5.512e-06 0.946 0.34666
## b1 3.609e+00 3.948e-01 9.142 2.3e-14 ***
## b2 -1.532e+00 4.930e-01 -3.108 0.00254 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.403 on 87 degrees of freedom
##
## Algorithm "port", convergence message: relative convergence (4)# best model (5) for LDB estimation
summary(Ml5)##
## Formula: LDB ~ b0 + b1 * D2 + b2 * H2
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b0 1.224e+00 6.071e-01 2.017 0.0468 *
## b1 2.484e-04 2.751e-05 9.033 3.85e-14 ***
## b2 -9.576e-03 4.617e-03 -2.074 0.0410 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.219 on 87 degrees of freedom
##
## Number of iterations to convergence: 1
## Achieved convergence tolerance: 6.089e-06Table 3. Goodness-of-fit measures of the best models for SDB, BDB, and LDB.
| Biomass Component | Col2 | Col3 | Col4 |
|---|---|---|---|
| Equation No. | 9 | 9 | 5 |
| a | 0.001 | 0.00001 | 1.224 |
| SE (a) | 0.0001 | 0.00001 | 0.6071 |
| b | 1.694 | 3.609 | 0.00025 |
| SE (b) | 0.091 | 0.395 | 0.00003 |
| c | 1.214 | -1.532 | -0.010 |
| SE (c) | 0.130 | 0.493 | 0.005 |
| R2 | 0.984 | 0.795 | 0.789 |
| RSE | 12.174 | 4.403 | 2.219 |
| AIC | 710.2 | 527.2 | 403.9 |
3.1.2. Seemingly Unrelated Regression
Models for different small biomass parts were fitted. Because total aboveground biomass is explicitly related to the biomass of tree components as their sum, it was not included in the set of equations. Estimates for all biomass components were obtained as a function of DBH and height (Table 4). Residual distribution for analyzed aboveground biomass models and the relationship between predicted and observed ADB are shown in Figure 3.
Table 4. Goodness-of-fit for the final models for biomass components. Parameters (a, b, c) with their standard errors (SE) and R2 = coefficient of determination, RSE = residual standard error.
| Biomass Component | Equation No. | a | SE (a) | b | SE (b) | c | SE (c) | R2 | RMSE |
|---|---|---|---|---|---|---|---|---|---|
| SDB | 9 | 0.0005 | 0.0001 | 1.620 | 0.084 | 1.344 | 0.118 | 0.984 | 12.25 |
| BDB | 9 | 0.000006 | 0.00005 | 3.337 | 0.330 | -1.063 | 0.400 | 0.793 | 4.43 |
| LDB | 5 | 0.440 | 0.435 | 0.000 | 0.000 | -0.004 | 0.004 | 0.784 | 2.24 |
#### visualization of model residuals and model fitting
# load prediction result
ADB_predicted <- read_excel("data/ADB_predicted.xlsx")
# residuals plot
p10 <- ggplot(data = ADB_predicted, aes(ADB_predicted$pridicted, ADB_predicted$residuals))+
geom_point()+theme_bw()+labs(title = "(a)", x="Predicted ADB (kg)", y="Residuals")
## ADB measured and predicted plot
# load prediction and reference data
ADB_compare <- read_excel("data/ADB_compare.xlsx")
p11 <- ggplot(data = ADB_compare, aes(ADB_compare$Measured,ADB_compare$predicted))+
geom_point()+labs(title = "(b)",x="Observed ADB (kg)", y="Predicted ADB (kg)")+
geom_smooth(method='lm',se = FALSE, color="black")+
# stat_poly_eq(formula = y ~ x,
# aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
# parse = TRUE)+
theme_bw()
grid.arrange(p10,p11,ncol=2, top='Figure 3. Distribution of residuals for analyzed aboveground biomass models (a) \n and the relationship between predicted and observed aboveground biomass (ADB).\n')## Warning: Use of `ADB_predicted$pridicted` is discouraged.
## ℹ Use `pridicted` instead.## Warning: Use of `ADB_predicted$residuals` is discouraged.
## ℹ Use `residuals` instead.## Warning: Use of `ADB_compare$Measured` is discouraged.
## ℹ Use `Measured` instead.## Warning: Use of `ADB_compare$predicted` is discouraged.
## ℹ Use `predicted` instead.## Warning: Use of `ADB_compare$Measured` is discouraged.
## ℹ Use `Measured` instead.## Warning: Use of `ADB_compare$predicted` is discouraged.
## ℹ Use `predicted` instead.## `geom_smooth()` using formula = 'y ~ x'
3.2. Biomass conversion and expansion factor
3.2.1. Model fitting
At the tree stand level, Scots pine biomass components were strongly correlated with mean H and DBH (Figure. 4, Table. 5). The correlation with age, basal area, and volume were weaker. BCEFs of young Scots pine tree stands differed across components and tree stand characteristics.
#### visualization of relationships between biomass component ~ stand characteristics
### SDB ~ stand characteristics
p16 <- ggplot(data = plot, aes(BEF_SDB,age))+
geom_point()+ labs(title = "SDB",x="",y="")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
p17 <- ggplot(data = plot, aes(BEF_SDB,BA))+
geom_point()+ labs(x="",y="")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
p18 <- ggplot(data = plot, aes(BEF_SDB,Dg))+
geom_point()+ labs(x="",y="")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
p19 <- ggplot(data = plot, aes(BEF_SDB, Hg))+
geom_point()+ labs(x="",y="")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
p20 <- ggplot(data = plot, aes(BEF_SDB, V))+
geom_point()+ labs(x="",y="")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()
### BDB ~ stand characteristics
p21 <- ggplot(data = plot, aes(BEF_BDB,age))+
geom_point()+ labs(title = "BDB",x="",y="")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),)
p22 <- ggplot(data = plot, aes(BEF_BDB,BA))+
geom_point()+ labs(x="",y="")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),)
p23 <- ggplot(data = plot, aes(BEF_BDB,Dg))+
geom_point()+ labs(x="",y="")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
p24 <- ggplot(data = plot, aes(BEF_BDB, Hg))+
geom_point()+ labs(x="",y="")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
p25 <- ggplot(data = plot, aes(BEF_BDB, V))+
geom_point()+ labs(x="",y="")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(axis.text.y=element_blank(),
axis.ticks.y=element_blank())
### LDB ~ stand characteristics
p26 <- ggplot(data = plot, aes(BEF_LDB,age))+
geom_point()+ labs(title = "LDB",x="",y="Age")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())+
scale_y_continuous(position = "right")
p27 <- ggplot(data = plot, aes(BEF_LDB,BA))+
geom_point()+ labs(x="",y="BA")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())+
scale_y_continuous(position = "right")
p28 <- ggplot(data = plot, aes(BEF_LDB,Dg))+
geom_point()+ labs(x="",y="Dg")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())+
scale_y_continuous(position = "right")
p29 <- ggplot(data = plot, aes(BEF_LDB, Hg))+
geom_point()+ labs(x="",y="Hg")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())+
scale_y_continuous(position = "right")
p30 <- ggplot(data = plot, aes(BEF_LDB, V))+
geom_point()+ labs(x="",y="V")+
geom_smooth(method = "lm",se = FALSE, color = "black")+
theme_bw()+
theme(axis.text.y=element_blank(),
axis.ticks.y=element_blank())+
scale_y_continuous(position = "right")
grid.arrange(p16,p21,p26,p17,p22,p27,p18,p23,p28,p19,p24,p29,p20,p25,p30,
ncol=3, left = textGrob("BEF",rot = 90, gp=gpar(fontsize=16)),
top=textGrob('Figure 4. Relationships between tree stand characteristics and BCEFs \n for biomass components: branches (BDB), foliage (LDB), and stem (SDB). \n',gp=gpar(fontsize=16)))## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Based on the goodness-of-fit measures and model parameters’ accuracy, the most appropriate equations to estimate BCEFs of stem, branches, and foliage were chosen separately (Table 5). Models number 4 and 1 showed the high goodness-of-fit in most cases. Notably, in terms of coefficient of determination, BCEFs of LDB were predicted much more accurately compared to BDB and SDB, with the R2 value is 0.973 and 0.789, 0.711, respectively.
Table 5. Goodness-of-fit of the best models for estimating BCEFs of SDB, BDB, and LDB.
| BEF | Predictor | Equation No. | a | SE (a) | b | SE (b) | R2 | RSE |
|---|---|---|---|---|---|---|---|---|
| SDB | Hg | 1 | 0.995 | 0.113 | -0.268 | 0.043 | 0.711 | 0.024 |
| BDB | Hg | 1 | 0.119 | 0.014 | -0.338 | 0.043 | 0.789 | 0.002 |
| LDB | Dg | 4 | -4.351 | 0.057 | -150.5 | 6.552 | 0.973 | 0.002 |
3.2.2. Seemingly Unrelated Regression
Models for BCEFs of different small biomass parts were fitted. The best-fitted models of different biomass parts were used as a base for Seemingly Unrelated Regression (SUR) (Table 6). This approach enabled fulfillment of the logical assumption that the sum of the estimated biomass values of tree parts matches the estimated total biomass, assuring additivity of the biomass equations.
Table 6. Goodness-of-fit of the final models for biomass components. Parameters (a, b, c) with their standard errors (SE) and R2 = coefficient of determination, RSE = residual standard error.
| BEF | Predictor | Equation No. | a | SE (a) | b | SE (b) | R2 | RSE |
|---|---|---|---|---|---|---|---|---|
| SDB | Hg | 1 | 1.037 | 0.116 | -0.284 | 0.042 | 0.709 | 0.0237 |
| BDB | Hg | 1 | 0.125 | 0.014 | -0.355 | 0.043 | 0.787 | 0.0024 |
| LDB | Dg | 4 | -4.378 | 0.057 | -153.859 | 6.447 | 0.972 | 0.0024 |
Residual distributions for analyzed BCEFs of biomass components are shown in Figure 5. The relationship between predicted and observed ADB (kg) using calculated BCEFs is shown in Figure 6. Figure 5 and Figure 6 show that predicted ADB calculated by BCEFs is relatively accurate, with the coefficient of determination is equal to 0.976 explaining almost entirely variation of measured ADB.
# load data
BEF_predicted <- read_excel("data/BEF_predicted.xlsx")
# residuals plot
p31 <- ggplot(data = BEF_predicted, aes(BEF_predicted$SDB_pre, BEF_predicted$SDB_re))+
geom_point()+theme_bw()+labs(title = "SDB", x="Predicted BEF", y="Residuals")+
theme(plot.title = element_text(hjust = 0.5))
p32 <- ggplot(data = BEF_predicted, aes(BEF_predicted$BDB_pre, BEF_predicted$BDB_re))+
geom_point()+theme_bw()+labs(title = "BDB", x="Predicted BEF", y="Residuals")+
theme(plot.title = element_text(hjust = 0.5))
p33 <- ggplot(data = BEF_predicted, aes(BEF_predicted$LDB_pre, BEF_predicted$LDB_re))+
geom_point()+theme_bw()+labs(title = "LDB", x="Predicted BEF", y="Residuals")+
theme(plot.title = element_text(hjust = 0.5))
p34 <- ggplot(data = BEF_predicted, aes(BEF_predicted$ADB_pre, BEF_predicted$ADB_re))+
geom_point()+theme_bw()+labs(title = "ADB", x="Predicted BEF", y="Residuals")+
theme(plot.title = element_text(hjust = 0.5))
grid.arrange(p31,p32,p33,p34,ncol=2, top=textGrob('Figure 5. Distribution of residuals for analyzed BCEFs of biomass components.\n', gp=gpar(fontsize=16)))## Warning: Use of `BEF_predicted$SDB_pre` is discouraged.
## ℹ Use `SDB_pre` instead.## Warning: Use of `BEF_predicted$SDB_re` is discouraged.
## ℹ Use `SDB_re` instead.## Warning: Use of `BEF_predicted$BDB_pre` is discouraged.
## ℹ Use `BDB_pre` instead.## Warning: Use of `BEF_predicted$BDB_re` is discouraged.
## ℹ Use `BDB_re` instead.## Warning: Use of `BEF_predicted$LDB_pre` is discouraged.
## ℹ Use `LDB_pre` instead.## Warning: Use of `BEF_predicted$LDB_re` is discouraged.
## ℹ Use `LDB_re` instead.## Warning: Use of `BEF_predicted$ADB_pre` is discouraged.
## ℹ Use `ADB_pre` instead.## Warning: Use of `BEF_predicted$ADB_re` is discouraged.
## ℹ Use `ADB_re` instead.
# ADB measured and predicted plot
BEF_ADB_compare <- read_excel("data/BEF_ADB_compare.xlsx")
ggplot(data = BEF_ADB_compare, aes(ADB,ADB_predicted))+
geom_point()+labs(x="Observed ADB (kg)", y="Predicted ADB (kg)")+
geom_smooth(method='lm',se = FALSE, color="black")+
theme_bw()+
labs(title = 'Figure 6. Relationship between predicted and observed ADB (kg) \n using calculated BCEFs.\n')## `geom_smooth()` using formula = 'y ~ x'
# stat_poly_eq(formula = y ~ x,
# aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
# parse = TRUE)+
4. Conclusion
Our study provided a comprehensive overview of the tree- and stand-level assessments of aboveground biomass estimation for Scots pine in the study area.
Elaborated empirical equations are built to estimate the value of the biomass of the tree’s above-ground component for Scots pine. Biomass models based on DBH and H are better than models based exclusively on DBH or H. The equations presented at this stage should only be used locally for the tree stands of the research area due to the provenance of the material for the investigations. The produced results can have considerable systematic errors if they are used outside of the study area, although the SUR approach has a relatively wide variety of applications. The accuracy of the produced formula must be assessed on an independent empirical material derived from various habitats and regions of the pine range before they can be used on a larger scale.
Stand-level equations based on BCEFs are provided based on the relationships between tree characteristics and biomass components. The calculated BCEFs may be a useful tool for assessing the biomass and carbon sequestration of the Scots pine stands. However, due to BCEFs are site-specific factors, they should only be applied in similar conditions. In this study, we suggest Dg and Hg might be the best predictors for the stands.
References
Bronisz, K., Bijak, S., Wojtan, R., Tomusiak, R., Bronisz, A., Baran, P., & Zasada, M. (2021). Seemingly Unrelated Mixed-Effects Biomass Models for Black Locust in West Poland. Forests, 12(3), 380. https://doi.org/10.3390/f12030380
Bronisz, K., Strub, M., Cieszewski, C., Bijak, S., Bronisz, A., Tomusiak, R., Wojtan, R., & Zasada, M. (2016). Empirical equations for estimating aboveground biomass of Betula pendula growing on former farmland in central Poland. Silva Fennica, 50(4). https://doi.org/10.14214/sf.1559
Ellenberg, H. (1988). Vegetation ecology of Central Europe (4th ed.). Cambridge University Press. Forest Data Bank. (14-Jun-21). https://www.bdl.lasy.gov.pl/portal/
Jagodziński, A., Dyderski, M., Gęsikiewicz, K., & Horodecki, P. (2018). Tree- and Stand-Level Biomass Estimation in a Larix decidua Mill. Chronosequence. Forests, 9(10), 587. https://doi.org/10.3390/f9100587
Jagodziński, A. M., Dyderski, M. K., Gęsikiewicz, K., Horodecki, P., Cysewska, A., Wierczyńska, S., & Maciejczyk, K. (2018). How do tree stand parameters affect young Scots pine biomass? – Allometric equations and biomass conversion and expansion factors. Forest Ecology and Management, 409, 74–83. https://doi.org/10.1016/j.foreco.2017.11.001
Neumann, M., Moreno, A., Mues, V., Härkönen, S., Mura, M., Bouriaud, O., Lang, M., Achten, W. M., Thivolle-Cazat, A., Bronisz, K., Merganič, J., Decuyper, M., Alberdi, I., Astrup, R., Mohren, F., & Hasenauer, H. (2016). Comparison of carbon estimation methods for European forests. Forest Ecology and Management, 361, 397–420. https://doi.org/10.1016/j.foreco.2015.11.016
Oleksyn, J., Reich, P. B., Chalupka, W., & Tjoelker, M. G. (1999). Differential Above- and Below-ground Biomass Accumulation of European Pinus sylvestris Populations in a 12-year-old Provenance Experiment. Scandinavian Journal of Forest Research, 14(1), 7–17. https://doi.org/10.1080/02827589908540804
Socha J, Wężyk P. (2004). EMPIRICAL FORMULAE TO ASSESS THE BIOMASS OF THE ABOVE-GROUND PART OF PINE TREES. http://www.ejpau.media.pl/volume7/issue2/forestry/art-04.html
Appendix I
Table 7. Scots Pine tree level data
data_tree
Appendix II
Table 8. Scots Pine plot level data
data_plot