Comparison of stem taper equations for eight major tree species in the Spanish Plateau

Aim of study: A stem taper function and a compatible merchantable volume system are compared to evaluate which provides a better description of the stem profile for the main species in central Spain. Area of study: This research was carried out in the region of Castile-Leon, located in Central Spain. Material and Methods: A total of 6,357 trees were selected for destructive sampling. All models were fitted using a first-order continuous autoregressive error structure to address the problem of autocorrelation. Main results: In terms of accuracy, the root mean square error (RMSE) in both models ranged from 0.75 to 2.72 depending on the species analyzed, presenting values similar to those reported in other studies. Small differences in the goodness-of-fit for both procedures were also found, and the Stud model provided better accuracy for 6 of the 8 species studied, with RMSE reductions of 0.5% to 8.6%. The RMSE obtained in the cross-validation phase was on average 1.22 times higher than what was obtained in the fitting phase. Research highlights: The non-linear extra sum of squares method indicated that the stem taper differs among the five softwood species and three hardwood species. In hardwoods, the first inflection point is lower than in softwoods (at around 5%) and the second inflection point is higher (at around 85%) than those of softwoods.


Introduction
Accurate predictions for wood products classified by merchantable size are a matter of interest for forest managers and forestry companies, in order to estimate the monetary value of some of the many commodities and services that forests provide to society. The accuracy of this estimate is directly related to its final use. Thus, while forest managers need it for planning and quantification of land use, forestry companies need it to assess the profitability of a harvest.
Nowadays, society demands multifunctionality from our forests, which increases the need for more detailed knowledge of the different products and services a forest can provide. Product classification thus becomes an important tool for assessing the role of forest and wood products, especially in light of climate changes due to carbon fixation and sequestration. Volume pre-2

Study area and data collection
This research was carried out in the region of Castile-Leon, located in Central Spain. The region covers approximately 9.4 million ha and is one of the most important areas for timber production in Spain. Predominant oak stands cover more than half of the area, while pine stands cover one third of the area. Altitude fluctuates between 110 and 2650 m above sea level. The climate in Castile-Leon is both continental and Mediterranean: average winter temperatures range between 4 and 7 °C, while average summer temperatures range from 19 to 22 °C; with three or four dry summer months that are typical of a Mediterranean climate. Average annual rainfall is only 450 to 500 mm, mostly in lower altitudes.
The data used in this study were collected in 242 public and private forests in Central Spain. Data from 1,844 Scots pines, 456 stone pines, 533 black pines, 1,715 Mediterranean maritime pines, 326 Spanish junipers, 302 Pyrenean oaks, 992 poplars and 189 beeches were used to test the statistical performance of the taper functions. Thus, a total of 6,357 trees were selected for destructive sampling using the protocol of Garber and Maguire (2003). The trees were felled from thinned and unthinned stands, which were subjectively selected to represent the existing range of site qualities. Before felling, two attributes were recorded for each sample tree: (i) diameter at breast height, D (to the nearest 0.1 cm); (ii) total tree height, H (to the nearest 0.01 m). Each sample tree was felled in a manner that minimized stem breakage. A tape measure was stretched along the bole and total height was recorded from the base of the stump to the top of the tree. The tree was divided into sections, and thin disks were removed for diameter measurements at the base of the stump; at 80 cm above the ground level; at breast height; and at a height midway between every other whorl pair above breast height (1-3 m intervals). For each disk, height (h, to the nearest 0.01 m) and diameter over bark (d, to the nearest 0.1 cm) were measured along two perpendicular axes. Log volumes in cubic meters were calculated using Smalian's formula. The topmost log was considered conical in shape. Total volume (V in m 3 ) was obtained by summing the log volumes. Between 3 and 40 disks were cut per tree, for a total of 87,568 disks. Summary statistics including the number of observations, arithmetic mean, standard deviation, and minimum/maximum values of the main variables are presented in Table 1. merchantable volume equations derived from volumeratio equations are very easy to use and develop, those obtained from taper functions are preferred nowadays; perhaps because they allow for estimation of diameter at a given height (Diéguez-Aranda et al., 2006).
Compatibility means that an integrated model can be obtained through summation of the differential model. Thus, for a given merchantable volume equation, there is an intrinsically defined compatible taper function (Clutter, 1980). This implies that integration of the taper function from the ground to the top of the tree would provide the appropriate volume, and subsequently the merchantable volume equation (Fang et al., 2000). Existing research indicates that segmented models appear to be more accurate than other model formulations (e.g., Corral-Rivas et al., 2007) for creating compatible systems.
To develop a taper function, diameter/height data pairs are required along the stem. Most taper functions can be included in the following groups: single and segmented taper models, trigonometric equations, and variable-form taper models. Stem taper functions are usually based on the diameter at breast height (dbh), and predict diameter over bark (Rojo et al., 2005;Crecente-Campo et al., 2009); though they can also predict diameter inside bark (Garber & Maguire, 2003;Calama & Montero, 2006) and bark thickness (Laasasenaho et al., 2005).
Forestry researchers in Spain have been developing taper equations since the 1970s. Recently, a considerable number of taper equations have been developed for particular regions and species: some for hardwoods but most for softwoods (Bravo et al., 2011). Due to their complicated formulations, most taper functions are implemented with specific software in order to estimate total and merchantable volume from inventory data (e.g. SiManFor; www.simanfor.es).
The objective of this study was to compare a stem taper function and a compatible merchantable volume system to ascertain which provides a better description of the stem profile, in order to obtain accurate partial or total stem volume estimates for the main species in the Spanish plateau (central Spain).   Figure 1 shows relative height against relative diameter for the data set used. The range of the data reflects the magnitude of variation in stem form among the sample trees. To detect possible anomalies in the data and increase the efficiency of the process, the systematic approach proposed by Bi (2000) was applied. This approach is flexible because no assumptions about the parametric form of the regression model are needed. In this case, a local quadratic fitting with a smoothing parameter of 0.3 was used for all components. The maximum number of extreme values cor-responded to Spanish juniper (around 20%). At the opposite end of the spectrum, extreme values for black pine and poplar were only 4% and 1%, respectively. In the remaining species, extreme values represented approximately 10% of the total. Most of these data points corresponded to stem deformations, large knots and other physical damage. Since taper functions are not intended for deformed stems, these data points (but not the whole tree) were excluded from further analysis. This approach corresponds to the LOESS procedure, using SAS/STAT software (SAS Institute Inc., 2010a). 4 The Fang system assumes three sections with a variable-form constant factor for each one. The expression of this model is as follows: where:

Functions selected for comparison
Numerous taper functions have been developed and many describe the diameter along the stem quite well. Among them, the segmented function of Fang et al. (2000), which we will refer to as the Fang system, and the Stud model, or Stud variable exponent function developed by Daquitaine et al. (1999), have shown very good results in many studies of several species in Spain (Rojo et al., 2005;Diéguez-Aranda et al., 2006;Barrio-Anta et al., 2007a;Rodriguez et al., 2010). These two outperformed other functions in preliminary analyses and were therefore selected for further analysis.
Variable-exponent taper equations describe the stem shape with a changing exponent from ground to top. This approach is based on the assumption that the stem form varies continuously along the length of a tree. The Stud model is basically an allometric function of the form d=u·(h/H) q , where u is an exponential function that describes the butt region and q is the exponent term describing the tree form. The parameters are dendrometrically and biologically interpretable: θ 1 and θ 2 describe the upper and middle stem, respectively. The other parameters pertain to the function u, which describes the butt region, where θ 3 refers to width, θ 4 refers to length and θ 5 to height. The expression of this model is: We followed the indications of Fang et al. (2000), for avoiding problems in the estimation of their system parameters when h = H; i.e., when d = 0. A small value, lower than the appreciation limit used in the data collection, was reassigned to diameters equal to zero. This approach does not significantly change the parameter estimates (Diéguez-Aranda et al., 2006;Corral-Rivas et al., 2007;Crecente-Campo et al., 2009).
Among the different options to estimate the parameters in the Fang systems, where the taper equation includes a total volume equation, in this study we prioritized the taper function, fitting it first and subsequently performing the predicted volume calculation from the estimation parameters obtained (Menéndez-Miguélez et al., 2014). In this way, we achieved a more accurate comparison between the estimations of the Fang system and the Stud model where the same variable (d) was optimized during the fitting process in both cases.

Model comparison and validation
Estimates of the different fitted models were compared by numerical and graphical analyses. Four goodness-of-fit statistics were used: the adjusted coefficient of determination (R 2 adj ), the mean bias error (BE), the root mean square error (RMSE), and the Bayesian Information Criterion (BIC).
The taper functions were also assessed using box plots for diameter residuals (d) by relative height along the stem (5%, 15%, 25%, and so on up to 95%). The same was done for volume residuals by diameter classes. These graphs, assessed by class for relative height or diameter at breast height, are very important for visualizing areas or tree size classes for which the functions provide especially good or especially poor predictions (Kozak & Smith, 1993).
An n-way cross-validation of each species was carried out, estimating the residual for one tree by excluding that tree every time. The RMSE and the adjusted model efficiency (MEF adj ), equivalent to the R 2 adj of the fitting phase) were calculated from the residuals.

Species differences in the taper equations
To evaluate whether the taper equations vary among the different species, the nonlinear extra sum of squares method was used (Bates & Watts, 1988). This is frequently applied for comparing different geographic regions (Crecente-Campo et al., 2009) or analyzing differences among species (Corral-Rivas et al., 2007;Rodríguez et al., 2010). The nonlinear extra sum of and a i , b i and p i are the parameters to be estimated. Fang et al. (2000) also derived a compatible model for merchantable volume (v, in m 3 ) and total volume (V, in m 3 ) by direct integration of the taper model. Their expressions are: Although Eq. (4) was used to develop the compatible Fang system, any other volume equation can be used as input into the system.

Model fitting
The models were fitted using least squares techniques. However, there are several problems associated with stem taper function analysis that violate the fundamental least squares assumptions of independence of errors: the two most important are multicollinearity and autocorrelation (Kozak, 1997). The condition number (CN) was used to evaluate the presence of multicollinearity among variables in the models analyzed. Since the database contains multiple observations for each tree, it is realistic to expect that the observations within each tree are spatially correlated. The error terms were modeled using a continuous autoregressive error structure (CAR(x)), which makes it possible to apply the model to irregularly spaced, unbalanced data and overcome possible autocorrelation from the longitudinal data sets used for model fitting (Gregoire et al., 1995). Autocorrelation plots are commonly used to check randomness in a data set by computing autocorrelations for data values at varying time lags. The presence of autocorrelation and the order of the CAR(x) were also assessed. Models were fitted using the SAS/ETS MODEL procedure (SAS Institute Inc., 2010b). 6 order continuous autoregressive error structure was required to model the inherent autocorrelation of the hierarchical data. An exception was found in the case of Spanish juniper in the Fang system, where the model did not converge. The autocorrelation may be explained by the effect of stand conditions (e.g., stand density, thinning effects, etc.) on stem form (Calama & Montero, 2006), because stand conditions have a great impact on crown length and stem diameter (Burkhart & Tomé, 2012). Figure 2 shows an example of the autocorrelation plots of the residuals obtained for Scots pine (sp=21).
All the parameters were significant at p<0.05 (Table  2), except in the case of Spanish juniper in the Stud model, where a 11 and a 4 were not significant. The detailed error analysis of diameter predictions showed a similar trend in all the fits evaluated (sixteen combinations of the residuals, two models for each species). Although the poplar analysis rendered poorer graphics than the other species, all cases indicated a random pattern of residuals around zero ( Figure 3). All species obtained good data fits in both models and both phases. The root mean square error obtained in the cross-validation phase ranged from 1.01 to 1.49 times higher and was 1.22 times higher on average than those obtained in the fitting phase (Table 3). More than 99% of the total variance of the diameter was explained for poplar plantations and black pine; more than 96% of total variance was explained for Spanish juniper. As expected, taper functions for stone pine gave worse predictions than equations obtained for other pines, primarily due to its high variability in bark thickness. The greatest RMSE was found in stone pine (RMSE S-TUD = 2.7241; RMSE FANG = 2.6939), while the largest bias occurred in beech (BE STUD = -0.1057; BE FANG = -0.1882). By contrast, smaller RMSE were found in poplar plantations (RMSE STUD = 0.7531; RMSE FANG = 0.7574) and smaller BE in Scots pine (BE STUD = -0.0689; BE FANG = -0.0314). In all species, the Stud model and Fang system showed a similar behavior pattern for the root mean square error against classes of diameter at breast height: accuracy decreased with increasing size class (Table 4). Multicollinearity in both models was moderate, as was inferred from the condition number (Table 3), with values ranging from 9 to 30 for the Stud model and 34 to 110 for the Fang system.
The box plots of diameter residuals against relative height classes (Figure 4) did not show any clear systematic trend that could describe deficient behavior in the models or any clear differences between the Stud model and the Fang system.
Results of the fitting process for full and reduced forms of the model are shown in Table 5. All of the 10 possible paired comparisons in softwoods and the 3 possible paired comparisons in hardwoods produced squares method is based on the likelihood-ratio test for detecting simultaneous homogeneity among parameters and requires the fitting of full and reduced models. The reduced model consists of the same set of parameters for all species and the full model incorporates the different sets of parameters for each species. The full model was obtained by expanding each global parameter to include an associated parameter and a dummy variable to differentiate species. If the F-test (Eq. 9 above) results revealed that there were no differences (p > 0.05) for different species, only a composite model fitted to the combined data was needed. The nonlinear extra sum of squares follows an F-distribution and uses the expression: where SSE R is the error sum of squares of the reduced model, SSE F is the error sum of squares of the full model, and df R and df F are the degrees of freedom of the reduced and full models, respectively.

Results
The models were first fitted without expanding the error terms to account for autocorrelation, thus a strong autocorrelation among all models was observed. A first-

Discussion
Numerous taper equations have been developed for all but two of the species analyzed. Stem taper functions for Spanish juniper and Pyrenean oak have not yet been developed. This paper provides a good starting point for fitting these species, for which only growth (Adame Note: Sp = species as defined in Table 1. Asterisk (*) indicates a not significant parameter estimate at a probability level of 5%. 8 (RMSE around 1.55). The stem form of stone pine has only been characterized in Spain by Calama and Montero (2006), who obtained a 60% lower RMSE (1.71 compared to 2.72 obtained in this study) because they analyzed the diameter inside the bark. There are few taper models for black pine, (Meridieu, 1998), and the results are similar to those observed in Central Spain (RMSE = 0.84). Maritime pine is a widely modeled species (Rojo et al., 2005). Accuracy in the Mediterranean variety (RMSE: 1.479), which is usu-et al., 2008) and biomass (Carvalho & Parresol, 2005) have been modeled so far. In terms of accuracy, the proposed models show that RMSE ranged from 0.75 to 2.72 depending on the species analyzed; a range that is similar to what has been reported in other softwood and hardwood studies. Taper functions for Scots pine have been widely studied in Spain (Diéguez-Aranda et al., 2006;Crecente-Campo et al., 2009) and Europe (Lappi, 1986;Petersson, 1999;Karlsson et al., 2002). All report small errors, similar to those obtained in this work    Stem taper equations in the Spanish Plateau Table 4. Goodness-of-fit for the species and models evaluated in the fitting phase and in the cross-validation phase. Note: Sp = species as defined in Table 1. Asterisk (*) indicates a not significant BE at a probability level of 5%. Note: Sp = species as defined in Table 1. Asterisk (*) indicates a not significant BE at a probability level of 5%. 10 ally measured in natural stands, is only 10% lower than in the Atlantic variety (RMSE: 1.642), which is usually measured in restocked stands. Hardwood species displayed noteworthy goodness-of-fit. Very high precision was obtained for Pyrenean oak (RMSE<1.2, on average 30% lower than other oaks, (Tarp-Johansen et al., 1997;Barrio-Anta et al., 2007a)), which explained more than 98% of the total variance of the dependent variables. There are many published papers on poplar plantations, due to their high commercial value (Roda,  Beach -Sp = 71

Fitting phase Cross-validation phase
Mediterranean maritime pine -Sp = 26 Stud model (Daquitaine et al., 1999) Fang system (Fang et al., 2000) Stud model (Daquitaine et al., 1999) Fang system (Fang et al., 2000) Stud model (Daquitaine et al., 1999) Fang system (Fang et al., 2000)  2001; Barrio-Anta et al., 2007b;Rodriguez et al., 2010). In this work, the poplar was the most successful species, probably due to minimal variability in its management and the fact that all data corresponded to the same clone. The results of the present work are slightly better than those reported in other published papers (about 10% more accurate than Barrio-Anta et al., 2007b). Finally, results for beech were similar to those observed in other European locations (Trincado et al., 1997;Stoltze, 2000).

11
Stem taper equations in the Spanish Plateau use of the models. For sections closer to the ground level, both models provided good estimates. Accurate diameter predictions in these sections are more relevant because the basal log is particularly important from a commercial point of view. According to Crecente-Campo et al., (2009), this behavior is common in stem taper predictions. All these statistics and plots revealed no clear advantage of the Stud model or the Fang system. However, Diéguez-Aranda et al., (2006) consider the Fang system to be advantageous because it is compatible with both a merchantable and a total volume equation. Additionally, a new or pre-existing volume equation can be used as input into the system, making its application more flexible. The non-linear extra sum of squares method indicated that the stem taper differs among the five softwood species and three hardwood species. All of the 10 possible paired comparisons in softwoods and the 3 possible paired comparisons in hardwoods produced significant F-values. In softwoods, the greatest differences (as inferred from the F-values) occurred between black pine and Spanish juniper (F-value = 477.94) while the smallest differences were found between black pine and Mediterranean maritime pine (F-value = 4.63). In hardwoods, the greatest differences appeared between beech and hybrid poplar (F-value = 505.73) while the smallest differences were observed between Pyrenean oak and beech (F-value=41.93). These results are probably due to the strong apical dominance of the hybrid poplar compared to the Pyrenean oak and the beech, and their systematic sylvicul- The results for Scots pine were similar to other findings (Crecente-Campo et al., 2009), where small differences in the goodness-of-fit for both procedures were also found. The Stud model provided lower errors for 6 out of the 8 species, with reductions of 0.004 to 0.111 cm in RMSE (or 0.5% to 8.6%). Both models estimated diameters along the tree reasonably well across the diameter classes that were sampled. It should be noted that the bias tends to increase with the diameter class (Diéguez-Aranda et al., 2006;Corral-Rivas et al., 2007). Cao et al., (1980), found that the variable exponent model represented stem shape quite accurately, especially in the high-volume butt region. The pattern of d residual plots against relative height classes (Figure 4) was similar to corresponding plots for other species (Garber & Maguire, 2003;Crecente-Campo et al., 2009). Both models tended to underestimate the diameters in the lower and upper sections, whereas the midsection diameters were overestimated. For relative heights of 0-10%, both models showed larger standard errors of the estimates than at other height intervals, except in poplar plantations. In the remaining sections, the standard errors of the estimates were more or less constant, depending on the species analyzed. Because stem analysis was usually stopped at a diameter of 7.5 cm, few measurements exist in the top sections and these results should be considered carefully. However, as the latter part of the stem accumulates least volume and is the least valuable, these results do not have a great impact on the overall performance and applied Table 6. F-test of the differences between species obtained with the Stud model (Daquitaine et al., 1999)  ture. All the pines studied were found to have the first inflection point at around 10% of total height, and the second inflection point at around 70% of total height (Table 2), with the exception of stone pine. Similar inflection points have been obtained for other pines in different locations around Spain (e.g., Crecente- Campo et al., 2009). In hardwoods, the first inflection point is lower than in softwoods (around 5% of total height) and the second is higher (around 85%), except in poplar plantations, where the second inflection point is similar to that of pines.