Introduction Top

Litsea glutinosa (Lour.) C. B. Rob. belongs to the Lauraceae family and is a small to medium-sized tree species that grows naturally throughout Asia (Heuze et al., 2015; Hinsinger & Strijk, 2016; Useful Tropical Plants, 2021). In the Central Highlands of Vietnam, this species is distributed in tropical evergreen broadleaf forests and is widely planted in combination with cassava (Manihot esculenta Crantz) in the agroforestry model on fallow land of ethnic minorities (Huy, 2009a, 2009b, 2014). In this agroforestry model, L. glutinosa provides valuable biomass to the industry, especially valuable stem bark, and contributes to ecosystem services such as sequestering CO² (Huy, 2009a, 2009b, 2014).

L. glutinosa is a fast-growing multi-purpose plant. Its bark contains essential oils used in the pharmaceutical industry (Heuze et al., 2015), industrial glues, and manufacturing paints (Tiwari et al., 2010). Bark extract is used as an ingredient in commercial cosmetic preparations such as a skin conditioner (Useful Tropical Plants, 2011) and incense for religious services. Its golden brown hardwood is used for furniture and paper manufacturing. Additionally, wood and bark contain gluten and are used as binders (Useful Tropical Plants, 2021). Its leaves are used as fodder for livestock. The medicinal value of L. glutinosa has been mentioned in many studies. Its root bark is anti-inflammatory (Wua et al., 2017), stem bark is an antidiarrheal (Sumithregowda et al., 2017), and L. glutinosa is also used to treat joint and back pain (Pandey & Mandal, 2012). Among the different components, the biomass of the stem bark of L. glutinosa has the highest value (Huy, 2009a, 2009b; Mulia & Nguyen, 2021). Therefore, L. glutinosa is also considered one of the most important non-timber forest tree species in many countries such as Vietnam and India (Mohammad et al., 2020).

The agroforestry model with L. glutinosa and cassava is easy to grow, does not require intensive care, and provides various products from forest trees and agricultural crops. Therefore, it is suitable for the farming practices of ethnic minorities in tropical highlands. On the other hand, planting this tree species in agroforestry models also plays an important role in absorbing CO² to mitigate the greenhouse effect (Mulia & Nguyen, 2021). Therefore, there is a demand to develop a modeling system for L. glutinosa species to estimate individual tree biomass and the carbon sequestration potential of this species in an agroforestry system.

Allometric equations have usually been developed to estimate the biomass of aboveground components, including stem, branches, leaves, bark, and the total aboveground biomass (AGB). Although biomass components are sometimes modeled independently of each other, the tree parts and the tree total biomass are biologically related (Huy et al., 2019). Therefore, component models are best fitted simultaneously with the model for AGB as a system of equations using Seemingly Unrelated Regression (SUR) (Parresol, 2001). SUR accounts for the cross-equation correlation and helps reduce the variability of the parameters and increase the reliability of the biomass estimation of tree parts and totals (Poudel & Temesgen, 2016; Kralicek et al., 2017; Huy et al., 2019; Trautenmüller et al., 2021). The relationship between biomass components and common dendrometric variables (diameter at breast height (D), and tree height (H)) is nonlinear and the residuals are heteroscedastic. Therefore, applying the weighted nonlinear seemingly unrelated regression (WNSUR) method can further improve the reliability of the estimates of AGB and its components (Trautenmüller et al., 2021).

The objectives of this study were to develop and cross-validate the simultaneous modeling systems for estimating components and total tree biomass and carbon of L. glutinosa in an agroforestry model with cassava. We hypothesized that simultaneous estimation of these attributes using WNSUR provides improved and reliable biomass and carbon estimates compared to independently fitted models. We used the dataset derived from a technical report (Huy, 2009a) and applied it to a new methodology to improve biomass and carbon estimates for this studied species in the agroforestry model

Material and methods Top

Study sites and agroforestry model studied

The study area is in the Central Highlands of Vietnam (Fig. 1). The average annual temperature is 21.60C and the average annual rainfall is 2,213 mm (Hydro-meteorological Station in the Central Highlands). Located at an altitude of 400-800 m above sea level, the soil types in the study area include red-brown soil on basalt, gray soil, and red-yellow soil on granite.

Figure 1. Location of the study in the tropical Central Highlands of Vietnam.

The studied agroforestry model includes the native multi-purpose tree species L. glutinosa combined with cassava. L. glutinosa is combined in different proportions with cassava ranging from 50% to 80% of the agroforestry area and is managed in a 6-7 year rotation cycle. Planting density (Nplant) ranged from 500 to 1967 plants ha^-1. L. glutinosa was grown from seeds collected at local tropical evergreen broadleaf forests in the first rotation, and in the following cycles, regenerated stems from coppice were used. The average number of stems plant^-1 (Nstemplant) in the second and third cycles ranged from 1 to 3.4. The number of stems ha^-1 (Nstem) ranged from 500 to 5900. After the third cycle, L. glutinosa was replanted from seeds (Table 1).

Sampling design, data collection, and variable calculation

Twenty-two 300-m² circular plots were sampled covering a full range of ages (A) (1-7 years), a range of planting density, and a range of rotation cycles 1-3 in the agroforestry model. For all trees in the sample plots D (cm), H (m), and the Nstemplant were recorded. In each sample plot, one selected tree having the same D as the quadratic mean diameter (Dg), was destructively sampled to obtain a tree dry biomass/carbon dataset and its components. Age distributions of all sampled plots and destructively sampled trees are shown in Fig. 2. Before felling the sample tree, D, H, and crown diameter (CD, m) in two cardinal directions, North-South and East-West, and the sampled tree’s age (A, year) were recorded. The crown area of the trees sampled (CA, m²) was calculated by the formula

, where CD (the average crown diameter) was obtained from two cardinal direction measurements. After felling, tree heights were remeasured and fresh weights of the stem, branches, leaves, and bark were recorded. The four sampled tree components including stem, bark, branches and leaves were separated and weighed for fresh biomass in the field using a scale with a precision of 0.01 kg. Stem wood and bark samples were obtained from the base, middle, and top sections of the sampled tree. Sampled branches included small and large branches, and the sampled leaves on these branches. For four sampled tree components, a total of 88 samples were collected for the analysis; 100-300 g of each sample was weighed for fresh biomass on site using an electronic scale with a precision of 0.01 g.

Figure 2. Tree age distributions in the 22 sampled plots (left) and destructively sampled trees (right).

The sample materials were dried at 105°C until a constant weight was attained. This provided the average fresh-to-dry mass ratio for all tree components to calculate the respective dry biomass – stem (Bst, kg tree^-1), branches (Bbr, kg tree^-1), leaves (Ble, kg tree^-1), and bark (Bba, kg tree^-1). Total tree aboveground biomass was computed as the sum of component biomass i.e., AGB = Bst + Bbr + Ble + Bba (kg tree^-1). The samples were analyzed after drying and the percentage of carbon was estimated by the Walkley & Black method (1934) and the carbon content of tree stem (Cst, kg tree^-1), branches (Cbr, kg tree^-1), leaves (Cle, kg tree^-1), bark (Cba, kg tree^-1), and the total tree aboveground carbon (AGC = Cst + Cbr + Cle + Cba, kg tree^-1) was calculated. Table 2 shows the summary statistics for each tree predictor and the destructively sampled tree response variables.

Statistical analysis

In this study, we compared two methods to fit systems of equations for AGB/AGC and their components: weighted nonlinear least-squares (WNLS) and weighted nonlinear seemingly unrelated regression (WNSUR) fit by the generalized least squares method.

— Relationship among tree biomass-carbon components and selection of predictors: The factor analysis method was performed to examine the relationship among the tree biomass and carbon components and select their predictors from D, H, A, and CA (Kim & Mueller, 1978; DeCoster, 1998). This method supports both principal components and classical factor analysis, which produces a linear combination of multiple quantitative variables and explains the largest percentage variation among those variables. The values of the variables were standardized by subtracting their means and dividing by their standard deviations.

Figure 3. Scatter plot of biomass and carbon vs. D for tree stem (Bst/Cst, kg), branches (Bbr/Cbr, kg), leaves (Ble/Cle, kg), bark (Bba/Cba, kg), and total tree aboveground biomass/carbon (AGB/AGC, kg).

— Model calibration – Independent fit: Fig. 3 shows that the relationship between tree biomass-carbon components, AGB and AGC vs. D conforms to the power law that was used in the study. Preliminary analysis showed that power models fit by the nonlinear method produced higher reliability than the log-linear model. This is consistent with previous studies (Huy et al., 2016a, 2016b, 2016c). Therefore, we used the nonlinear method to fit biomass-carbon modeling systems, and the heteroscedasticity in residuals was accounted for by using appropriate weighting (Davidian & Giltinan, 1995; Picard et al., 2012; Huy et al., 2016a, 2016b, 2016c, 2019). The WNLS models were fitted using nls function in the statistical software R (R Core Team, 2021). The model forms used were as follows:

where Y_i is the Bst/Cst, Bbr/Cbr, Ble/Cle, Bba/Cba or AGB/AGC in kg for the i^th sampled tree; a and b are the parameters of the model; X_i is the predictor(s) selected by factor analysis such as a combination of D (cm), H (m), A (year), D²H for the i^th sampled tree; and ε_i is the random error associated with the i^th sampled tree. The weighting variable is 1/D^δ or 1/(D²H)^δ and ^δ is selected in a range of ±2 (Picard et al., 2012).

— Model calibration – Simultaneous estimation: Independently fitted component models do not ensure that the AGB/AGC calculated as the sum of predicted tree component models is the same as the AGB/AGC estimated from independently developed AGB/AGC models (Sanquetta et al., 2015; Affleck & Dieguez-Aranda, 2016; Poudel & Temesgen, 2016; Gonzalez-Benecke et al., 2018; Huy et al., 2019; Trautenmüller et al., 2021). The SUR can solve that limitation by allowing simultaneous estimation of the tree biomass or carbon component and AGB/AGC (Parresol, 2001). Additionally, different weighting factors can be used for each equation to account for heteroscedasticity. The WNSUR models were fitted using SAS procedure Proc Model (SAS Inst., 2014) with the generalized least squares method. The modeling systems for the component and total tree biomass and carbon had the following general forms:

where Bst/Cst, Bbr/Cbr, Ble/Cle, Bba/Cba and AGB/AGC are biomass/carbon of stem, branches, leaves, bark, and total in kg, respectively; a_i and b_i are parameters of the power model i(i = 1, 2, 3, 4 for the stem, branches, leaves, and bark respectively); Xij is the predictor variables (D, H, D²HWWD) for the i^th equation and the j^th predictor; and ε_i is the residuals for the i^th equation (i = 1, 2, 3, 4, 5). The weighting variable is 1/D^δ or 1/(D²H)^δ and δ is selected in a range of ±2 (Picard et al., 2012).

Cross-validation

Leave-one-out cross validation (LOOCV) is a special case of K-fold cross-validation, where K is set to the number of samples in the dataset. LOOCV has the maximum computational cost. It requires that one model be created and evaluated for each sample in the dataset (Cheng et al., 2017). LOOCV was applied for the cross-validation of modeling systems. The dataset was split into two parts with n-1 samples for model development and one for cross-validation. The process was repeated for all the observations, and the error statistics were computed at each iteration and averaged over the number of iterations.

The Akaike Information Criterion (AIC) (Akaike, 1973) compares and selects the best model. The adjusted R² describes the variability in the dependent variable explained by the set of predictors accounting for the number of predictors used in the model, whereas the diagnostic plots are used to examine residuals for any departure from model assumptions. The lower the AIC and the closer the adj. R² is to 1, the better the model. At the same time, cross-validation errors based on each sample such as bias (%), root mean square error (RMSE, kg tree^-1), and mean absolute percent error (MAPE, %) were used and averaged over the n realizations of LOOCV. The models with the smallest values of cross-validation errors were preferred.

where n is the number of realizations (number of samples); and y_i and ŷ_i are observed and predicted Bst/Cst, Bbr/Cbr, Ble/Cle, Bba/Cba and AGB/AGC for the i^th realization, respectively.

Bias and MAPE statistics estimated using Eqs. (8) and (10) have been the most widely used by many authors (e.g., Chave et al., 2005; Basuki et al., 2009; Huy et al., 2019, 2022). In some cases, the errors based on Eqs. (8) and (10) tend to be biased on negative errors, y_i < ŷ_i than in positive errors; therefore, the average systematic error (ASE) and the mean percent standard error (MPSE) were used (Zeng et al., 2017). The only difference among these metrics is the denominator; that is, instead of the observed value, the predicted value of the response is used in the denominator (Zeng et al., 2017; Huy et al., 2022).

ResultsTop

Relationships among tree biomass, carbon components, and predictors

Factor analysis was carried out with two groups of eight variables each. Group 1 included tree biomass components (Bst, Bbr, Ble, and Bba) and tree-level predictors (D, H, A, and CA) and group 2 included tree carbon components (Cst, Cbr, Cle and Cba) and tree predictors (D, H, A and CA). Factor analysis provided a small number of factors that account for most of the variability in the eight variables and supported an examination of the relationship between biomass and carbon components of trees. In both groups, two factors were extracted since two factors had eigenvalues greater than or equal to 1.0. Together, they accounted for 89.24% (group 1) and 89.30% (group 2) of the variability in the original data.

According to the results of the factor analysis, there were strong relationships between tree biomass and carbon components and tree predictors (D, H, A). They were placed close together in Fig. 4. They had a high coefficient (over 0.83) while CA had the lowest coefficient (below 0.30), which was shown in the first common factor equations (Factor 1) below. CA was excluded from further analysis because of its negligible influence on the biomass and carbon of the tree components.

Figure 4. The plot of Factor 1 vs. Factor 2 of the Factor Analysis. Left: Tree biomass components (Bst, Bbr, Ble and Bba) along with tree predictors (D, H, A, and CA). Right: Tree carbon sequestration components (Cst, Cbr, Cle and Cba) and tree predictors (D, H, A, and CA). The factor plots show the location of each variable; the variables furthest from the reference lines at 0 make the largest contribution to the factors, and the variables located close to each other are closely related.

The first common factor equations for tree biomass/carbon components are:

$\mathrm{Factor\; 1}=0.96088✕\mathrm{Bst}+0.95130✕\mathrm{Bba}+0.93313✕\mathrm{Ble}+0.92899✕\mathrm{Bbr}+0.82619✕A+0.97296✕D+0.92957✕H+0.27349✕\mathrm{CA}$ (11)

 

$\mathrm{Factor\; 1}=0.96183✕\mathrm{Cst}+0.953593✕\mathrm{Cba}+0.939325✕\mathrm{Cle}+0.927541✕\mathrm{Cbr}+0.825706✕A+0.972641✕D+0.928428✕H+0.272406✕\mathrm{CA}$ (12)

Independently fit models

The results for selecting independent models using the WNLS method for tree biomass and carbon components and AGB/AGC vs. different combinations of variables D, H, and D²H by LOOCV are presented in Tables 3 and 4. We also observed that tree age (A) was marginally significant for all components and total (p > 0.05). The selected models from the independent fitting of component and total AGB and AGC were as follows:

$\mathrm{Bst/Cst}=a✕{\left(D^{2}H\right)}^b$ (13)

 

$\mathrm{Bbr/Cbr}=a✕D^b$ (14)

 

$\mathrm{Ble/Cle}=a✕D^b$ (15)

 

$\mathrm{Bba/Cba}=a✕{\left(D^{2}H\right)}^b$ (16)

 

$\mathrm{AGB/AGC}=a✕{\left(D^{2}H\right)}^b$ (17)

The plots of observed vs. fitted values and weighted residuals vs. the fitted values (Fig. 5) depict a good model fit for independent fitting of component and total biomass. The results also showed that the model selection retained identical predictors or a combination of predictors when independent tree biomass and carbon components models were fitted (Table 3 and Table 4).

Figure 5. Plots of selected weighted non-linear models of tree biomass components Bst, Bbr, Ble and Bba (biomass of stem, branches, leaves, and bark, respectively) and total tree aboveground biomass (AGB) developed separately: Right: Fitted vs. Observed values; Left: Weighted residuals vs. Fitted values.

Table 3. Leave-One-Out Cross Validations (LOOCV) to select separate weighted nonlinear equations along with different combinations of tree predictor(s) for biomass components including Bst, Bbr, Ble and Bba (biomass of stem, branches, leaves and bark, respectively) and total tree aboveground biomass (AGB).

Model form	Weight variable	AIC	Adj. R2	Bias (%)	RMSE (kg)	MAPE (%)
Bst = a × Db	1/D	12.9	0.925	9.0	0.293	22.9
Bst = a × Db × Hc	1/D	3.0	0.956	5.9	0.269	19.3
Bst = a × (D2H)b	1/(D2H)0.2	3.8	0.958	6.4	0.238	18.2
Bst = a × Db × Ac*	1/D	14.1	0.923	4.0	0.306	23.3
Bst = a × Db × Hc × Ad*	1/D	4.7	0.953	5.4	0.286	19.7
Bst = a × (D2H)b×Ac*	1/(D2H)	1.2	0.945	-1.7	0.245	19.4
Bbr = a × Db	1/D	-14.1	0.842	-8.3	0.141	29.7
Bbr = a × Db × Hc*	1/D	-15.2	0.867	-6.1	0.143	30.0
Bbr = a × (D2H)b	1/(D2H)	-10.7	0.767	-13.9	0.161	34.5
Bbr = a × Db × Ac*	1/D	-12.2	0.832	-8.8	0.151	30.9
Bbr = a × Db × Hc* × Ad*	1/D	-13.3	0.859	-6.5	0.150	31.1
Bbr = a × (D2H)b×Ac*	1/(D2H)	-8.8	0.747	-14.7	0.173	36.2
Ble = a × Db	1/Dδ	1.1	0.732	-11.4	0.214	34.5
Ble = a × Db × Hc*	1/D	2.5	0.721	-11.5	0.240	37.5
Ble = a × (D2H)b	1/(D2H)	8.9	0.702	-18.1	0.213	39.6
Ble = a × Db × Ac*	1/D	2.9	0.720	-11.4	0.229	36.0
Ble = a × Db × Hc* × Ad*	1/D	4.3	0.710	-12.1	0.259	39.6
Ble = a × (D2H)b×Ac*	1/(D2H)	10.5	0.668	-18.7	0.234	42.0
Bba = a × Db	1/D	-15.8	0.797	-7.7	0.148	30.0
Bba = a × Db × Hc*	1/D	-18.9	0.837	-7.3	0.124	24.0
Bba= a × (D2H)b	1/(D2H)	-29.9	0.844	-6.4	0.119	24.0
Bba = a × Db × Ac	1/D	-20.1	0.840	-3.9	0.135	29.9
Bba = a × Db × Hc × Ad	1/D	-26.9	0.897	-2.1	0.125	27.1
Bba = a × (D2H)b×Ac	1/(D2H)	-33.4	0.889	-4.1	0.115	26.1
AGB = a × Db × Hc*	1/D	42.2	0.934	1.8	0.640	21.3
AGB = a × (D2H)b	1/(D2H)0.2	42.7	0.934	1.8	0.589	20.1
AGB = a × Db × Ac*	1/D	45.6	0.922	1.5	0.590	21.6
AGB = a × Db × Hc × Ad*	1/D	43.8	0.932	2.0	0.695	22.7
AGB = a × (D2H)b×Ac*	1/(D2H)	45.9	0.912	-5.2	0.684	23.0

In LOOCV, the dataset was split into n-1 samples used for developing models, the remaining one sample used for validation, calculation of AIC, Adj. R², Bias, RMSE, MAPE; finally, all of those statistics averaged over n realizations. *: parameter with p > 0.05. In bold, models selected based on LOOCV statistics and diagnostic plots.

Table 4. Leave-One-Out Cross Validations (LOOCV) to select separate weighted nonlinear equations along with different combinations of tree predictor/s for carbon components including Cst, Cbr, Cle and Cba (carbon sequestration of stem, branches, leaves and bark, respectively) and total tree aboveground carbon (AGC).

Model form	Weight variable	AIC	Adj. R2	Bias (%)	RMSE (kg)	MAPE (%)
Cst = a × Db	1/D	-18.6	0.924	3.5	0.135	22.2
Cst = a × Db × Hc	1/D	-26.2	0.950	4.9	0.134	19.4
Cst = a × (D2H)b	1/(D2H)	-30.1	0.946	-1.6	0.118	19.6
Cst = a × Db × Ac*	1/D	-17.7	0.925	3.0	0.146	22.8
Cst = a × Db × Hc × Ad*	1/D	-24.9	0.949	4.3	0.142	19.8
Cst = a × (D2H)b×Ac*	1/(D2H)	-29.4	0.943	-1.8	0.121	19.2
Cbr = a × Db	1/D	-44.1	0.830	-9.0	0.067	29.8
Cbr = a × Db × Hc*	1/D	-45.3	0.859	-6.7	0.069	30.6
Cbr = a × (D2H)b	1/(D2H)	-40.9	0.752	-14.2	0.077	34.4
Cbr = a × Db × Ac*	1/D	-42.2	0.819	-9.5	0.072	31.2
Cbr = a × Db × Hc* × Ad*	1/D	-43.5	0.850	-7.1	0.074	31.9
Cbr = a × (D2H)b×Ac*	1/(D2H)	-38.9	0.735	-15.1	0.083	36.4
Cle = a × Db	1/D	-29.9	0.741	-10.9	0.100	33.9
Cle = a × Db × Hc*	1/D	-28.5	0.735	-11.3	0.113	37.2
Cle = a × (D2H)b	1/(D2H)	-21.6	0.715	-18.2	0.101	39.6
Cle = a × Db × Ac*	1/D	-28.0	0.732	-10.8	0.108	35.5
Cle = a × Db × Hc* × Ad*	1/D	-26.8	0.726	-11.8	0.123	39.5
Cle = a × (D2H)b×Ac*	1/(D2H)	-19.9	0.685	-18.9	0.112	42.4
Cba = a × Db	1/D	-48.9	0.798	-7.5	0.066	29.0
Cba = a × Db × Hc*	1/D	-51.7	0.834	-7.0	0.056	23.2
Cba= a × (D2H)b	1/(D2H)	-62.3	0.842	-6.3	0.053	22.9
Cba = a × Db × Ac	1/D	-52.1	0.831	-4.2	0.062	29.1
Cba = a × Db × Hc × Ad	1/D	-57.4	0.882	-2.7	0.060	27.0
Cba = a × (D2H)b×Ac	1/(D2H)	-65.5	0.878	-4.3	0.054	26.0
AGC = a × Db	1/D	13.2	0.922	0.8	0.273	20.8
AGC = a × Db × Hc*	1/D	12.2	0.929	1.3	0.316	21.2
AGC= a × (D2H)b	1/(D2H)0.2	13.0	0.929	0.9	0.292	20.1
AGC = a × Db × Ac*	1/D	15.0	0.918	0.9	0.293	21.7
AGC = a × Db × Hc* × Ad*	1/D	14.0	0.926	1.4	0.344	22.6
AGC = a × (D2H)b×Ac*	1/(D2H)	15.6	0.910	-5.4	0.331	23.5

In LOOCV, the dataset was split into n-1 samples used for developing models, the remaining one sample used for validation, calculation of AIC, Adj. R² Bias, RMSE, MAPE; finally, all of those statistics averaged over n realizations. *: Parameter with p > 0.05. In bold, models selected based on LOOCV statistics and diagnostic plots.

The carbon fraction of the tree species studied

The results of calculating the average carbon fraction (CF = Carbon/Biomass weight) for tree biomass components and total and the variation within the 95% confidence range are shown in the box plot in Fig. 6.

Figure 6. Box plot of carbon fraction (CF) by total tree aboveground biomass (AGB) and its components. Points and numbers within the boxplots represent the means of CF for the tree components and the total. Sample size n = 22.

Simultaneous modeling

The tree age was excluded from further analysis based on the results of independently fit models. Sixteen WNSUR equation systems were developed using D, H, and their combinations as predictors to predict AGB, AGC, and their components (Tables S1 and S2 [suppl]). Unlike the independent model fit, the results showed that the model selection retained the different predictors or a combination of predictors when these models were fit simultaneously. Based on the errors statistics of LOOCV, we found two optimal combinations of tree predictors for simultaneous modeling systems for tree biomass and carbon components and the totals are as follows:

$\mathrm{Bst}=a_1✕D^{\mathrm{b1}}$ (18)

 

$\mathrm{Bbr}=a_2✕D^{\mathrm{b2}}$ (19)

 

$\mathrm{Ble}=a_3✕D^{\mathrm{b3}}$ (20)

 

$\mathrm{Bba}=a_4✕{\left(D^{2}H\right)}^{\mathrm{b4}}$ (21)

 

$\mathrm{AGB}=\mathrm{Bst}+\mathrm{Bbr}+\mathrm{Ble}+\mathrm{Bba}=a_1✕D^{\mathrm{b1}}+a_2✕D^{\mathrm{b2}}+a_3✕D^{\mathrm{b3}}+a_4✕{\left(D^{2}H\right)}^{\mathrm{b4}}$ (22)

 

$\mathrm{Cst}=a_1✕{\left(D^{2}H\right)}^{\mathrm{b1}}$ (23)

 

$\mathrm{Cbr}=a_2✕D^{\mathrm{b2}}$ (24)

 

$\mathrm{Cle}=a_3✕{\left(D^{2}H\right)}^{\mathrm{b3}}$ (25)

 

$\mathrm{Cba}=a_4✕{\left(D^{2}H\right)}^{\mathrm{b4}}$ (26)

 

$\mathrm{AGC}=\mathrm{Cst}+\mathrm{Cbr}+\mathrm{Cle}+\mathrm{Cba}=a_1✕{\left(D^{2}H\right)}^{\mathrm{b1}}+a_2✕D^{\mathrm{b2}}+a_3✕{\left(D^{2}H\right)}^{\mathrm{b3}}+a_4✕{\left(D^{2}H\right)}^{\mathrm{b4}}$ (27)

The parameters of the two selected simultaneous modeling systems that were fit by the WNSUR method using the entire dataset are presented in Table 5. In addition to the simplicity of application, since measuring tree H is costly, two modeling systems for simultaneously predicting tree biomass and carbon were created with only D (Table 6).

Table 5. Model equations and estimated parameters of the selected modeling systems for simultaneous estimation of AGB/AGC and its components using the WNSUR method (based on the entire dataset).

Prediction of tree AGB or AGC and its components	Parameters	Estimate ± Std. Error	RMSE (kg tree-1)	Adj. R2
AGB and its components
Bst = a1×Db1	a1	0.03203 ± 0.01450	0.349	0.936
	b1	2.50566 ± 0.24550
Bbr = a2×Db2	a2	0.04571 ± 0.02050	0.173	0.865
	b2	1.75896 ± 0.24880
Ble = a3×Db3	a3	0.07274 ± 0.03410	0.238	0.774
	b3	1.57179 ± 0.25980
Bba = a4×(D2H)b4	a4	0.01545 ± 0.00766	0.156	0.867
	b4	0.80636 ± 0.09310
AGB = Bst + Bbr + Ble + Bba			0.692	0.934
AGC and its components
Cst = a1×(D2H)b1	a1	0.01399 ± 0.00367	0.134	0.956
	b1	0.91022 ± 0.05050
Cbr = a2×Db2	a2	0.01514 ± 0.00685	0.087	0.847
	b2	1.97080 ± 0.24890
Cle = a3×(D2H)b3	a3	0.05770 ± 0.02160	0.113	0.782
	b3	0.46021 ± 0.04530
Cba = a4×(D2H)b4	a4	0.00702 ± 0.00342	0.071	0.866
	b4	0.80219 ± 0.09160
AGC = Cst + Cbr + Cle + Cba			0.328	0.934

Bst/Cst, Bbr/Cbr, Ble/Cle, Bba/Cba and AGB/AGC are biomass/carbon of stem, branches, leaves, bark and total aboveground biomass/carbon, respectively. All parameters with p < 0.05.

Table 6. Model equations and estimated parameters for a sole predictor of the diameter at breast height (D) for the simultaneous estimation of AGB/AGC and its components using the WNSUR method (based on the entire dataset).

Prediction of tree AGB or AGC and its components	Parameters	Estimate ± Std. Error	RMSE (kg tree-1)	Adj. R2
AGB and its components
Bst = a1×Db1	a1	0.04884 ± 0.01690	0.351	0.935
	b1	2.27942 ± 0.18910
Bbr = a2×Db2	a2	0.05699 ± 0.02530	0.175	0.863
	b2	1.64333 ± 0.24450
Ble = a3×Db3	a3	0.11233 ± 0.04200	0.240	0.772
	b3	1.34080 ± 0.20790
Bba = a4×Db4	a4	0.02612 ± 0.01230	0.178	0.827
	b4	1.96718 ± 0.26420
AGB = Bst + Bbr + Ble + Bba			0.729	0.916
AGC and its components
Cst = a1×Db1	a1	0.02224 ± 0.00540	0.165	0.935
	b1	2.29456 ± 0.13520
Cbr = a2×Db2	a2	0.02892 ± 0.01330	0.086	0.853
	b2	1.59982 ± 0.25500
Cle = a3×Db3	a3	0.04956 ± 0.01770	0.113	0.782
	b3	1.38225 ± 0.20110
Cba = a4×Db4	a4	0.01286 ± 0.00613	0.081	0.826
	b4	1.91162 ± 0.26470
AGC = Cst + Cbr + Cle + Cba			0.352	0.923

Bst/Cst, Bbr/Cbr, Ble/Cle, Bba/Cba and AGB/AGC are biomass/carbon of stem, branches, leaves, bark and total aboveground biomass/carbon, respectively. All parameters with p < 0.05.

Both WNLS and WNSUR addressed the issue of heteroscedasticity of the residuals. The WNLS residuals fluctuated uniformly along with the fitted values of the models, whereas the WNSUR model reduced the variability of the residuals significantly compared to the WNLS model (Fig. 7). The result of comparison of LOOCV statistics between two methods of independent WNLS and simultaneous WNSUR to fit selected AGB and AGC models is shown in Table 7.

Figure 7. Plots of fitted vs. observed AGC (left) and weighted residuals vs. fitted AGC (right) for the model fitted with the entire dataset. Comparison of two methods of independent weighted non-linear (WNLS) and simultaneous weighted non-linear SUR (WNSUR).

Table 7. Model equations and estimated parameters for a sole predictor of the diameter at breast height (D) for the simultaneous estimation of AGB/AGC and its components using the WNSUR method (based on the entire dataset).

Prediction of tree AGB or AGC	Method to fit the model	Selected model form	Bias (%)	RMSE (kg tree-1)	MAPE (%)
AGB	Independent WNLS	AGB = a × (D2H)b	1.8	0.589	20.1
	Simultaneous WNSUR	AGB = Bst + Bbr + Ble + Bba = a1×Db1 + a2×Db2 + a3×Db3 + a4×(D2H)b4	2.2	0.564	18.8
AGC	Independent WNLS	AGC = a × (D2H)b	0.9	0.292	20.1
	Simultaneous WNSUR	AGC = Cst + Cbr + Cle + Cba = a1×(D2H)b1 + a2×Db2 + a3×(D2H)b3 + a4×(D2H)b4	15.5	0.455	16.2

Bst/Cst, Bbr/Cbr, Ble/Cle, Bba/Cba and AGB/AGC are biomass/carbon of stem, branches, leaves, bark and total aboveground biomass/carbon, respectively.