Introduction

Tree total height (H) and diameter at breast height (DBH) are the most commonly used variables in forest management and modelling, and therefore the most widely measured in forest inventories. The basis of growth and yield systems depends on accurate measurements of these variables (^{Misik, Antal, Kárász, & Tóthmérész, 2016}; ^{Sharma, 2016}), as they are used in the simulation of tree growth (^{Mehtätalo, de-Miguel, & Gregoire, 2015}); furthermore, H is needed to estimate volume, biomass, carbon content and economic value (^{Gonzalez-Benecke et al., 2018}; ^{MacPhee et al., 2018}). The H measurement can be costly and time consuming, so it is taken on a subsample in forest inventory plots, while DBH is easily measured (^{Adame, del Rio, & Cañellas, 2008}). In the sampling plots it is common to measure DBH from all trees and some of them are selected to measure H, then the rest are generally estimated with H-DBH equations (^{Santos, Terra, Chaer, & Monte, 2018}; ^{Sharma, 2016}) due to the exhibited correlation between these variables (^{Calama & Montero, 2004}; ^{Zang, Lei, & Zeng, 2016}). These equations are often used to minimize costs associated with forest inventories and to reduce problems related to measurement error (^{MacPhee, Kershaw, Weiskittel, Golding, & Lavigne, 2018}).

The H-DBH relationship has been studied with several approaches and types of datasets. Of the main approaches, local or generalized equations have been used in forestry research. A local equation has the form H = f (DBH, ω), while a generalized equation has the form H = f (DBH, BA, HD, ω) with plot or stand variables such as basal area (BA), dominant height (HD) and density (N); in both cases, ω is the parameter vector. In the last two decades, local and generalized equations have been studied in many species, and the most used approach has been nonlinear mixed effects models (^{Calama & Montero, 2004}; ^{Crecente-Campo, Tomé, Soares, & Diéguez-Aranda, 2010}; ^{Gómez-García, Diéguez-Aranda, Castedo-Dorado, & Crecente-Campo, 2014}; ^{Kearsley et al., 2017}; ^{MacPhee et al., 2018}; ^{Saunders & Wagner, 2008}; ^{Sharma, 2016}). In addition, quantile regression, additive generalized model or Dummy variable approaches have been used (^{Duan, Gao, Wang, & Fu, 2018}; ^{Zang et al., 2016}). Similarly, in the Mexican research, local and general equations have been reported for species and species groups (^{Corral-Rivas, Alvarez-Gonzalez, Crecente-Campo, & Corral-Rivas, 2014}; ^{Crecente-Campo, Corral-Rivas, Vargas-Larreta, & Wehenkel, 2014}; ^{Vargas-Larreta, Castedo-Dorado, Alvarez-Gonzalez, Barrio-Anta, & Cruz-Cobos, 2009}). The datasets are often based on H-DBH measurements from forest inventories, measurement and remeasurement plots, or permanent plots. In most studies, generalized equations perform better than local equations; when a local equation is expanded with stand variables, the application may occasionally be complex, because such variables are not always available.

In the development of H-DBH models, they should capture the full range of population variation, and datasets should take this range with an appropriate sampling design and size (^{MacPhee et al., 2018}). In this study, equations of outside-bark DBH as a function of inside-bark DBH were developed and the H-DBH relationship was extended to the H-DBH, H-A and DBH-A relationships by using age (A). Also a prediction and projection equation systems were fitted for Pinus arizonica Engelmann and Pinus durangensis Martinez with Dummy variable approach and using longitudinal measurements of stem analysis data.

Materials and methods

Area of study and data

Data come from 46 and 66 stem analysis trees of P. arizonica (Pa) and P. durangensis (Pd), respectively, collected in 2006 in northwest Durango, Mexico. The forest region of study was the ejido San Diego de Tezains, geographically located at 24° 48′ 16.98′′ - 25° 13′ 47.25′′ LN and 105° 53′ 9.81′′-106° 12′ 52.58′′ LO (Figure 1). The total area is 61 089 ha, of which 28 636 ha are using for timber production (^{Quiñonez-Barraza, Zhao, De los Santos Posadas, & Corral-Rivas, 2018}). The main climate is cold temperate, the average annual temperature varies from 5 to 12 °C and the average annual rainfall is 840 mm (^{García, 2004}; ^{Quiñonez-Barraza et al., 2019}). About 85 % of the stands are uneven-age and mixed and they are represented by the genus Pinus, Quercus, Juniperus, Cupressus, Pseudotsuga, Arbutus and Alnus. The stands have been managed according to the continuous cover forestry system with individual selection and natural regeneration, or the rotation forest management system characterized by three or four thinnings and a shelterwood cut (^{Pukkala & Gadow, 2011}; ^{Quiñonez-Barraza et al., 2018}). The cutting cycle is 15 years (^{Quiñonez-Barraza et al., 2018}).

Figure 1 Study area and location of sampled trees from Pinus arizonica (Pa) and Pinus durangensis (Pd) in northwest Durango, Mexico. 

Sampled trees were randomly selected in mixed-species stands to cover canopy types, densities, heights, DBHs, stand structures, ages, and productivity (Figure 1). All trees were cut as close to the ground as possible; each fallen tree was sectioned and measured at 1.3 m, and subsequently every 2 m to the tree top. The DBH, commercial H and number of growth rings were measured in each cross-section. Growth in diameter was examined with the width of the annual growth rings in the DBH cross-section. In the stem analysis trees, 601 and 760 repeated measurements of H-DBH, H-A, and DBH-A were considered for Pa and Pd, respectively. A summary of the individual tree dataset and the measurements of the stem analysis variables for each species is shown in Table 1. Figure 2 illustrates the relationship of outside-bark DBH (DBHob) as a function of inside-bark DBH (DBHib) corresponding to the cross-section at 1.3 m and Figure 3 shows the trends for the H-DBHob, H-A and DBHob-A relationships. The H corresponding to an age for a specific cross section was estimated with the ^{Carmean method (Carmean, 1972}; ^{Dyer & Bailey, 1987}). With the data from stem analysis, four relationships were tested and modeled for each species: (1) DBHob-DBHib, (2) H-DBHob (H-DBH), (3) H-A in the cross section of the DBH, and (4) DBHob-A (DBH-A) in the cross section of the DBH.

Table 1 Variables of sampled trees (1) and stem analysis (2) of Pinus arizonica and Pinus durangensis in mixed-species stands of Durango, Mexico. 

Data	Variable	Pinus arizonica					Pinus durangensis
		n	Min	Max	Mean	SD	n	Min	Max	Mean	S
1	H (m)	46	8.3	30.5	18.1	4.7	66	8.7	26.5	16.1	4.6
	DBH (cm)		14.0	51.0	27.2	10.3		13.0	45.0	24.6	8.9
	A (years)		39.0	190.0	76.7	29.8		35.0	172.0	75.2	30.7
2	H (m)	601	1.3	30.5	11.4	6.5	760	1.3	26.5	9.9	5.7
	DBHib (cm)		0.0	44.0	14.5	9.6		0.0	43.6	12.4	8.4
	DBHob (cm)		0.0	48.7	16.3	10.7		0.0	46.5	14.2	9.1
	A (years)		3.0	190.0	43.1	27.8		2.0	172.0	42.5	27.6
	A-DBH (years)		0.0	181.0	34.9	27.9		0.0	167.0	34.7	27.6

n = number of observations, Min = minimum, Max = maximum, Mean = media, SD = standard deviation, H = height, DBH = diameter at breast height, DBHib = outside-bark DBH, DBHob = inside-bark DBH, A = number of rings on each tree section, A-DBH = age on DBH cross section (i.e., number of rings on 1.3 m cross section).

Figure 2 Outside-bark DBH (DBHob) and inside-bark DBH (DBHib) in cross section at 1.3 m for Pinus arizonica (Pa) and Pinus durangensis (Pd). 

Figure 3 Distribution of H-DBH, H-A and DBH-A relationships (H = height, DBH = diameter at breast height and A = age) for Pinus arizonica (Pa) and Pinus durangensis (Pd), from repetitive measurements of stem analysis. 

Results and discussion

The preliminarily adjusted equations modeled the DBHob as a function of the DBHib for each species and the potential was better for both species. The results of this equation allowed the modeling of H-DBH and DBH-A relationships with repetitive measurements of DBHib measured in the DBH cross section. Table 2 shows the estimated parameters and the statistics of fitting of the potential, linear and exponential equations. The values of R ² with the potential equation were 0.9804 and 0.9664 for Pa and Pd, respectively. The linear equation was similar to the potential, but the intercept parameter was significantly equal to zero (α = 0.05) for Pa, while the exponential equation showed the poorest fitting.

The application of the H-DBH or DBH-A equations should use measurements of the DBHob, because it is usually measured in forest inventories and used to predict H, which reduces the problems associated with measurement error (^{MacPhee et al., 2018}; ^{Misik et al., 2016}). If DBH measurements are derived from stem analysis data, each DBHib measurement should be transformed to DBHob. Then, the H-DBH relationship can be established to predict H, which is one of the most important characteristics in forest structure (^{Mehtätalo et al., 2015}; ^{Wang et al., 2017}). The H and DBH variables can be predicted and projected as functions of the A of the DBH cross section and used in the characterization of stand productivity (^{Duan et al., 2018}; ^{Fu et al., 2018}). The estimated parameters in the three relationships studied are listed in Table 3. All parameters, both common and specific, were significantly different from zero at a significance level of 5 %.

Fitting statistics such as R², RMSE, AIC and E are shown in Table 4. The variance-covariance matrices for the prediction and projection equations, for each relationship and species, and the descriptive statistics for the specific adjusted parameters are found in Appendix B1. For the H-DBH relationship, the prediction model had R² values of 0.8210 and 0.8521 and RMSE values of 2.7238 m and 2.1717 m, for Pa and Pd, respectively. The projection equation showed better results than the prediction equation. The R² values were 0.9611 and 0.9595, while the RMSE values were 1.2212 m and 1.0865 m, for Pa and Pd, respectively. In the H-A equations, the results were similar among species; the R² values were 0.7710 and 0.7086 with the prediction equation and 0.9505 and 0.9435 with the projection equation for Pa and Pd, respectively. In the DBH-A relationship, the prediction equation showed R2 values of 0.7338 and 0.7200, while the projection equation had values of 0.9669 and 0.9560 for Pa and Pd, respectively.

The DVA ensured common and tree-specific parameters in the projection equations. This procedure is similar to the non-linear mixed effect model (NLMEM) approach for the H-A relationship. The DVA also considers a specific or random parameter for each tree (^{Diéguez-Aranda et al., 2006}; ^{Wang et al., 2008}). The H-DBH, H-A and DBH-A relationships were fitted as equation systems of prediction and projection with the compatible growth rate parameter (α₁). The projection equations presented specific parameters for each tree (summary information is found in Appendix B2). The GADA equation used (^{Quiñonez-Barraza et al., 2015}) allowed the projection of exogenous variables with common and tree-specific parameters for the established relationships.

The values of SEE and E and the relative ones (RSEE and RE) for the prediction and projection equations by relative H classes, in the relations H-DBH and H-A, are showed in Table 5, and by relative DBH classes for DBH-A are shown in Table 6. Pa had higher values of SEE and E than Pd for all data and by relative classes of H or DBH. The largest values of SEE and E were observed in the first relative classes of H and DBH for both species.

The percentages of gain in R ² values between the projection and prediction equations were 14.6 % and 11.2 % and the percentages of reduction in RMSE of 123.0 % and 99.9 % for Pa and Pd, respectively. In contrast, ^{Vargas-Larreta et al. (2009)} reported gains in R² of 1.2 % and 2.4 % and reductions in RMSE of 3.72 % and 6.09 % for Pinus cooperi C. E. Blanco between the NLMEM and OLS procedures, respectively, with the generalized equation based on the Chapman-Richards equation and the inclusion of the quadratic mean diameter. ^{Corral-Rivas et al. (2014)} determined R² gains of 14 % and 2 % and RMSE reductions of 33 % and 5.2 % with NLMEM and OLS between generalized and local models based on the Chapman-Richards equation, respectively, for the genus Pinus. Similarly, ^{Crecente-Campo et al. (2014)} indicated an RMSE reduction of 14.9 % between general and local equations fitted by OLS for all species combined, as well as a 7.5 % reduction between a general and local model fitted by NLMEM. Therefore, the local equations developed in this study perform better than the equations reported for the mixed and uneven-aged forests of Durango, Mexico. The H-DBH relationship fitted by NLMEM usually considers the variability between stand-stand, plot-plot, or ecoregion-ecoregion (^{Calama & Montero, 2004}; ^{Crecente-Campo et al., 2014}; ^{Mehtätalo et al., 2015}; ^{Paulo, Tomé, & Tomé, 2011}; ^{Rijal, Weiskittel, & Kershaw, 2012}); however, in this study, the variability was assumed between tree and tree.

The residual and predicted values for the prediction and projection equations in the H-DBH, H-A, and DBH-A relationships are shown in Figures 4 and 5. In all cases, the graphs showed the residual trajectories for the growth series of each tree.

Figure 4 Residual versus predicted values for the H-DBH (Height-Diameter at breast height) and H-A (Height-Age) relationship of the prediction (left) and projection (right) equations for Pinus arizonica (Pa) and Pinus durangensis (Pd). 

Figure 5 Residual versus predicted values for the H-A (height-age) and DBH-A (Diameter at breast height-Age) relationships of the prediction (left) and projection (right) equations for Pinus arizonica (Pa) and Pinus durangensis (Pd). 

The three relations of prediction and projection performed well for Pa and Pd and the inverted equations (Appendix B3) formed a global system of equations. Those inverted equations hold the four criteria cited by ^{Bi, Fox, Li, Lei, and Pang (2012)}, for the H-DBH relationship and the extended H-A and DBH-A relationships. These criteria are: (1) the H-DBH, H-A, and DBH-A functions are reversible; (2) the inverted functions are continuous and monotonically increasing over a specific range of H and DBH; (3) the DBH is equal to zero when H is equal to DBH in the inverted function, and A in the cross section of the DBH is equal to zero when H or DBH are equal to A in the inverted function; (4) the inverted functions have an inflection point that is consistent with biological expectations.

The DVA approach is similar to NLMEM; common parameters can be compared with fixed parameters and specific parameters with random parameters (^{Duan et al., 2018}; ^{Wang et al., 2008}). The H-A, DBH-A and H-DBH relationships may be interpreted as a site productivity index for even-age or uneven-age forests (^{Duan et al., 2018}; ^{Fu et al., 2018}; ^{Wang et al., 2017}). However, the usefulness of this index as a measure of site quality for uneven-age and mixed-species stands should be questioned and perhaps rejected, since age classes are diverse within and between mixed stands (^{Duan et al., 2018}; ^{Wang, 1998}). Consequently, the objectives of this study were to develop and illustrate prediction and projection equation systems for Pa and Pd H-DBH-A relationships growing in natural and mixed stands. Therefore, in the H-DBH relationship, the theoretical assumption described by ^{Fang and Bailey (1998)} was guaranteed in the prediction and projection equations, in which the DBH of a tree should be equal to zero, when the H is equal to 1.3 m. This theoretical assumption was followed for the H-A relationship; the A of a tree in the DBH cross section was zero when the H was equal to 1.3 m. On the other hand, for the DBH-A relationship; the A of the tree was zero when the DBH was zero as well. In all cases, the A of the trees was taken as the number of rings in the cross section corresponding to the DBH. Thus, individual tree local equations can be applied to predict the variables of interest in a wide range of DBHs and Hs measurements, and in a variety of ecological conditions and silvicultural treatments (^{Gonzalez-Benecke et al., 2018}). Additionally, a generalized equation may improve these applications in diverse stand-specific characteristics (^{Gonzalez-Benecke et al., 2014}).

The least-squares procedure proposed by ^{Wang et al. (2008)} could be used to predict the specific parameter for an H or DBH from a new growth series in the projection equation, similar to the search for random parameters in the NLMEM procedure (^{Crecente-Campo et al., 2014}). The stem analysis data guaranteed criterion 3 noted by ^{Bi et al. (2012)} for the Al parameter in both the prediction and projection equations and the corresponding inverted equations. In cases where the dataset does not cover this property; the lower asymptote should be guaranteed with the addition of the minimum measurement of DBH or A on the right side of the equation.

Conclusions

In this study, the H-DBH-A (Height-Diameter at breast height-Age) relationships were introduced as systems of independent prediction and projection equations for the species Pinus arizonica and Pinus durangensis. The developed equations allowed the prediction and projection of Hs and DBHs as a function of endogenous variables. With the inverted equations, a global system of equations was integrated when the measured variables were the DBHob (inside-bark diameter at breast height) and A (Age). The apparently unrelated regression method allowed the adjustment of the H-DBH, H-A and DBH-A relationships as prediction and projection equations with the focus on Dummy Variables Approach and common and specific parameters. In such a context, equation systems can be used as a powerful biometric tool in forest management or as input variables in growth and yield models; furthermore, they can maximize the cost effectiveness of timber volumes in the mixed-species forests of Durango, Mexico.