A comparison between traditional ordinary least-squares regression and three methods for enforcing additivity in biomass equations using a sample of Pinus radiata trees

Background: Additivity has long been recognised as a desirable property of systems of equations to predict the biomass of components and the whole tree. However, most tree biomass studies report biomass equations fitted using traditional ordinary least-squares regression. Therefore, we aimed to develop models to estimate components, subtotals and aboveground total biomass for a Pinus radiata D.Don biomass dataset using traditional linear and nonlinear ordinary leastsquares regressions, and to contrast these equations with the additive procedures of biomass estimation. Methods: A total of 24 ten-year-old trees were felled to assess above-ground biomass. Two broad procedures were implemented for biomass modelling: (a) independent; and (b) additive. For the independent procedure, traditional linear models (LINOLS) with scaled power transformations and y-intercepts and nonlinear power models (NLINOLS) without y-intercepts were compared. The best linear (transformed) models from the independent procedure were further tested in three different additive structures (LINADD1, LINADD2, and LINADD3). All models were evaluated using goodness-of-fit statistics, standard errors of estimates, and residual plots. Results: The LINOLS with scaled power transformations and y-intercepts performed better for all components, subtotals and total above-ground biomass in contrast to NLINOLS that lacked y-intercepts. The additive model (LINADD3) in a joint generalised linear least-squares regression, also called seemingly unrelated regression (SUR), provided the best goodnessof-fit statistics and residual plots for four out of six components (stem, branch, new foliage and old foliage), two out of three subtotals (foliage and crown), and above-ground total biomass compared to other methods. However, bark, cone and bole biomass were better predicted by the LINOLS method. Conclusions: SUR was the best method to predict biomass for the 24-tree dataset because it provided the best goodnessof-fit statistics with unbiased estimates for 7 out of 10 biomass components. This study may assist silviculturists and forest managers to overcome one of the main problems when using biomass equations fitted independently for each tree component, which is that the sum of the biomasses of the predicted tree components does not necessarily add to the total biomass, as the additive biomass models do. New Zealand Journal of Forestry Science KC et al. New Zealand Journal of Forestry Science (2020) 50:7 https://doi.org/10.33494/nzjfs502020x90x E-ISSN: 1179-5395 published on-line: 25/10/2020 © The Author(s). 2020 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Research Article Open Access Pinus radiata D.Don, native to California, is a widely planted commercial tree species in the Southern Hemisphere, including New Zealand, Australia, Chile, Spain and South Africa (Lavery & Mead 2000; Mead Introduction Forests play a vital role in the carbon cycle to mitigate climate change by accumulating and sequestering atmospheric carbon dioxide (CO2) (Houghton 1991).


Introduction
Forests play a vital role in the carbon cycle to mitigate climate change by accumulating and sequestering atmospheric carbon dioxide (CO 2 ) (Houghton 1991).
Keywords: above-ground; additive; biomass; linear; nonlinear; radiata pine; seemingly unrelated regression 2013). It is grown primarily for timber production, as this species is versatile, fast-growing, and has a wide range of end uses (Lavery & Mead 2000;Lewis et al. 1993;Rogers 2002;Sutton 1999;Toro & Gessel 1999). The global plantation area of P. radiata is now more than 4.2 million hectares (Mead 2013). In New Zealand, it is the predominant planted species, and accounts for about 90% of the total 1.7 million hectares of forest plantations (Nixon et al. 2017). Plantation forests in New Zealand have not only been recognised as providing financial returns from traditional wood products, but also as providing environmental services by accumulating biomass and sequestering a substantial amount of carbon. To quantify such benefits, a precise biomass model with a required level of accuracy is essential.
Tree biomass estimation is required by scientists and practitioners alike as a surrogate of ecosystem production, product outturn and carbon accounting, among others. Biomass modelling is important for estimating carbon sequestration of forest ecosystems, as individual tree biomass or its components is aggregated to yield the stand biomass (Zheng et al. 2015). Allometric models are commonly used to assess the biomass accumulated in forests. Allometric relationships can be developed from destructive sampling by using several forms of regression equations. Generally, biomass equations are fitted in linear form using logarithmic transformation of B and D of the form, B = aD b , where B is the biomass of the tree, or its components, and D is the diameter of the tree (Baskerville 1972;Beauchamp & Olson 1973;Canadell et al. 1988;Santa Regina et al. 1997;Sprugel 1983). Clutter et al. (1983) explained various linear and nonlinear additive regression models to estimate the biomass of an individual tree or its components. The additivity of biomass equations has long been recognised as a desirable property, so that predictions of tree components added together equals predictions of total tree biomass (Cunia & Briggs 1984, 1985Parresol 1999Parresol , 2001. Three procedures for forcing additivity have been proposed (Cunia & Briggs 1985;Parresol 1999): (a) adding the best regression functions of the components' biomass to determine the total biomass regression function; (b) using the same independent variables for each component; and (c) using joint generalised least-squares regression, also known as seemingly unrelated regression (SUR), in which statistical dependencies among sample data are accounted for by forcing constraints on the regression coefficients. These three procedures have been extensively applied for estimating tree biomass around the world (Canga et al. 2013;Návar, González et al. 2004). The additive procedure in SUR has been mostly used for biomass modelling of single species (Cunia & Briggs 1984, 1985Green & Reed 1985;Parresol 1999;Zheng et al. 2015).
Over the last 50 years, a substantial number of biomass studies for P. radiata have been undertaken in New Zealand (Beets & Madgwick 1988;Beets et al. 2007;Beets & Pollock 1987;Cromer et al. 1985;Madgwick 1983Madgwick , 1985Madgwick , 1994Madgwick et al. 1977;Mead et al. 1984;Moore 2010;Webber & Madgwick 1983;Will 1964). Previous studies for P. radiata in New Zealand aimed to find the best biomass equations using various functional linear and nonlinear forms, with models generally fitted separately for each individual biomass component and for the whole tree. Separately calculated biomass equations ignore inherent correlation among the component equations measured on the same tree (Kozak 1970;Parresol 1999). Simultaneous fits with related equations using additive procedures have greater statistical efficiency, as they take into account statistical dependencies among biomass components in parameter estimation recorded from the same tree (Bi et al. 2010;Bi et al. 2004;Carvalho & Parresol 2003;Parresol 1999Parresol , 2001. Two country specific systems of additive biomass equations were developed for P. radiata using routinely measured stand variables from Australia and New Zealand (Bi et al. 2010). It has been noted that prediction accuracy varies across methodological differences and uncertainties associated over a range of stand variables (Bi et al. 2010;Moore 2010). As there is uncertainty about how to better meet additivity requirements, this study was undertaken to compare traditional linear and nonlinear ordinary least-squares regressions, and additive procedures in the estimation of tree component and total biomass for a dataset composed of 24 trees of P. radiata.

Study site and experiment
This study was carried out in the Canterbury region of New Zealand, planted with P. radiata in 2005, in a forestry trial designed to test the effect of stocking, genetics, fertiliser application, and follow-up weed control treatment on productivity and wood quality (Mason 2008). The site is located at latitude 43° 37.2′ S and longitude 172° 20.4′ E, and about 45 m above sea level on a flat landscape. The site has a mean annual air temperature between 11 and 13 °C with a monthly minimum (July) of −2 to +4 °C and a monthly maximum (January) of 20 to 23 °C (Macara 2016). Annual rainfall is about 618 mm with a monthly range between 38 and 68 mm (Macara 2016). The experiment consisted of 48 permanent plots with a randomised complete block split-split design, with the arrangement of factors within four complete blocks (Mason 2008). During the summer of 2015 to 2016, 24 ten-year-old trees of P. radiata were harvested and measured from six plots of the trial, and within each plot four trees were felled. These plots consisted of three levels of stocking (625, 1250 and 2500 stems ha −1 ), two levels of follow-up weed control treatment (herbicide and no chemical treatment) and two clones (1 and 2).

Biomass data
Trees were felled at ground level. The over-bark diameter of each tree at breast height was recorded at 1.4 m. Total tree height was measured from ground level to the tip of the tree bole. For each tree, the components were separated into stem, branch, bark, foliage, and cones. Needles and twigs less than 1 cm in diameter were considered foliage, and this was separated into "new" and "old" foliage. The total fresh mass of all components including subsamples were measured immediately after felling, using a portable balance. All the cones and small branches were weighed separately. The logs were separated into small pieces and weighed fresh in the field. A subsample of stem discs with bark (cut at the 1.4 m section and every 2 m upwards in the stem) and subsamples of all other components, were weighed to determine fresh weight in the field. Subsamples were dried in an oven at 70 °C until constant mass was achieved, and then this weight was recorded. Dry mass of each component was calculated as the fresh mass recorded in the field for that component multiplied by the ratio of subsample dry to fresh mass (Eq. 1): (1) where Y is the total dry mass (kg), DW and FW refers to the sub sampled dry and fresh mass (kg) respectively, TFW is the total fresh mass (kg), and i is the tree component such as stem, bark, branch, new foliage, old foliage and cones.
Descriptive statistics of the trees including components, sub-total and above-ground total biomass are shown in Table 1. The notations and definitions used in this manuscript are explained in the Abbreviations section.

Variance stabilisation
Biomass data generally exhibit non-constant variance in model residuals (Parresol 1993(Parresol , 2001. When developing predictive equations, variance can be stabilised either by providing a weight function or by using transformations (Parresol 1993(Parresol , 2001. Curvilinearity and heterogeneity in variance of all linear models were reduced by transforming the response as well as explanatory variables using scaled-power transformations (Eq. 2), widely known as Box-Cox transformations (Box & Cox 1964). The predicted values of these models were back transformed to the original form using Eq. 3. A similar variance stabilisation technique was implemented by Zheng (2015) while using the additive procedure of biomass modelling for Quercus variabilis in northern China. (2) where Y (λ) is the transformed variable, and λ is a coefficient of the transformed variable that varies normally between −3 and +5 (Cook & Weisberg 2009), Y' is the back-transformed variable. A λ term is chosen to make the frequency distribution of each variable as close to normal as possible, thus promoting linear relationships and stabilising variance.

Model assessment and evaluation
In this study, a dataset consisting of 24 trees was used to evaluate the fitting bias, precision, and validity of models using the following goodness-of-fit statistics: root mean square error (RMSE), mean absolute bias (MAB), mean prediction error (MPE), residual standard error (RSE), coefficient of variation (CV), coefficient of determination (R 2 ), index of agreement (IOA), and Akaike information criterion (AIC). Models were considered better with small AIC, RMSE, MAB, MPE, RSE, and CV of the residuals,  1: Descriptive statistics for the 24-felled trees used for developing regression models, and components, sub-total and above-ground total biomass. and large R 2 and IOA. The interpretation of these fitting statistics can be found in Von Gadow and Hui (2001) and Goicoa et al. (2011). In addition, model performance was assessed by residual plots and histograms of residuals.

Modelling procedure
In this research, two procedures were implemented to estimate components, subtotals and above-ground total biomass: (1) independent; and (2) additive. All models were fitted to estimate biomass in terms of kg tree −1 .

Independent procedure for biomass estimation
In this procedure, biomass equations were fitted independently using traditional linear ordinary leastsquares regressions with scaled power transformations and y-intercepts (denoted as, LINOLS; Eq. 4) and nonlinear ordinary least-squares power equations that lacked y-intercepts (denoted as, NLINOLS; Eq. 5). The mathematical specifications of these models are as follows (Parresol 1999;Zeng 2011;Zianis et al. 2005). (4) where f l (X l , β l ) is the regression function for the above ground biomass or one of its components, X l are tree dimension variables such as D, H and CrL (l = 1, 2, . . . . , p) while β l denote the regression coefficients. Each component equation contained its own independent variables. All components, subtotals and AGT biomass equations were fitted separately using the lm and nls function of R statistical software (R Core Team, 2018), for linear and nonlinear regressions, respectively.

Additive procedure of biomass estimation
In this procedure, biomass equations were fitted based on three additive procedures, described and compared by Parresol (1999Parresol ( , 2001. The additivity requirement to estimate total tree biomass is ensured by (a) adding the separately calculated best regression functions of each component, (b) using the same independent variables for each component, and (c) using joint generalised leastsquares methods, also known as seemingly unrelated regression (SUR). In SUR, statistical dependencies among components are forced by constraining regression coefficients (Cunia & Briggs 1985;Parresol 1999). In this study, four restrictions were provided for the SUR model: (1) foliage; (2) crown; (3) bole; and (4) AGT, as illustrated in Figure 1. For example, foliage biomass is the sum of NF, OF and cone biomass (Eq. 6). Mathematically, the additive system of biomass equations in additive error terms with cross-equation correlation is specified in Eq. 6 where Ŷ i represents the predicted biomass of a given component and f i (X i , β i ) is a regression function for the biomass component, (i = cone, new foliage, old foliage, branch, bark and stem, foliage, crown, bole and AGT biomass). The residual is ε ij for the i th equation and j is an index for component. All additive biomass equations were fitted in the R statistical software (R Core Team 2018) using the systemfit package (Henningsen & Hamann 2007).
In the first additive procedure, the additivity was ensured by adding individually calculated best regression functions of each component to give a total biomass regression function (Cunia & Briggs 1985;Parresol 1999). The best regression functions obtained from the independent procedure of biomass modelling that were fitted separately for each component given in Table A.1 were used. The additive structure of this model, denoted as LINADD1, is specified in Eq. 7. In the second additive procedure, additivity was implemented by using the same explanatory variables for each component. For this, the most frequent independent variable (D) was selected from the best linear regression function as it was best fitted for stem, bark, foliage, bole, and AGT (Table A.1). Using D as an independent variable (6) for all components, the additive structure of the model, denoted as LINADD2, is specified in Eq. 8.
In the third additive procedure, we used different explanatory variables in a joint generalised linear leastsquares regression, known as SUR (Cunia & Briggs 1985;Parresol 1999). For this, best-fitted explanatory variables from the independent procedure of biomass modelling were used for stem, cone, branch, NF, and OF (Table A. 1). We used the second-best regression D 2 H as an independent variable for bark (data not shown). The additive structure of the model, denoted as LINADD3, is specified in Eq. 9.

Comparison of fitted equations for components, subtotals and AGT
Tested LINOLS and NLINOLS equations with their bestfit results are given in Table A.1, and fitted statistics with their regression estimates are presented in Table A.2. We attempted to take into account follow-up herbicide, stocking, and clone factors into all models as dummy variables. These were found to be non-significant (P>0.05) so were discarded from all subsequent modelling. In comparison, LINOLS provided relatively higher R 2 values than NLINOLS for all, except for branch and cone biomass (Table A.2). However, plotting residuals with predicted values and with other variables demonstrated that NLINOLS regression was unsuitable for these two components (data not shown). Therefore, overall, the best fitted LINOLS model according to goodness-of-fit statistics and residual plots were Eq. (i) for stem, bark, foliage, bole and AGT biomass, Eq. (ii) for cone biomass, Eq. (iii) for branch biomass, Eq. (viii) for NF and crown biomass, and Eq. (ix) for OF biomass (Table  A. 1). Finally, these selected LINOLS models were further tested in the additive process of biomass estimations. The estimated coefficients for six components, three subtotals, and AGT using four methods (LINOLS, LINADD1, LINADD2 and LINADD3) are presented in Table 2 and their goodness of fit statistics are given in Table 3. The distribution of residuals with predicted values of the fitted best models for six components, three subtotals, and AGT are given in Fig. 2. The LINADD3 fitted in SUR was considered best to predict stem (Eq. 10), branch (Eq. 11), NF (Eq. 13), OF (Eq. 14) biomass as it provided the better-fitting statistics when compared to the other three equations (Table 3). For stem, LINADD3 simultaneously decreased the RMSE, RSE, and CV by 0.1%, MPE by 0.2% while R 2 increased by 0.005%, compared to the other three equations. For branches, LINADD3 provided a marginal decrease in the goodness-of-fit statistics (e.g. RMSE, RSE, and CV by 1%) in contrast to LINOLS and LINADD1. For NF, LINADD3 model recorded a decrease in fitting statistics (e.g. RMSE by 3.7%, RSE by 1.43%), in contrast to the other three    equations. For OF, the LINADD3 provided a decrease in RMSE (by 10.9%), MAB (by 12.7%), MPE (by 20.3%), RSE (by 9.6%) and CV (by 10.9%), and an increase in R 2 (by 1.2%) and IOA (by 0.3%), in contrast to the other three equations (Table 3). The LINOLS was considered the best to predict bark (Eq. 12), cone (Eq. 15) and bole (Eq. 12) biomass exhibiting relatively good fit statistics as compared to SUR models (  (Table 3).
The additive equation fitted in SUR was considered the best to predict foliage (Eq. 16) and crown (Eq. 17) biomass as LINADD3 provided better fit statistics ( Table  3). For foliage, LINADD3 provided a decrease in RMSE, MAB, MPE and CV by 2.1%, 1.6%, 4.1% and 2.1%, respectively, and R 2 by 0.3% and IOA by 0.1% compared to other three equations (Table 3). For crown, using LINADD3, average fitting statistics decreased by 21.5% for RMSE, 16.7% for MAB, 36.3% for MPE, 18.1% for RSE, and 21.5% for CV; and increased R 2 by 7% and IOA by 1.7%, in contrast to other three equations. The LINADD3 fitted in SUR (Eq. 19) was considered the best to predict AGT biomass as it provided a decrease in RMSE by 10.1%, MAB by 9.3%, MPE by 18.7, RSE by 5.8% and CV by 10.1%; and an increase in R 2 by 1.2% and IOA by 0.3%, in contrast to the other three methods.

Discussion
In this study, follow-up herbicide, stocking, and clone factors fitted into models as dummy variables were all non-significant. The possible reason for this insignificance could be that the site or silvicultural effects KC et al. New Zealand Journal of Forestry Science (2020)   The solid black horizontal line across zero represent baseline and the dotted red line is LOESS curve.
were explained by the tree size variables because the silvicultural treatments may simply speed up/slow the growth of the trees (Beets & Pollock 1987;Moore 2010), as opposed to changing their allometry. The biomass of P. radiata plantation changes with a wide range of climatic, edaphic, silvicultural and genetic factors (Beets & Pollock 1987;Bi et al. 2010;Moore 2010). Another reason of this lack of significance could be the small sample size (Duncanson et al. 2015) as the parameterisation of allometric equations depends significantly on the size of sample.
In this study, the best model of each component selected from independent procedures provided the logical base equations for further tests in the additivity of biomass equations. The same approach was also used by other researchers (Magalhães & Seifert 2015;Návar et al. 2002) to utilise additive properties in biomass modelling. Biomass additivity reduces the discrepancy between the sum of predicted values for components and those for a total tree (Kozak 1970), and it has long been documented as a desirable property of systems of equations to predict total tree biomass (Bi et al. 2004;Parresol 2001). Three procedures were implemented for the additivity in the biomass model (Parresol 1999(Parresol , 2001: (1) using a separately calculated best linear function of the biomass of the components (best linear functions were D, D and H, D, D and CrL, D and CrL 2 , and D 2 for stem, branch, bark, NF, OF, and cone biomass, respectively); (2) using the most frequently observed predictor (D) as the same independent variable for all components; and (3) using different independent variables for each component by forcing four linear restrictions on the regression coefficients, the SUR technique. The additivity of biomass equations to predict biomass of components, subtotal, and AGT has been explained in some other studies (Carvalho & KC et al. New Zealand Journal of Forestry Science (2020) (Beets & Madgwick 1988;Beets & Pollock 1987), and nutrients and silvicultural practices (Beets & Madgwick 1988;Madgwick 1985;Mead et al. 1984) may influence the biomass allometry.

Conclusions
Two procedures for biomass modelling were compared, namely, independent and additive. For the independent procedure of biomass modelling, LINOLS models with scaled power transformations and y-intercepts and NLINOLS power models that lacked y-intercepts were compared for six components, three subtotals, and AGT biomass. The LINOLS models with scaled power transformations and y-intercepts provided superior results in contrast to NLINOLS power models without y-intercepts. The best-fitted component equations from LINOLS models were further tested in an additive procedure. A system of additive equations in SUR with different independent variables for each component (LINADD3) showed better performance than LINOLS, LINADD1, and LINADD2. Besides, the linear SUR model provided comparatively unbiased estimates for stem (Eq. 10), branch (Eq. 11), NF (Eq. 13), OF (Eq. 14), foliage (Eq. 16), crown (Eq. 17), and AGT (Eq. 19), while LINOLS showed comparatively better fitting statistics for bark (Eq. 12), cone (Eq. 15) and bole biomass (Eq. 18) for the dataset of this study. Since seven out of ten biomass components were well fitted with SUR that provided lower variance by taking account of the existence of correlations among residuals of the component equations, we suggest that SUR could be a superior method for fitting biomass equations.
A linear SUR model (LINADD3) of this study provided better results than LINOLS, LINADD1, and LINADD2, in terms of goodness-of-fit statistics, standard error of estimates and residual plots. LINADD3 fitted in SUR was superior to the other two additive models since it considered the correlation between each component equation, and provided greater statistical efficiencies (Carvalho & Parresol 2003). In contrast to our results, a study reported that the additive model (denoted as CON) that used the same independent variable for all components, similar to our LINADD2 model, was statistically superior to the linear and nonlinear SUR model with the different independent variables in parameter restriction (Magalhães & Seifert 2015). However, the authors (Magalhães & Seifert 2015) indicated that the CON method had the limitation that it did not take into account the correlations among plant parts.
Applying SUR to the system of additive equations with the same explanatory variables for each component does not provide precise estimation of biomass (Srivastava & Giles 1987). Therefore, the LINADD3, a SUR model that consisted of different explanatory variables for each component is consistent with that of Srivastava and Giles (1987), which was more effective than the other two additive models. Model LINADD1 also consisted of the same explanatory variables for two-component equations such as stem and bark, and LINADD2 consisted of the same independent variable for all component equations. Therefore, LINADD1 and LINADD2 were not effective compared to LINADD3. In addition, the individually calculated best equations for each component (LINOLS) provided the least efficient biomass estimates for all components except for bark, cone and bole biomass, compared with the linear SUR model (LINADD3). Researchers recommended using SUR to estimate biomass as it provides greater statistical efficiency than separately calculated equations for each component (Bi et al. 2010;Bi et al. 2004;Kozak 1970;Návar et al. 2002;Návar, Méndez et al. 2004;Parresol 2001).
Although this study focused on testing different procedures to fit biomass equations and highlighted the importance of SUR, there could be concerns over the applicability of the resulting models given the size of the dataset from the trial. While applying these models, it is advisable to consider that the dataset used was relatively small with only 24 trees of the same age class, with D (8.2 to 28 cm), H (8.85 to 13.77 m), and CrL (0.2 to 6 m). Small sample sizes may provide biased estimates, as the allometric parameters are sample size dependent (Duncanson et al. 2015). In addition, the models developed were based on only two out of five clones, selected across a range of treatment plots. Extrapolation should be avoided as uncertain prediction errors are expected from the selected models. The site characteristics (Duncanson et al. 2015;Madgwick 1994;Mason 2008), growth inputs (Dong et al. 2015