Eucalyptus grandis and Eucalyptus dunnii are the most planted tree species in Uruguay. Anticipating information about the quantity and quality of wood is important for managing intensive forest plantation. The estimate of merchantable and total wood volume is an essential tool in forest planning and management. The aim of this study was to evaluate four systems of taper and merchantable volume that consisted in a taper, a merchantable volume and a total tree volume function. A modified secondorder continuous autoregressive error structure corrected the inherent serial autocorrelation of different observations in one tree. Taper and volume equations were fitted simultaneously after autocorrelation correction by full information maximum likelihood method. The segmented system proposed by Fang et al. (2000) produced the best fit as it explained more than 98% of the taper, merchantable volume and total volume variability for both species. In addition, precision of the segmented system was compared with and without incorporating stand density as a variable. Results of this analysis showed that for E. grandis, the predictive accuracy of the model was improved by including the stand density variable, whereas for E. dunnii this variable was not statistically significant. This modelling framework provides an improvement in taper and tree volume predictions for E. dunnii and E. grandis in Uruguay. The possibilities offered by this methodology could be of interest for its application in countries where fast growing plantations are managed.
The total volume of individual trees is commonly estimated using models that take into account the tree diameter and total height, but this method does not allow to estimate the distribution of products according to their final commercial use. To solve this problem, stem taper function equations are usually developed (Kozak 1988, Bi 2000). The use of taper functions ensures an accurate description of the diametric profile of the stem and facilitates the estimation of merchantable logs distribution at any arbitrary diameter or length (Teshome 2011).
Taper functions are equations that describe the diameter narrowing rate along the stem (Gray 1956), which can be expressed as a function of the total height of the tree (H), and the diameter at the breast height (D  Clutter 1980). Taper equations utilize: (i) the diameter (with or without bark) at any point of the stem; (ii) the height at a certain diameter; and (iii) the commercial volume (with or without bark) at a given commercial height or minimum diameter from any stump height (Kozak 2004). According to Prodan et al. (1997), the use of those functions allows to assess the combination of wood products that can be obtained from a stand, and they are usually applied along with simulations of different silvicultural regimes. The integration of the taper function from ground level to any height will provide an estimate of the merchantable volume at that height (Kozak 1988, Bi 2000).
These functions can also be compatible with volume equations that are used to estimate individual tree volumes; ideally, the volume computed by integrating the taper equation from the ground to the tree top should be equal to that calculated using the total volume equation (Clutter 1980). Following DiéguezAranda et al. (2006), in the development of a compatible volume system, the most common approach is to express the α coefficient of a simplified Spurr’s volume equation (Spurr 1954) in terms of the taper equation coefficients, or vice versa, by using a compatibility relationship. This ensures that the taper function and the volume equation are analytically consistent (Sharma & Oderwald 2001). Another approach involves the development of a system that ensures compatibility between the taper function and a preestablished volume equation, by including this function in the mathematical expression of the taper function (Fang et al. 2000). Examples of volume equations used for this purpose are those proposed by Schumacher & Hall (1933) and Spurr (1954).
Compatible systems have been usually fitted using two different approaches: ordinary least squares (OLS) and simultaneously, using seemingly unrelated regression (SUR) or full information maximum likelihood (FIML  DiéguezAranda et al. 2006, RachidCasnati et al. 2014). In the classic OLS adjustment, it must be decided whether the error is minimized in the taper equations or in the merchantable volume. Then, the selected equation is fitted and the specific compatibility relationship is used to obtain the parameters of the nonfitted equation. On the contrary, the simultaneous fitting of the complete systems using FIML or SUR presents consistent estimators for all the equations of the system (Fang et al. 2000).
Data for adjusting taper models include multiple height and diameter measurements from each of the individuals in the sample, which determines a longitudinal data structure (Lindstrom & Bates 1990). Simultaneously, the variability of records within each tree is lower than the variability between trees; therefore, the assumption of independence is not fulfilled. Problems associated to the use of longitudinal data to taper functions fitting are related to multicollinearity, autocorrelation, and heteroscedasticity. In the presence of autocorrelation, the estimated parameters are not those of minimum variance since the mean square error of the model underestimates the variance of the error term, thus invalidating the significance contrasts. Although the estimators obtained in the regression adjustment remain unbiased, they are not the most efficient (Kozak 1988). For a correct analysis of longitudinal data, numerous studies propose to assume a particular error structure, expressing this autocorrelation among the errors as a stationary autoregressive model CAR(x) of order x (Garber & Maguire 2003). The advantages of this methodology are: (i) greater efficiency in obtaining parameters estimators of the model; (ii) adequate estimators for the standard errors of parameters; (iii) problems of lack of data or different number of measurements between individuals are solved; and (iv) it can be used even when the number of measurements is large compared to the number of individuals (Zimmerman & NúñezAntón 2001). According to Teshome (2011) these systems should be designed to be simple, flexible and easy to handle so that they can accurately estimate merchantable volume of any log length or to any specific diameter, based solely on easily measurable individual parameters such as diameter at breast height and tree height.
Tree growth follows biological limitations; stem growth depends on complex multifactorial and sitespecific factor combinations. Tree growth as well as stem shape depends on several intercorrelated stand specific factors, e.g., available resources, stand conditions, and structure attributes, which might potentially challenge the model’s predictive accuracy. This concern has been addressed in several studies by including stand density information in taper equations in order to improve model performance (Sharma & Zhang 2004, Sharma & Parton 2009). Calama & Montero (2006) pointed out that the high level of random variability with respect to stem form indicated the presence of factors that acted both at plot (silvicultural treatment) and tree levels (competition).
Jacobs et al. (2020) determined that when decreasing stand density, the lower part of the stem is conic instead of cylindrical. The lowering of stand density can cause a decrease in form factors, an increase in taper rates, and a decline in slenderness. Social status also affected stem shape: dominant trees increased in stem taper, while dominated trees in closed stands produced smaller crowns with less tapered stem (Jacobs et al. 2020).
To date, compatible systems of taper and volume equations are not widely used in Uruguay. The development of more accurate taper and volume models for intensive Eucalyptus plantations will significantly contribute to improve forest yield and planning. Therefore, the aim of the present study was to simultaneously fit a system of equations for predicting taper, merchantable volume and total tree volume, for Eucalyptus grandis (Hill ex. Maiden) and Eucalyptus dunnii (Maiden) growing in different zones of Uruguay. Specific objectives were: (i) to adjust four widely used stem taper profile models for predicting merchantable volume; (ii) to select and validate the best system; (iii) to assess the effect of stand density on stem taper and volume predictions. The present study involved simultaneous fitting, while addressing the autocorrelation, multicollinearity among independent variables, and heteroscedasticity.
Materials and methodsStudy site
Data was collected in the trials of the National Institute of Agricultural Research (INIA) in 2015 (PROBIO 2015), located in three different zones in Uruguay: northwest (32° 18′ 18″ S, 57° 44′ 06″ W, next to Guichon city, Paysandú); central east (33° 01′ 43″ S, 55° 31′ 08″ W, next to Sarandí del Yi, Durazno) and central west (33° 21′ 06″ S, 56° 41′ 17″ W, near Trinidad city, Flores  Fig. 1).
The study zones have a temperate subtropical climate, with a mean annual temperature of 18 °C (12 °C in the coldest month, 24 °C in the warmest month) and mean annual rainfall between 1300 and 1400 mm (Castaño et al. 2011). According to the classification of the National Commission for Agroeconomic Studies of the Land (CONEAT), Guichon predominant soils are planosols, with low fertility, horizon A 4050 cm depth, weak structure, low level of organic matter, slopes of 13%, sandy texture, medium to low risk of drought, imperfect drainage, moderatelyslow to slow permeability and good rooting ability. Trinidad soils are Vertisols, with medium fertility, horizon A 5060 cm depth, weak structure, medium level of organic matter, slopes of 25%, loam franc structure, low risk of drought, moderate drainage, moderately permeability and good rooting ability. Sarandí del Yi presented soils corresponding vertisols ruptic and Brunosols with a horizon A 5060 cm depth, sandyloam, loam franc texture, low fertility, low risk of erosion, moderate slopes (13%), weak structure, low organic matter, imperfect drainage and good rooting capacity.
Data collection
A total of 60 trees of Eucalyptus dunnii and 113 trees of E. grandis were felled and measured for diameters at breast height (D), total height (H), and diameters (d_{i}) at the corresponding heights (h_{i} = 0.2, 0.7, 1.3 m; and beyond this point, every 1.0 m up to total height  Tab. 1). The stump height was fixed at nearly 20 cm from the ground. Partial and total volumes of the stem (with and without bark), used for fitting the system, were estimated by Smalian method (Prodan et al. 1997). The top section volume was estimated using the volume formula for a cone. Plots of relative height (h_{i}/H) against relative diameter (d_{i}/D) for each species were carefully examined to detect possible anomalies, such as knots, bark defects, etc. (Fig. 2).
Compatible system fitting
Four compatible taper and volume systems widely used in forestry were selected from literature (Martin 1981, DiéguezAranda et al. 2006, Tang et al. 2016), to test and evaluate the taper and stem volume of both species (see Tab. S1 in Supplementary material). The compatible systems analysed had two components, a taper function and a total volume equation. In addition, the system proposed by Fang et al. (2000) includes also a merchantable volume equation to or above any specified height or diameter limit.
Model fitting
Compatible taper and volume equations were fitted using a simultaneous method, where all parameter must be expanded through a compatibility relationship. This method reduces the total system squared error (SAS Institute 2004), minimizing prediction errors in diameter at different heights and volume (DiéguezAranda et al. 2006). Parameter estimation was carried out using the “proc model” of SAS/ETS^{®} (SAS Institute 2004). Systems were simultaneously fitted using the FIML method, that properly accounts for crossequation correlation and is appropriate for systems of simultaneous equations (SAS Institute 2004). The full information of maximum likelihood estimation in the “proc model” procedure was considered. Simultaneous fitting of the systems applied requires an equal number of observations for all the equations. However, the number of observations in each equation was not equal, as there was more than one diameter observation for each tree but only one observation for total stem volume. Therefore, the structure of the data used (DiéguezAranda et al. 2006, Tang et al. 2016) includes the total volume and diameter of each tree and the inverse of the number of observations in each tree (1/n_{i}) used as a weight when fitting total volume.
Multicollinearity, autocorrelation and heteroscedasticity
Main challenges when fitting stem taper and volume equations lay in the violation of the fundamental least squares assumption of independence and equal distribution of errors with zero mean and constant variance. The most important issues are multicollinearity, autocorrelation and heteroscedasticity, which affect the efficiency of the regression coefficients while the leastsquares estimates remain unbiased and consistent.
Although the existence of multicollinearity does not seriously affect the predictive ability of the model (Kozak 1988), we evaluated the presence of multicollinearity between variables in the models analysed through the condition number (CN). According to Belsley (1991), when 5 ≤ CN ≤ 10, multicollinearity is not a serious issue. If 30 ≤ CN ≤ 100, problems associated with multicollinearity may arise, and if CN>100 severe problems of multicollinearity exist. Thus, for CN values greater than 30, corrections must be made.
The data used in the study includes hierarchical relationships, with multiple observations for each tree; thus, autocorrelation within the residuals for each individual might be expected. To demonstrate the autocorrelation, residuals were analysed graphically and by assessing the DurbinWatson statistic (Dw  HernándezRamos et al. 2017). Dw close to 2 allows the nonrejection of the null hypothesis (of nonexistence of autocorrelation); Dw lower than 2 indicates positive autocorrelation of order one, while values between 2 and 4 indicate negative autocorrelation of order one (Dufour & Dagenais 1985).
One of the general methods proposed to deal with continuous, unbalanced, multilevel longitudinal data is to model the correlation structure directly. In the present study we considered a secondorder continuous autoregressive error structure CAR (2), where h accounted for the distance between measurements, and their relative position on the stem was used to overcome the inherent autocorrelation of the longitudinal data. This structure was applied to data from the same sample unit measured at intervals that are not constant in time or space and with a different number of measurements in each sample unit, which is characteristic of taper functions. The CAR(x) error structure was programmed to enable dynamic updating of the residuals.
To account for secondorder autocorrelation, we used a CAR (2) model form that expands the error terms in the following way (Zimmerman & NúñezAntón 2001  eqn. 1):
where e_{ij} is the jth ordinary residual on the ith individual (i.e., the difference between the observed and the estimated diameters of the ith tree at jth height measurement), d_{1}=1 for j>k, k=1, 2 and it is zero for j≤k, ρ_{k} are the autoregressive parameters to be estimated, h_{ij}h_{ijk }are the distances separating the jth from the jkth observations within each tree, h_{ij}>h_{ijk}
When data from different sample units with different characteristics are used, a nonconstant the variance of the errors (σ_{ι}^{2}) is expected which leads to heteroscedasticity. To detect the deviation of the homoscedasticity of the residuals, the BreuschPagan test was applied. To remove the restriction of the variance heterogeneity of errors, since variances are unknown, an alternative is to graphically identify those independent variables (or combinations thereof) that cause such variation. Neter et al. (1996) suggested using the following formula (eqn. 2):
\sigma_ {i}^2 = \left (X_{i} \right )^k
Prodan et al. (1997), working with compatible taper and volume functions, proposed a power function of the variable D^{2}_{i}Ht_{i}, (the square of the tree diameter multiplied by the total height tree) as a candidate for (eqn. 3):
\sigma_{i} ^2 = \left (D_{i}^2Ht_{i} \right )^k
To estimate the kvalue (power term), we used the method suggested by Harvey (1976, in Parresol 1993), which consists of using the errors obtained with the unweighted model as a dependent variable in the power model of error variance.
The k parameter of eqn. 3 was estimated through linear regression using “proc reg” (SAS Institute 2012) following the methodology described above. According to DiéguezAranda et al. (2006), the weighting factor for heteroscedasticity 1/(D^{2}_{i} Ht_{i})^{k} along with the correction for the special structure of the data, was multiplied and programmed in the model procedure of SAS/ETS^{®} (SAS Institute 2012) by specifying the following (eqn. 4):
where n_{i} is the number of observations in each tree, and res are the residuals.
The use of models with random effects parameters can handle autocorrelation and heteroscedasticity, but increases operating costs (Liu et al. 2020). This methodology involves additional measurements to estimate random effects parameters, as upper stem diameterheight measurements (calibration) for improving diameter prediction (Liu et al. 2020).
Model comparison and validation
The accuracy and precision of each model were compared through graphical and numerical analysis of their residuals. The statistical criteria applied were: the adjusted coefficient of determination (R²adj, eqn. 5), coefficient of determination for nonlinear regression (R^{2}, eqn. 6), root of the mean squared error (RMSE, eqn. 7), and mean absolute bias (MAB, eqn. 8  Prodan et al. 1997). These expressions are summarized as follows:
where r^{2 }the correlation coefficient between the measured and estimated values (y_{i}, hat{y}_{i}, respectively), of the dependent variable, n is the total number of observations and p is the number of equation parameters.
The validation process was applied to assess whether the goodnessoffit also reflects in the quality of predictions (Myers 1990), which influences the selection of the best model. The use of cross validation in forestry is a common practice (Hirigoyen et al. 2018), and it is a method for model selection considering the predictive ability of the model. This is performed using the same data used for fitting (Myers 1990), and consists of the calculation of the residuals of the ith observation using parameters estimated by fitting the equation with all data except the ith observation. This process is known as leaveoneout crossvalidation, LOOCV (Kirschen et al. 2000).
The sum of squares of the eliminated residues is called PRESS (Predicted Residual Sum of Squares  Picard & Cook 1984), and it is used to calculate the selection criteria or RMSE for crossvalidation (RMSE_{cv}  eqn. 9, eqn. 10):
A close agreement between RMSE_{cv} and RMSE indicates that the model is not overfitting the data and has good predictive value (AndersonSprecher 1994). The Model Efficiency (ME) represents the proportion of variability observed in the original data that is explained by the model, and it varies between 0 (without adjustment) and 1 (perfect fit; Vanclay & Skovsgaard 1997  eqn. 11):
Plots of residuals versus predicted values, as well as versus dependent and independent variables, were examined for detecting bias. Statistics and the graphical analysis of residuals are essential to describe the goodnessoffit of a model and to understand which of the systems better represents observations, and therefore as criteria for recommending a model.
Model ranking
A classical procedure of ranking m models is to assign numbers (1, 2, 3, …, m) for model comparison. In this respective model order the exact place of a model with reference to other models is not known. Poudel & Cao (2013) proposed a method where relative ranks were developed to display the relative positions of the models. This method was used to obtain the specific and relative position of each model, where the relative rank of method i is defined as follows (eqn. 12):
where R_{i} is the relative rank of model i (i =1, 2, 3, …, m), S_{i} is the goodnessoffit statistics produced by model i, S_{min }is the minimum value of S_{i}, and S_{max }is the maximum value of S_{i}. Assigned numbers by respective order of the model show that the best and the worst methods have relative ranks of 1 and m. Because the magnitude (and not only the order) of the S_{i}’s are taken into consideration, this ranking method provides more information than the common ordinal ranks.
In our study the above ranking system was applied to R^{2}, RMSE and MAB statistics for diameter, merchantable and total volume to calculate the mean rank. The overall rank was determined by taking the mean of the ranks (MR) for these variables.
Stand density effects
The stand density variable was tested into the selected system taper equations to improve the fitting. Liu et al. (2020) in a study on Larix gmelinii showed the viability of including ecologically parameters (as stocking degree) and improving accuracy of taper equations by random effects.
Stand density was introduced through mathematical relationship between the resulting coefficients of volume equation and the density classes. For example, Schumacher & Hall (1933) volume model included in Fang et al. (2000) or its modified version in Sharma & Oderwald (2001) for a compatible system (v = a_{0 }dbh^{a1 }H^{a2}), has three parameters. Parameter a_{0} represents the estimated asymptotic volume at a determined dbh and H. Selecting parameter a_{0} as expanded parameter in the model can improve the estimation and indicates that the stand density has a significant effect on dbh and H individual volume relationship. The suitable taper system with stand density variable was fitted and evaluated by the whole fit data set and cross validation.
We adjusted model compatible system taper and volume for E. grandis and E. dunnii with stand density into parameters used in the volume compatible equations (Tab. 2). The resultant parameter prediction equation for predicting a_{0} can be given by (eqn. 13 to eqn. 16):
where N is the stand density, and a_{01} and a_{02} are the parameters to be estimated. Equations (13 to 16) were substituted into selected model system to examine the effect of stand density on tree taper and volume. We named this modified version of the model as “Modified system”, to differentiate it from the original (“Original system”).
ResultsModel fitting and selection
Parameter shared by both the taper and volume equations estimates and their approximated standard errors for each compatible volume system are listed in Tab. 2. The goodnessoffit statistics for both species are shown in Tab. 3. All systems were initially fitted without expanding the error term to account for autocorrelation. As expected, because of the hierarchical nature of the data, the values of the DurbinWatson (Dw) test for all functions of both species indicated that the residual trend could be described as a function of lag1 and lag2residuals within the same tree. The error term was, therefore, expanded as a first order and a second order autoregressive structure. The correlation trend disappears with a firstorder autoregressive error structure which was sufficient to eliminate the autocorrelation in the system for both species (Tab. 3).
Considering goodnessoffit statistics and MR, the best models were Max & Burkhart (1976) and Fang et al. (2000) for both species. These equations exhibited the lowest values of RMSE (<0.8 cm), MAB (<0.62 cm), and the highest R^{2}adj (98.8%) for E. grandis and E. dunnii.
Concerning the total stem volume, the MR suggests that the models fall into two groups separated by a large gap. Fang et al. (2000), Max & Burkhart (1976) and Kozak et al. (1969) models were included into a group with lower MR for E. grandis (MR=1.1) and E. dunnii (MR=1.8). Conversely, Sharma & Oderwald (2001) appeared to be inadequate for both species (MR=4). R^{2}adj, RMSE and MAB show the same trend: all models except Sharma & Oderwald (2001) were competitive and small differences were found. The total volume variance was similar for all the models for both species. However, the Fang et al. (2000) model had the best MR for both species and was therefore selected. Based on the best fit, 99% of the total variation in predicting upper stem diameters and total volume was explained; RMSE was lower than 1 cm and 0.035 m^{3} in predicting stem diameter and total volume for both species. The model had an average absolute bias less than 0.0025 m^{3 }for volume over bark and 0.61 cm for diameters.
In addition, the Fang et al. (2000) model allowed the inclusion of a function to estimate merchantable volume; thereby, we readjusted a system with three equations: taper, merchantable volume (integration of taper equation to any height limit), and total stem volume (integration of taper equation to total tree height). The coefficient estimates for complete system after correction for autocorrelation, their approximated standard errors, and the goodnessoffit statistics of the taper function, merchantable stem volume and total stem volume are given in Tab. 4.
For Fang et al. (2000) with three equation, mean average bias (MAB) were 0.610 and 0.659 cm for taper equation; 0.003 m^{3} and 0.016 m^{3 }for total volume and 0.031 m^{3 }and 0.022 m^{3 }for merchantable volume for E. grandis and E. dunnii, respectively. Scatter plots of studentized residuals against predicted and observed values against predicted values for d_{ub}, merchantable volume and total volume are shown in Fig. 3. Residuals did not show a trend of increasing error variances, and the zerostudentized residuals cross the centre of the data points. This suggests that the taper and volume function were appropriately identified, and the error structure of the model is associated with the equal error variance. The multicollinearity was inferred by CN, 46 for E. grandis and 56 for E. dunnii.
The precision of predictions obtained by the Fang et al. (2000) model through the LOOCV procedure was evaluated in terms of model efficiency (ME) and root mean square error (RMSE_{cv}) of estimated diameters, merchantable volume, and total stem volume. For E. grandis, ME was 0.98, 0.98 and 0.87, and RMSE_{cv }was 0.803 cm, 0.037 m^{3}, 0.024 m^{3} for taper, total volume and merchantable volume, respectively. For E. dunnii, ME was 0.95, 0.98 and 0.87, and RMSE_{cv} was 1.03 cm, 0.028 m^{3}, 0.03 m^{3} for taper, total volume and merchantable volume, respectively. Mean average bias (MAB) values were 0.708 and 0.877 cm for taper equation, 0.026 m^{3 }and 0.025 m^{3 }for total volume, and 0.014 m^{3 }and 0.019 m^{3 }for merchantable volume for E. grandis and E. dunnii, respectively.
Model development and evaluation with density variable
When comparing the quotient between RMSE_{cv} and RMSE of the selected model, we expected values lower or close to 1, indicating good predictive value of fitted models. For E. dunnii this was confirmed since the RMSE ratio was always lower than 1 (range: 0.800.88) for taper, total and merchantable volume. However, we found quotient values of 1.14 for total volume and 1.7 for merchantable volume for E. grandis, which indicates a possible overfitting.
The parameters p_{1} and p_{2} of Fang et al. (2000) segmented model are the inflection points and marked differences between stems. To check the usefulness of an adjustment by species, the values of these parameters were compared with those obtained in a general model. These models were compared through R^{2}adj, RMSE and AIC values.
The calculated parameter p_{1} was 0.037, 0.031, 0.039, while p_{2} was 0.868, 0.941, 0.908 for E. dunnii, E. grandis and both species together, respectively. We obtained values of R^{2}adj = 0.98, 0.98, 0.96; RMSE = 0.8867 cm, 0.712 cm, 1.235 cm; ACI = 6099, 6093, 7269 for E. dunnii, E. grandis and both species together, justifying the adjustment by species.
The stand density range significantly affected the compatible system for E. grandis. The stand density in this study was 390, 453, 535, 850, 930 and 1090 tree ha^{1}. However, for E. dunnii there were not significant effects, probably because the study only covered a limited range of stand densities (930, 990, 1010 and 1100 tree ha^{1}) for this species.
Equations (13 to 16) were fitted and evaluated to the entire data set composed of different stand density by MR (calculated on R^{2}adj, RMSE and MAB values). Eqn. 13 obtained the first place into mean ranking (R^{2}adj = 0.98, 0.99, 0.99; RMSE = 0.895 cm, 0.034 m^{3}, 0.031 m^{3}; MAB = 0.668 cm, 0.024 m^{3}, 0.025 m^{3 }for taper, total and merchantable volume, respectively), followed by eqn. 14. Regarding bias, eqn. 13 had the lowest values: 0.18 cm, 0.0025 m^{3} and 0.021 m^{3}.
In reference to the Original system, Modified system improves in terms of R^{2}adj, RMSE and MAB for total and merchantable volume, while for the taper equation the goodnessoffit was almost invariant (Tab. 5). The effect of stand density was analysed visually by generating tree profiles, total and merchantable volume using Original and Modified systems for an average tree with d = 22.7 cm and H = 26.6 m at different stand densities (Fig. 4). Box plots of residuals for both systems, Modified and Original, for total volume and diameter against relative height class and relative diameters class are reported in Fig. 5.
Data for total volume was grouped into four diameter classes (20, 25, 30, and 35 cm). Original system has an overestimation trend while Modified system improved results for all classes, particularly at high diameters. For diameter data, height classes were 18, 25, 35 and 45 m, both Original and Modified system showed an overestimation for diameters in all classes. However Modified system residuals have a homogeneous distribution around zero.
DiscussionOriginal system model
After estimating the fit and crossvalidation statistics of the four fitted models, the Fang et al. (2000) model was selected and used in subsequent analysis. To ensure numeric consistency, a simultaneous fitting procedure based on the FIML estimation method was used. This method has been successfully used in several studies because it optimizes parameters of taper and commercial volume while it minimizes and homogenizes standard deviation. Results for the taper and volume equations fitted for both species indicated that Fang et al. (2000) is the most precise and least biased model. This model had the first place in the model ranking. Previous research found that the modification of the Muhairwe variable exponent model was the best model to describe stem shape to E. grandis (Gomat et al. 2011). Alternatively, Methol (2008) adjusted the segmented model of Max & Burkhart (1976) for both species. This model came second in our study.
The system of Fang et al. (2000) is widely used for the description of stems form and estimated volume, and its implementation has been reported in biometric systems for a wide range of species and countries (DiéguezAranda et al. 2006, Tamarit Urías et al. 2014, Tang et al. 2016, HernándezRamos et al. 2017). This model has the merits of being flexible and analytically integrable into a system providing high accurate estimates of stem volume for any stem segment. RachidCasnati et al. (2014) adjusted for E. grandis the explicit model for total volume of SchumacherHall, and this model is integrated to the adjusted model in this study. The quality of fit in terms of ME was similar in both studies, but the RSME was slightly lower in RachidCasnati et al. (2014), and the parameters values were similar, with small differences attributable to the simultaneous adjustment of the models.
Fang et al. (2000) segmented model presents two inflection points on the relative height of the stem: p_{1} and p_{2} are the relative heights from the ground where the points are assumed by the model to occur, dividing the stem in three sections. The first inflection point (as percentage of the total height), where the neiloid changes to a paraboloid, and the second inflection point, at which the change from the paraboloid to a cone occurs. Inflection points for the selecting system occur at 3.1% and 94% for E. grandis and 3.7% and 86% for E. dunnii. The first inflection point were a very small percentage of total height in both species. These results are similar to those reported by previous studies, e.g., for E. grandis p_{1} was 3.6% in Ozçelik & Göçeri (2015), although p_{2} was lower (45%) than that observed in this study. The presence of multicollinearity (correlations among variables) in the model has been examined by CN, its value for both species was in range of 30100. This values are tolerated when models contain polynomial and crossproduct terms (Belsley 1991). The equation systems of Fang et al. (2000) showed weak multicollinearity for both Eucalyptus species, and Mora et al. (2013) stated that multicollinearity does not seriously affect the predictive ability of the model. By comparing RSME_{cv} through LOOCV, and RMSE of the simultaneous adjustment, it was evidenced that the values were very close (Tab. 5), and this indicates that the model is not overfitting the data and has a good predictive performance.
Modified system model
Sharma & Parton (2009) improved fit statistics and predictive precision by including stand density information into a variable taper equation. These studies identified tree competition as a relevant stand characteristic responsible for volume miscalculation. Liu et al. (2020) assessed taper models and included both fixedeffects parameters, showing that age and stocking degree had significant effects on tree stem and improving the predictive taper equations. Traditional and simple volume functions do not consider the dependency of volume on competition or stem form. They generally involve dbh and H as predictors (Spurr 1954).
Including age and stand density in models, part of the site variability that affects the growing conditions of each individual tree is taken into account. This has proven to enhance predictions of the stem taper (Liu et al. 2020).
In our study the stand density had a statistically significant influence in volume estimation for E. grandis, while for E. dunnii the limited range of stand density compared to E. grandis could be too narrow to properly explore effects on volume.
Jacobs et al. (2020) pointed out that the effect of competition was significant on both stem volume and stem taper in Picea abies. Sharma & Parton (2009) reported that the difference in bole diameter between trees at lower and higher stand densities diminished increasing stand density, the density treatments only significantly affected the stem form of the smallest trees. According to Jacobs et al. (2020), there is a lack of knowledge on how the stem form of trees with the same dbh and H is affected by the individual competition situation inside the stand. However, for the average tree no differences in stem profile were observed, while for merchantable and total volume the differences were considerable.
Results of the Original Fang et al. (2000) model for merchantable volume showed an underestimation for densities greater than 400 tree h^{1}. The effect is stronger in the tree top section (0.8 to 1.0 relative height). For the total volume the trend in underestimation for all relative height classes is more evident. This means that the Original Fang et al. (2000) model tends to underestimate the total volume regardless of the total height of the tree, while the underestimation of the merchantable volume is limited to the upper part of the tree.
Our results did not show difference in stem form by density range. Similar results were found by Ahnlund et al. (2014) with Pinus sylvestris, where there were no significant differences in stem form for trees with dbh > 5 cm between density treatments. In a study on Eucalyptus in Congo with density range of 5001300 tree h^{1}, Gomat et al. (2011) reports that planting density had no significant influence on the stem profile. Unlike the findings of Jacobs et al. (2020), who reports density effects in stem form, the butt logs of trees of the same size decreased in taper given increasing competition, and stem taper of the top log decreased with increasing competition, while the stem taper of the mid log increased with competition.
The inflection points for Modified system do not vary significantly compared to the Original system. These results show that a segmented model with two inflection points is appropriate to describe the stems of Eucalyptus species, and that stand density has no effects. Gomat et al. (2011) studied factors affecting stem taper variation in Eucalyptus, reporting that the tree shape markedly changes between 1 and 2 years and is roughly stabilised afterwards. The plantations in this study had range between 717 years, therefore pronounced changes are not expected in the inflection points.
The overall form of the trees remains the same irrespective of the stand age and density, although slight volume differences can be observed (Fig. 5).
Results showed that stand density had an effect on volume parameters. Based on the fitting, validation statistics, and graphical analysis, the Fang et al. (2000) system with stand density information is recommended for E. grandis. Thus, a system including stand density results in more accurate estimates for E. grandis, as stand density affects the volume and shape of stems. However, for E. dunnii, the inclusion of stand density was not significant and the variable was not included in compatible volume and taper systems.
Conclusions
The compatible volume system of Fang et al. (2000) is the most adequate for describing the stem profile and predicting stem volume of E. dunnii and E. grandis. A modified second order continuous autoregressive error structure corrected the inherent serial autocorrelation of different observations within trees. In E. grandis, the model fit and the prediction performance of the taper and volume compatible systems were greatly improved by including stand density within the volume equation.
For E. dunnii no improvement was observed probably due to a narrow population range in the dataset. The inclusion of stand density in the model improved the efficacy of the system for E. grandis, while for E. dunnii no significant effects of stand density were found. Including stand density in taper and volume equations has practical implications, since for the same diameter and height a range of values are obtained depending on the population of the stand. This improves the characterization of stands associated to different management practices and subsequently assortment of products.
Acknowledgments
The authors thank the Instituto Nacional de Investigaciones Agropecuarias (INIAUruguay) for supporting fieldwork and the INIA Scholarship for PhD studies. We are particularly grateful for the support of Leonidas Carrasco Letelier, Roberto Scoz, Demian Gomez, Gonzalo Martinez, Jorge Basso (INIA) and Juan Gabriel Álvarez González (USC). We acknowledge the institutional support of the University of Cordoba  Campus de Excelencia CEIA3. We also thank the ERSAF group. We thank Dr. David Walker for his revisions of the manuscript versions, and the anonymous referees for their comments and corrections.
ReferencesAhnlund K, Ulvcrona T, Nilsson U, Lundmark TStand density and fertilization effects on aboveground allocation patterns and stem form of Pinus sylvestris in young stands. Scandinavian Journal of Forest Research 29: 197209.2014AndersonSprecher RModel comparisons and r^{2}. American Statistician 48: 113117.1994Belsley DAConditioning diagnostics: collinearity and weak data in regression. Operational Research Society 44: 8888.1991Bi HTrigonometric variableform taper equations for Australian eucalypts. Forest Science 46: 397409.2000Calama R, Montero GStand and treelevel variability on stem form and tree volume in Pinus pinea L.: a multilevel random components approach. Investigación Agraria: Sistemas y Recursos Forestales 1 (5): 2441.2006Castaño JP, Giménez A, Ceroni M, Furest J, Aunchayna R, Bidegain MCaracterización agroclimática del Uruguay 19802009. [Agroclimatic characterization of Uruguay 19802000]. Serie Técnica INIA 193. Montevideo, Uruguay, pp. 33. [in Spanish]2011Clutter JDevelopment of taper functions from variabletop merchantable volume equations. Forest Science 26: 117120.1980DiéguezAranda U, CastedoDorado F, AlvarezGonzález JG, Rojo ACompatible taper function for Scots pine plantations in northwestern Spain. Canadian Journal of Forest Research 36: 11901205.2006Dufour JM, Dagenais MGDurbinWatson tests for serial correlation in regressions with missing observations. Journal of Econometrics 27: 371381.1985Fang Z, Borders BE, Bailey RLCompatible volumetaper models for loblolly and slash pine based on a system with segmentedstem form factors. Forest Science 46: 112.2000Garber SM, Maguire DAModeling stem taper of three central Oregon species using nonlinear mixed effects models and autoregressive error structures. Forest Ecology and Management 179: 507522.2003Gomat HY, Deleporte P, Moukini R, Mialounguila G, Ognouabi N, Saya AR, Vigneron P, SaintAndre LWhat factors influence the stem taper of Eucalyptus: growth, environmental conditions, or genetics? Annals of Forest Science 68 (1): 109120.2011Gray HRThe form and taper of foresttree stems. Imperial Forestry Institute, University of Oxford, UK, pp. 79.1956HernándezRamos J, HernándezRamos A, GarcíaMagaña JDJ, GarcíaCuevas X, GarcíaEspinoza GG, Muñoz Flores HJ, OlveraDelgadillo EHSistema compatible de ahusamientovolumen comercial para plantaciones de Pinus greggii Engelm. en Hidalgo, México [Compatible tree taper and volume system for commercial plantations of Pinus greggii Engelm in Hidalgo, Mexico]. Revista Mexicana de Ciencias Forestales 8: 5970. [in Spanish]2017Hirigoyen A, Franco J, DiéguezAranda UModelo dinámico de rodal para Eucalyptus globulus (L.) en Uruguay [Dynamic stand model for Eucalyptus globulus (L.) in Uruguay]. Agrociencia Uruguay 22: 6380. [in Spanish]2018Jacobs M, Rais A, Pretzsch HAnalysis of stand density effects on the stem form of Norway spruce trees and volume miscalculation by traditional form factor equations using terrestrial laser scanning (TLS). Canadian Journal of Forest Research 50: 5164.2020Kirschen RH, O’Higgins EA, Lee RTThe Royal London space planning: an integration of space analysis and treatment planning. American Journal of Orthodontics and Dentofacial Orthopedics 118: 456461.2000Kozak A, Munro DD, Smith JHGTaper functions and their application in forest inventory. The Forestry Chronicle 45: 278283.1969KozakA variableexponent taper equation. Canadian Journal of Forest Research 18: 13631368.1988Kozak AMy last words on taper equations. Forestry Chronicle 80: 507515.2004Lindstrom MJ, Bates DMNonlinear mixed effects models for repeated measures data. Biometrics 46 (3): 673.1990Liu Y, Yue C, Wei X, Blanco JA, Trancoso RTree profile equations are significantly improved when adding tree age and stocking degree: an example for Larix gmelinii in the Greater Khingan Mountains of Inner Mongolia, northeast China. European Journal of Forest Research 139: 443458.2020Martin AJTaper and volume equations for selected Appalachian hardwood species. Research Paper NE490, USDA Forest Service, Northeastern Forest Experiment Station, Broomall, PA, USA, pp. 22.1981Max T, Burkhart HESegmented polynomial regression applied to taper equations. Forest Science 22: 283289.1976Methol RSAG Eucalyptus: sistema de apoyo a la gestión de Eucalyptus orientadas a la producción de celulosa en Uruguay [SAG Eucalyptus: Eucalyptus management support system aimed at the production of cellulose in Uruguay]. INIA Serie Técnica 173, Montevideo, Uruguay, pp. 26. [in Spanish]2008Mora B, Wulder MA, White JC, Hobart GModeling stand height, volume, and biomass from very high spatial resolution satellite imagery and samples of airborne LIDAR. Remote Sensing 5: 23082326.2013Myers RHClassical and modern regression with applications. Duxbury press, Belmont, CA, USA, vol. 2, pp. 473.1990Neter J, Kutner MH, Nachtsheim CJ, Wasserman WApplied linear statistical models. Irwin, Chicago, IL, USA, pp. 318.1996Ozçelik R, Göçeri MFCompatible merchantable stem volume and taper equations for eucalyptus plantations in the eastern Mediterranean region of Turkey. Turkish Journal of Agriculture and Forestry 39: 851863.2015Parresol BRModeling multiplicative error variance  An example predicting tree diameter from stump dimensions in bald cypress. Forest Science 39: 670679.1993Picard RR, Cook RDCrossvalidation of regression models. Journal of the American Statistical Association 79 (387): 575583.1984Poudel KP, Cao QVEvaluation of methods to predict Weibull parameters for characterizing diameter distributions. Forest Science 59: 243252.2013PROBIOMejoramiento en la calidad de la información vinculada con la utilización de la biomasa forestal [Improvement in the quality of information related to the use of forest biomass]. Producción de Electricidad a partir de Biomasa  PROBIO, Montevideo, Uruguay, pp. 157. [in Spanish]2015Prodan M, Peters R, Cox F, Real PMensura Forestal [Forest Measurement]. GTZ/IICA, Serie Investigación y Educación en Desarrollo Sostenible, San José, Costa Rica. pp. 561. [in Spanish]1997RachidCasnati C, Mason E, Woollons R, Resquin FVolume and taper equations for P. taeda (L.) and E. grandis (Hill ex. Maiden). Agrociencia Uruguay 18: 4760.2014SAS InstituteSAS/ETS 9.1 user’s guide. SAS Institute Inc., Cary, NY, USA, pp. 4420.2004SAS InstituteSAS/STAT user’s guide, version 9.1. SAS Institute Inc., Cary, NY, USA, pp. 128.2012Schumacher F, Hall FLogaritmic expression of timbertree volume. Journal of Agricultural Research 47: 719734.1933Sharma M, Oderwald RGDimensionally compatible volume and taper equations. Canadian Journal of Forest Research 31: 797803.2001Sharma M, Zhang SYVariableexponent taper equations for jack pine, black spruce, and balsam fir in eastern Canada. Forest Ecology and Management 198 (13): 3953.2004Sharma M, Parton JModeling stand density effects on taper for jack pine and black spruce plantations using dimensional analysis. Forest Science 55: 268282.2009Spurr SSimplified computation of volume and growth. Journal of Forestry 52: 914922.1954Tamarit Urías JC, De Los Santos Posadas HM, Aldrete A, Valdez Lazalde JR, Ramírez Maldonado H, Guerra De La Cruz VSistema de cubicación para árboles individuales de Tectona grandis L.f. mediante funciones compatibles de ahusamientovolumen [Volume estimation system for individual Tectona grandis L.f. trees through compatible taper/volume functions]. Revista Mexicana de Ciencias Forestales 5: 5874. [in Spanish]2014Tang X, PérezCruzado C, Fehrmann L, AlvarezGonzález JG, Lu Y, Kleinn CDevelopment of a compatible taper function and standlevel merchantable volume model for Chinese fir plantations. PLoS One 11 (1): e0147610.2016Teshome TCompatible volumetaper equations for predicting merchantable volume to variable merchantable limits for Cupressus lusitanica, Ethiopia. SINET: Ethiopian Journal of Science 28: 1522.2011Vanclay JK, Skovsgaard PEvaluating forest growth models. Ecological Modelling 98: 112.1997Zimmerman DL, NúñezAntón VParametric modelling of growth curve data: an overview. Test 10 (1): 173.2001
Location of the study area (yellow points) and soils forest classification sites in Uruguay.
Relative diameter (quotient between an objective diameter and the diameter at breast height) against relative height (quotient between height to an objective diameter and total tree height) for Eucalyptus grandis (a) and E. dunnii (b) in Uruguay.
Studentized residuals and observed values against predicted values for: (a) diameters under bark (d_{ub}, cm); (b) merchantable volume (m^{3}); and (c) total volume (m^{3}) obtained with the Fang et al. (2000) model for E. dunnii (firsts row) and E. grandis (seconds row) in Uruguay. (R^{2}): coefficient of determination.
Profiles (A), total volume (B) and merchantable volume (C) generated by Original (points) and Modified systems, for the average tree (d = 22.7 cm, H = 26.6 m) at different stand densities (blue line: 400 trees ha^{1}; green line: 850 trees ha^{1}; red line: 1200 trees ha^{1}; black line: 1600 trees ha^{1}) for Eucalyptus grandis in Uruguay.
Residuals boxplot of estimated total volume and estimated diameters by relative diameters classes of the Fang et al. (2000) Original and Modified systems. The boxes represent the interquartile ranges. The maximum and minimum total volume and diameters prediction errors are represented respectively by the upper and lower whiskers.
Summary statistics of the data for taper and volume of Eucalyptus grandis and E. dunnii in Uruguay. Variables and abbreviations: diameter over bark (d_{ob}, cm); diameter under bark (d_{ubi}, cm); diameter at breast height (D, cm); total height (H, m); total tree volume over bark (V_{ob}, m^{3}); total tree volume under bark (V_{ub}, m^{3}). (N): number of trees; (Sd): standard deviation; (Min): minimum; (Max): maximum.
Species
Variable
N
Mean
Sd
Max
Min
E. dunnii
D (cm)
60
18.21
3.55
24.40
9.90
H (m)
60
23.72
3.21
30.75
15.72
V_{ob }(m^{3})
60
0.28
0.13
0.65
0.05
V_{ub }(m^{3})
60
0.24
0.11
0.55
0.04
d_{obi }(cm)
1285
12.33
5.78
33.10
0.00
d_{ubi }(cm)
1285
11.32
5.31
31.83
0.00
E. grandis
D (cm)
113
22.71
5.82
36.90
10.95
H (m)
113
27.56
5.72
40.90
16.60
V_{ob }(cm^{3})
113
0.56
0.39
1.83
0.06
V_{ub }(cm^{3})
113
0.51
0.36
1.69
0.06
d_{obi }(cm)
2667
15.64
7.23
38.45
0.00
d_{ubi }(cm)
2667
14.80
6.95
33.50
0.00
Values of the estimated parameters for simultaneous fitting for Eucalyptus grandis and Eucalyptus dunnii. Approximated standard errors are in brackets. All parameters significant at pvalue < 0.0001.
Species
Volume system
a_{0}
a_{1}
a_{2}
b_{1}
b_{2}
b_{3}
b_{4}
p_{1}
p_{2}
E. grandis
Fang et al. (2000)
0.00005 (3.34E7)
2.136 (0.0014)
0.761 (0.0014)
3.88E6 (3.69E8)
0.00003 (1.03E7)
0.0004 (0.0002)

0.027 (0.0002)
0.946 (0.0054)
Kozak et al. (1969)



2.9181 (0.008)
1.0254 (0.004)




Max & Burkhart (1976)

0.694 (0.026)
0.038 (0.0006)
2.667 (0.111)
1.213 (0.062)
0.887 (0.073)
524.4 (15.31)


Sharma & Oderwald (2001)
0.00003 (7.58E9)
2.207 (0.0003)







Fang et al. (2000)
0.00002 (8.64E7)
1.860 (0.0063)
1.194 (0.0095)
4.20E6 (4.91E8)
0.00003 (1.55E7)
0.00006 (8.44E6)

0.034 (0.0004)
0.851 (0.006)
E. dunnii
Kozak et al. (1969)



2.431 (0.023)
0.789 (0.012)




Max & Burkhart (1976)

0.701 (0.104)
0.0413 (0.001)
2.261 (0.236)
1.002 (0.133)
0.410 (0.146)
503.38 (25.86)


Sharma & Oderwald (2001)

2.243 (0.0004)
0.00003 (4.67E8)






Goodnessoffit statistics of compatible volume systems for Eucalyptus grandis and E. dunnii in Uruguay. (R^{2}adj): adjusted coefficient of determination; (RMSE): root mean squared error; (MAB): mean absolute bias; (Dw): DurbinWatson statistic; (MR): mean rank.
Species
System
Taper equation
Dw
Volume equation
RMSE (cm)
R^{2}adj
MAB (cm)
MR
RMSE (m^{3})
R^{2}adj
MAB (m^{3})
Dw
MR
E. grandis
Fang et al. (2000)
0.790
0.99
0.610
2.2
1.000
0.035
0.99
0.024
1.99
1.000
Kozak et al. (1969)
1.789
0.93
1.255
2.08
2.000
0.035
0.99
0.024
1.98
1.100
Max & Burkhart (1976)
0.840
0.98
0.620
2.2
1.002
0.035
0.99
0.024
1.98
1.075
Sharma & Oderwald (2001)
1.146
0.97
0.816
1.9
1.382
0.036
0.99
0.025
1.96
4.000
E. dunnii
Fang et al. (2000)
0.787
0.99
0.612
2.0
1.000
0.021
0.98
0.015
1.9
1.000
Kozak et al. (1969)
1.802
0.88
1.352
2.0
2.000
0.022
0.97
0.015
2.09
1.789
Max & Burkhart (1976)
0.898
0.97
0.624
2.2
1.001
0.023
0.97
0.015
2.0
1.711
Sharma & Oderwald (2001)
1.042
0.96
0.800
1.9
1.253
0.024
0.97
0.018
2.1
4.000
Values of the estimated parameters for simultaneous fitting of Fang et al. (2000) model. Approximated standard errors are in brackets. Goodnessoffit statistics for taper, total volume and merchantable volume for Eucalyptus grandis and E. dunnii in Uruguay are reported. (R^{2}adj): adjusted coefficient of determination; (RMSE): root mean squared error.
Equation
Parameters
Estimates
E. grandis
E. dunnii

a_{0}
0.00004(2.4E7)
0.00002(5.4E7)
a_{1}
2.09(0.001)
1.87(0.005)
a_{2}
0.862(0.001)
1.233 (0.007)
b_{1}
4.2E06(3.3E5)
4.3E6(4.5E8)
b_{2}
0.00003(5.2E5)
0.00003(1.3E7)
b_{3}
0.0002(0.00008)
0.00018(0.00002)
p_{1}
0.031(0.0002)
0.037(0.0004)
p_{2}
0.941(0.011)
0.868(0.0033)
Taper
RMSE (cm)
0.71
0.87
R^{2}adj
0.98
0.98
Total Volume
RMSE (m^{3})
0.04
0.02
R^{2}adj
0.99
0.98
Merchantable Volume
RMSE (m^{3})
0.04
0.03
R^{2}adj
0.98
0.97
Values of the estimated parameters for the Modified model (Fang et al. 2000), approximated standard errors (Estd), and goodnessoffit statistics for taper, total volume and merchantable volume for Eucalyptus grandis in Uruguay. (R^{2}adj): adjusted coefficient of determination; (RMSE): root mean squared error; (MAB): mean absolute bias.
Equation
Parameters
Estimate
Estd
pvalue

p_{1}
0.032
0.00026
<0.0001
p_{2}
0.942
0.00959
<0.0001
a_{01}
0.000023
3.0E7
<0.0001
a_{02}
0.162
0.0045
<0.0001
a_{1}
2.079
0.0015
<0.0001
a_{2}
0.953
0.0019
<0.0001
b_{1}
4.26E06
3.4E8
<0.0001
b_{2}
0.000032
8.0E8
<0.0001
b_{3}
0.000401
1.5E4
0.0086
Taper
RMSE (cm)

0.795

R^{2}adj

0.98

MAB (cm)

0.668

Total volume
RMSE (m^{3})

0.034

R^{2}adj

0.99

MAB (cm)

0.024

Merchantable volume
RMSE (m^{3})

0.031

R^{2}adj

0.99

MAB (cm)

0.025

Tab. S1  Compatible taper and volume systems tested for Eucalyptus grandis and Eucalyptus dunnii plantations in Uruguay.