Local and generalized height-diameter models with random parameters for mixed, uneven-aged forests in Northwestern Durango, Mexico

Abstract

Background

We used mixed models with random components to develop height-diameter (h-d) functions for mixed, uneven-aged stands in northwestern Durango (Mexico), considering the breast height diameter (d) and stand variables as predictors.

Methods

The data were obtained from 44 permanent plots used to monitor stand growth under forest management in the study area.

Results

The generalized Bertalanffy-Richards model performed better than the other generalized models in predicting the total height of the species under study. For the genera Pinus and Quercus, the models were successfully calibrated by measuring the height of a subsample of three randomly selected trees close to the mean d, whereas for species of the genera Cupressus, Arbutus and Alnus, three trees were also selected, but they are specifically the maximum, minimum and mean d trees.

Conclusions

The presented equations represent a new tool for the evaluation and management of natural forest in the region.

Background

Most forests in Durango State (Mexico) are comprised of a mixture of species of the genera Pinus and Quercus with an irregular distribution of trees of all size classes. However, species of the genera Arbutus and Juniperus are also found in most of these forests (Wehenkel et al. 2011). These forests, which cover an area of 5.4 million ha, are considered as the primary forest reserve at a national level, and they provide almost a quarter of the national forest production in Mexico (SRNyMA 2006). The forests also play an important role in providing environmental services, such as protection against soil erosion, biodiversity conservation, carbon capture and protection of water reserves; they also provide recreational areas and represent an important source of income for their owners and local inhabitants.

Forest management requires prediction tools that provide detailed information about the development of mixed, uneven-aged stands. Growth and production models are the most commonly used tools for this purpose. When the breast height diameter (d) and total height (h) are known, application of these models is relatively easy (Sharma and Parton 2007). Measuring diameter is simple, accurate and inexpensive, whereas measuring height is relatively more complex, time-consuming and expensive. Therefore, height-diameter (h-d) functions are often utilized, so that the height of an individual tree can be predicted only from the diameter. These relationships are also very useful for estimating individual volume, site index and for describing growth and production in forest stands over time when the height is not measured (Curtis 1967).

Most h-d functions have been developed for forest plantations (e.g. Soares and Tomé 2002; López Sanchez et al. 2003). However, the relationship between the diameter and height of a tree varies between stands (Calama and Montero 2004) because it depends on stand characteristics such as density and site index (Sharma and Zhang 2004). Moreover, the h-d relationship also varies over time within the same stand (Curtis 1967). Such considerations indicate that stand variables should be used to construct generalized functions that represent all possible conditions in forest stands (Temesgen and Gadow 2004). This is particularly important in mixed, uneven-aged stands in which different species, ages, structures and levels of competition coexist (Vargas-Larreta et al. 2009).

The hierarchical structure of h-d data (i.e. trees grouped in plots and plots grouped in stands) results in a lack of independence between measurements because the observations in each sampling unit will be correlated (Gregoire 1987). Mixed models have been successfully used to address this type of problem (e.g. Lappi 1997; Calama and Montero 2004; Castedo Dorado et al. 2006). This approach simultaneously estimates fixed parameters (parameters that are common to the entire population) and random parameters (parameters that are specific to each plot) within the same model and enables the variability between plots of the same population to be modelled.

The objectives of this study were as follows: i) to compare different local h-d equations for the mixed, uneven-aged forests in north-western Durango; ii) to develop new generalized h-d equations for different groups of species based on the best local model previously fitted iii) to use the local and the generalized equations to study the capacity of mixed models to explain the variability in the h-d relationship; and iv) to determine the most suitable size and type of sample for calibrating the functions fitted with mixed models.

Methods

Study area

The study was carried out in the Ejido San Diego de Tezains, Municipality of Santiago Papasquiaro, Durango State, Mexico (between 105° 53′ 36″ and 106° 12′ 40″ W and 24° 48′ 16″ and 25° 13′ 32″ N). The predominant vegetation in the area is mixed, uneven-aged forests of Pinus and Quercus. The altitude above sea level of the study area varies between 1,400 and 3,000 m. The prevailing climate is temperate: the annual precipitation ranges between 800 and 1,100 mm and the mean annual temperature varies between 8°C in the highest elevations and 24°C in the lowest elevations (García 1981).

Data

The data were obtained from 44 permanent plots used to monitor the growth and production of the forests in the Ejido San Diego de Tezains. These plots, which were established in 2008, were selected with the aim of representing all types of vegetation, site qualities and diameter distributions in managed stands. The plots, of size 50x50 m, are distributed under a systematic grid sampling approach that varies between 3 and 5 km, and will be remeasured at 5 year intervals. We recorded the following main variables: number of trees, species code, breast height diameter at 1.3 m (d, cm), total tree height (h, m), azimuth (°) and radius (m) from the centre of the plot (point where the diagonals cross) towards all trees of breast height diameter ≥5 cm.

The database included 25 species, which were classified on the basis of their growth patterns into the following 13 groups for posterior analysis: 1 (Pinus arizonica), 2 (P. ayacahuite), 3 (P. durangensis), 4 (P. herrerae), 5 (P. lumholtzii), 6 (P. teocote), 7 (P. douglasina), 8 (Quercus sideroxyla), 9 (other species of Quercus: Q. arizonica, Q. mcvaughii, Q. durifolia, Q. crassifolia, Q. jonesii, Q. rugosa and Q. laeta), 10 Pinus species (all species of the genus Pinus [codes 1 to 7]), 11 Quercus species (all species of the genus Quercus [codes 8 and 9]), 12 other conifer species (Juniperus deppeanna, J. durangensis and Cupresus lusitanica) and 13 other broadleaf species (Arbutus arizonica, A. bicolor, A. madrensis, A. tesselata, A. xalapensis and Alnus firmifolia).

We examined the distribution of the pairs of h-d data for each species or group graphically to identify any possible anomalies. As extreme data points were observed, a systematic approach, similar to the one proposed by Bi (2000) for detecting abnormal data points, was applied to increase the efficiency of the process. A local quadratic equation with a smoothing parameter of 0.25 (selected after iterative fitting and visual examination of the smoothed curves for different smoothing parameters overlaid on the data), was fitted for each of the species or group. In this approach, the number of extreme values accounted for about 1% for all species together, which were excluded from the database used for fitting the equations. The main descriptive statistics for the breast height diameter and the total height of the main groups that included the species under study are shown in Table 1.

Table 1 Summary statistics of the database used in fitting the h-d equations

The following stand variables were calculated from the trees registered in each plot: number of trees per hectare (N, trees ha^-1), stand basal area (G, m² ha^-1), mean square diameter (d _g , cm), dominant height (estimated as the mean height of the 100 largest diameter trees per hectare, independently of the species [H₀, m]), dominant diameter (estimated as the mean diameter of the 100 largest diameter trees per hectare, independently of the species [D₀, cm]) and Hart’s index (%) estimated as follows: $\mathit{HI}=\raisebox{1ex}{$10000$}\!\left/ \!\raisebox{-1ex}{$\sqrt{\mathit{N}}*{\mathit{H}}_0$}\right.$.

Comparison of equations

We selected a total of 27 local equations (Huang et al. 2000) for data fitting. We also studied the relationship between the stand variables and the parameters of the local equations that best described the h-d relationship, with the aim of improving the accuracy of the equation and developing new generalized functions.

For preliminary selection, we used ordinary non-linear least squares (ONLS) to fit each of the local equations to the data from the 13 established groups, with the MODEL procedure in SAS/ETS^® statistical software package (SAS Institute Inc 2008). We evaluated the goodness of fit of the models by graphical analysis and by considering the following statistics, calculated from the residuals: root mean square error (RMSE), the coefficient of determination (R²), bias, and Bayesian information criterion (BIC; Schwarz 1978). We used the following formulae to calculate these statistics:

$$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle {\sum }_{\mathit{i}=1}^{\mathit{i}=\mathit{n}}{\left({\mathit{h}}_{\mathit{i}}-{\widehat{\mathit{h}}}_{\mathit{i}}\right)}^2}}{\mathit{n}-\mathit{k}}}$$

(1)

$${\mathit{R}}^2=\left[1-\frac{{\displaystyle {\sum }_{\mathit{i}=1}^{\mathit{i}=\mathit{n}}{\left({\mathit{h}}_{\mathit{i}}-{\widehat{\mathit{h}}}_{\mathit{i}}\right)}^2}}{{\displaystyle {\sum }_{\mathit{i}=1}^{\mathit{i}=\mathit{n}}{\left({\mathit{h}}_{\mathit{i}}-\overline{\mathit{h}}\right)}^2}}\right]$$

(2)

$$\mathit{Bias}=\frac{{\displaystyle {\sum }_{\mathit{i}=1}^{\mathit{i}=\mathit{n}}\left({\mathit{h}}_{\mathit{i}}-{\widehat{\mathit{h}}}_{\mathit{i}}\right)}}{\mathit{n}}$$

(3)

$$\mathrm{BIC}=\mathit{n}ln\left({\displaystyle {\sum }_{\mathit{i}=1}^{\mathit{i}=\mathit{n}}{\left({\mathit{h}}_{\mathit{i}}-{\widehat{\mathit{h}}}_{\mathit{i}}\right)}^2}/\left(\mathit{n}-\mathit{k}\right)\right)+\mathit{k}\cdot ln\left(\mathit{n}\right)$$

(4)

where h _i , ${\widehat{\mathit{h}}}_{\mathit{i}}$ and $\overline{\mathit{h}}$ are the observed and estimated heights and the mean of the observed heights, respectively; n is the number of observations used in the fitting; k is the number of parameters in the equation, and ln is the natural logarithm.

As each local equation has different strengths and weaknesses, which may lead to different goodness-of-fit results for each group of species, we used a Qualification Index (QI _t ) to evaluate the goodness of fit by considering the values of R² (with high values representing good fits), Bias (with low absolute values representing good fits) and RMSE and BIC (with low values representing good fits). For this index, a value of 1 is assigned to the equation that was best for each group of species and a value of 0 to the others. The qualifications obtained for each equation and statistics were then summed as follows: $\mathit{Q}{\mathit{I}}_{\mathit{total}}={\displaystyle \sum _{\mathit{i}}{\displaystyle \sum _{\mathit{j}}\mathit{Q}{\mathit{I}}_{\mathit{ij}}}}$; where QI_ij is the qualification for the j-th goodness of fit criterion in the i-th group of species.

For the local equation for which the QI _total was highest for the defined groups, we used graphical analysis and the CORR procedure in SAS (SAS Institute Inc 2008) to analyse the relation between each of the parameters and the main stand variables, with the aim of testing different forms of generalized equations.

Effect of mixed models

The h-d observations made in plots and stands may be highly correlated, thus violating the principle of independence of error terms (Calama and Montero 2004). One procedure used to deal with correlated observations is to fit mixed models, in which the variability between the sampling units can be explained by including random parameters, which are estimated at the same time as the fixed parameters (Lappi 1997; Calama and Montero 2004).

Basically, the parameter vector of a non-linear mixed model can be defined as follows (Pinheiro and Bates 1998):

$${\mathbf{\Phi }}_{\mathit{j}}={\mathbf{A}}_{\mathit{j}}\mathbf{\lambda }+{\mathbf{B}}_{\mathit{j}}{\mathbf{b}}_{\mathit{j}}$$

(5)

where Φ _j is the parameter vector r × 1 (where r is the total number of parameters in the model) specified for the j-th plot, λ is the vector p × 1 of the common fixed parameters for the whole population (p is the number of fixed parameters in the model), b _j is the vector q × 1 of the random parameters associated with the j- th plot (q is the number of random parameters in the model), A _j and B _j are matrices of size r × p and r × q for specific and random effects for the j-th plot, respectively.

The basic theory of non-linear mixed models says that the residual vector (${\widehat{\mathbf{e}}}_{\mathit{ij}}$) and the random effects vector (b _j ) are often assumed to be uncorrelated and normally distributed with mean zero and variance-covariance matrices R _j and D, respectively. The residual vector represents within subject (e.g., plot) variability and the random effects vector represents between subject variability (Littell et al. 1996).

We constructed the non-linear mixed effects model by selecting the local and generalized equations that yielded the best fits for the species groups defined using the NLMIXED procedure in SAS/ETS (SAS Institute Inc 2008). We tested different combinations of fixed and random parameters and compared the fitting statistics (RMSE, R², Bias and BIC), to determine which parameter(s) should be considered mixed.

Calibration

The inclusion of random parameters in h-d equations leads to two possible situations as regards prediction of the height of trees within a stand (Vonesh and Chinchilli 1997): i) a population mean response (PMR) when only diameters are measured (and the stand variables are included in the model in the case of generalized models) and the vector of random parameters is assumed to have an expected value of E(b _j ) = 0; and ii) a calibrated response, when the height of a subsample of m _j trees is measured along with diameter measurement in each new plot j (and the stand variables in the case of generalized models) and is subsequently used to calculate the specific random parameters of the new sampling units (Calibrated Response; CR), i.e. vector b _j , expressed as follows (Vonesh and Chinchilli 1997):

$${\widehat{\mathbf{b}}}_{\mathit{j}}\approx \widehat{\mathbf{D}}{\widehat{\mathbf{Z}}}_{\mathit{j}}^{\mathrm{T}}{\left({\widehat{\mathbf{R}}}_{\mathit{j}}+{\widehat{\mathbf{Z}}}_{\mathit{j}}\widehat{\mathbf{D}}{\widehat{\mathbf{Z}}}_{\mathit{j}}^{\mathrm{T}}\right)}^{-1}{\widehat{\mathbf{e}}}_{\mathit{ij}}$$

(6)

where $\widehat{\mathbf{D}}$ is the matrix q × q of variances-covariances associated with the random parameters (q = number of random parameters included in the model), which is common to all plots and is estimated in the general model fitting procedure; ${\widehat{\mathbf{R}}}_{\mathit{j}}$ is the m _j × m _j estimated matrix of variances-covariances of the error term; ${\widehat{\mathbf{e}}}_{\mathit{ij}}$ is the residuals vector m × 1, the components of which are obtained as the difference between the observed height of each tree and the value predicted using the model with fixed parameters only; and ${\widehat{\mathbf{Z}}}_{\mathit{j}}$ is the matrix m × q of the partial derivatives of the random parameters evaluated in ${\widehat{\mathbf{b}}}_{\mathit{j}}=0$.

Two sampling options were considered for selecting the subsample of trees to measure within each for calibration of the local and generalized equations:

• (i) CR1: Measuring the total height of between 1 and 5 randomly selected trees within each plot that are close (± 10% ) to the mean breast height diameter.

• (ii)  CR2: Measuring the total height of the tree of mean breast height diameter, or measuring the height of two trees – the mean and minimum breast height diameters, or measuring the height of three trees – the mean, minimum and maximum breast height diameters within each plot.

We evaluated these two alternatives in terms of the previously defined goodness-of-fit statistics (RMSE, R² and Bias), which we compared with the statistics obtained for the equations fitted by the ONLS and NLMIXED procedures.

Results and discussion

Local equations

In the comparison of the goodness-of-fit statistics for the local h-d equations fitted to the data for the 13 previously defined groups of species, the Bertalanffy-Richards equation (Bertalanffy 1949; Richards 1959) consistently yielded the highest R² and lowest RMSE values; however, the equations proposed by Stage (1975) and Meyer (1940) yielded the lowest values of BIC, which gives preference to models containing few parameters over those containing several parameters (Table 2).

Table 2 Qualification index for the 8 best local equations for the 13 groups of species

The Bertalanffy-Richards equation yielded the highest R² values for 6 of the 13 species groupings and the lowest RMSE values for 4 of the groups. Finally, comparison of the BIC values indicated that this was the preferred equation only for the Quercus spp. grouping.

In selecting the best local equation, we also examined graphs of the residuals, the significance of the parameters and the mean bias for each equation. In this respect, the Bertalanffy-Richards equation was the preferred model. The final structure of the local Bertalanffy-Richards equation used was as follows:

$$\mathit{h}=1.3+{\mathit{b}}_0{\left(1-exp\left(-{\mathit{b}}_1\cdot \mathit{d}\right)\right)}^{{\mathit{b}}_2}$$

(7)

where b ₀ -b ₂ are equation parameters and the rest of variables as defined in the data section.

The value of the goodness-of-fit statistics for Eq. (7) and the species groups are shown in Figure 1.

Figure 1

Values of R ² (solid line), RMSE (dashed line), and mean bias (bottom) obtained by fitting Eq. ( 7 ) to the h-d relationship for the 13 groups of species considered.

Considering that the fitting statistics for the different broad groupings (Pinus species, Quercus species, other conifer species and other broadleaf species) are similar to those obtained for each individual species and that some parameters were not significant in the individual fits for some of the species, we decided to use 4 different local equations, one for each broad grouping studied.

The parameters in Eq. (7) and the goodness-of fit statistics, obtained by the ONLS method for the 4 broad groupings, are shown in Table 3.

Table 3 Estimated parameters and fitting statistics obtained for the local model with and without mixed effects for the groups of species considered

Generalized equations

On relating the parameters in Eq. (7) to the stand variables, we found that in the 4 groups that included all species, parameter b ₀ , representative of the asymptote of Eq. (7), was positively correlated (almost 52%) with H₀and D₀, whereas parameter b ₁ , representative of “scale” was only positively correlated by more than 50% with Hart’s index (HI) in the other conifer species and other broadleaf species and, finally, parameter b ₂ , representative of “shape” was also positively correlated (~ 51%) with d _g and N for the four groups that include all species.

To develop a generalized equation from Eq. (7), we tested various combinations of the stand variables to improve their efficacy in the fit, taking into account the previously mentioned correlations; this resulted in Eqs (8) and (9). Equation 8 yielded the best fitting statistics for the Pinus and Quercus groups, whereas Eq. (9) yielded the best fitting statistics for the species included in other conifer species and other broadleaf species.

$$\mathit{h}=1.3+{\mathit{b}}_0\cdot {{\mathit{H}}_0}^{{\mathit{b}}_1}{{\left(1-exp\left(-{\mathit{b}}_2\cdot \mathit{d}\right)\right)}^{{{\mathit{b}}_3}^{\left(\raisebox{1ex}{$\mathit{N}$}\!\left/ \!\raisebox{-1ex}{${\mathit{d}}_{\mathit{g}}$}\right.\right)}}}^{{\mathit{b}}_4}$$

(8)

$$\mathit{h}=1.3+{{\mathit{H}}_0}^{{\mathit{b}}_0}{{\left(1-exp\left(-{\mathit{b}}_1\cdot \mathit{HI}\cdot \mathit{d}\right)\right)}^{{}^{\left(\raisebox{1ex}{$\mathit{N}$}\!\left/ \!\raisebox{-1ex}{${\mathit{d}}_{\mathit{g}}$}\right.\right)}}}^{{\mathit{b}}_2}$$

(9)

where b ₀ -b ₂ are equation parameters and the rest of variables as defined in the data section.

On comparing the goodness-of-fit statistics for Eqs (8) and (9) with those obtained when fitting the 30 generalized equations used in other studies (e.g., López Sanchez et al. 2003; Sharma and Zhang 2004; Sharma and Parton 2007), we found that the value of RMSE for Eq. (8) was slightly lower than those obtained by fitting the above-mentioned generalized equations, whereas the value of the RMSE for Eq. (9) was only lower than some of these. However, the latter equations included parameters that were not significant at the 0.05 level. Another advantage of Eq. (9) is that the value of the BIC was lower than that obtained for the 30 generalized equations compared. Analysis of the residuals also revealed that there were no anomalies associated with Eqs (8) and (9) that would indicate non-compliance with the underlying hypotheses of normality, homogeneity of variance or independence of errors. We therefore decided to use Eqs (8) and (9) as generalized models for the four groups that included all species considered. The values of the parameters of Eqs (8) and (9) obtained for fitting each group of species are shown in Table 4. The signs and values of all parameters are consistent with their biological interpretation and visual examination of the graphs of the h-d relationship indicates that its performance was consistent with the theory of growth.

Table 4 Estimated parameters and fitting statistics obtained for the generalized model with and without mixed effects for the groups of species considered

The generalized h–d equations selected in this study included dominant stand height. This represents an advantage over equations that include the mean height because less effort is required in conventional inventories to estimate the dominant height than the mean height of the stand (López Sanchez et al. 2003). These functions also include the density of the stand in terms of number of trees per unit of area and mean square diameter. Stand density is the most obvious factor affecting the h-d relationship in a stand (Zeide and VanderSchaaf 2002); in other words, trees of the same diameter are generally taller in denser stands.

Various stand variables have been proposed as predictors of the h-d relationship: stand age (Curtis 1967; Soares and Tomé 2002; López Sanchez et al. 2003); crown competition index (Temesgen et al. 2007); geographic variables (Schmidt et al. 2011); and wind speed (Meng et al. 2008). Although the inclusion of other variables may improve the predictive capacity of the selected functions, this requires great sampling effort and limits the practical application of the functions and therefore we did not take such variables into account.

Effect of mixed models

Parameters b ₀ , b ₁ and b ₂ in Eqs (7) and (9) determine the asymptote, the scale and the shape of the h-d curves, respectively, whereas in Eq. (8), the parameters b ₀ and b ₁ define the asymptote, b ₂ is the scale parameter and b ₃ and b ₄ define the shape of the curve. In fitting Eq. (7) to the data from each plot, we found that the parameter that affected the asymptote was the most variable, followed by the scale parameter. Therefore, in a first step, we fitted Eqs (7), (8) and (9) to the h-d data by considering the parameters that define the asymptote and scale as mixed, in other words, with a random parameter added. Pinheiro and Bates (1998) obtained similar results and found that the best results were obtained when the asymptote and scale of Eq. (7) were considered as random parameters. In most cases, the mixed model did not converge, so that we tested the inclusion of only mixed parameters associated with those parameters that define the asymptote of the h-d curve until reaching convergence. Similar results have been reported by Sharma and Parton (2007) and Vargas-Larreta et al. (2009). Some of the parameters were scaled so that all were of the same order of magnitude and to prevent instability in the fitting function (Calama and Montero 2004). The expressions of the mixed models finally obtained are Eqs (10) to (15):

$${\widehat{\mathit{h}}}_{\mathit{ij}}=1.3+\left({\mathit{b}}_0+{\mathit{u}}_{\mathit{j}}\right){\left(1-exp\left(-\left({\mathit{b}}_1+{\mathit{v}}_{\mathit{j}}\right){\mathit{d}}_{\mathit{ij}}\right)\right)}^{{\mathit{b}}_2}+{\mathit{e}}_{\mathit{ij}}$$

(10)

$${\widehat{\mathit{h}}}_{\mathit{ij}}=1.3+\left({\mathit{b}}_0+{\mathit{u}}_{\mathit{j}}\right){\left(1-exp\left(-{\mathit{b}}_1\cdot {\mathit{d}}_{\mathit{ij}}\right)\right)}^{{\mathit{b}}_2}+{\mathit{e}}_{\mathit{ij}}$$

(11)

$${\widehat{\mathit{h}}}_{\mathit{ij}}=1.3+\left({\mathit{b}}_0+{\mathit{u}}_{\mathit{j}}\right)\cdot \mathit{H}{\mathit{o}}^{{\mathit{b}}_1}{{\left(1-exp\left(-{\mathit{b}}_2\cdot {\mathit{d}}_{\mathit{ij}}\right)\right)}^{{{\mathit{b}}_3}^{\left(\raisebox{1ex}{$\mathit{N}$}\!\left/ \!\raisebox{-1ex}{${}_{{\mathit{d}}_{\mathit{g}}}$}\right.\right)}}}^{{\mathit{b}}_4}+{\mathit{e}}_{\mathit{ij}}$$

(12)

$${\widehat{\mathit{h}}}_{\mathit{ij}}=1.3+{\mathit{b}}_0\mathit{H}{\mathit{o}}^{\left({\mathit{b}}_1+{\mathit{u}}_{\mathit{j}}\right)}{{\left(1-exp\left(-{\mathit{b}}_2\cdot {\mathit{d}}_{\mathit{ij}}\right)\right)}^{{{\mathit{b}}_3}^{\left(\raisebox{1ex}{$\mathit{N}$}\!\left/ \!\raisebox{-1ex}{${}_{{\mathit{d}}_{\mathit{g}}}$}\right.\right)}}}^{{\mathit{b}}_4}+{\mathit{e}}_{\mathit{ij}}$$

(13)

$${\widehat{\mathit{h}}}_{\mathit{ij}}=1.3+\mathit{H}{\mathit{o}}^{\left({\mathit{b}}_0+{\mathit{u}}_{\mathit{j}}\right)}{{\left(1-exp\left(\left(\raisebox{1ex}{${\mathit{b}}_1$}\!\left/ \!\raisebox{-1ex}{$100$}\right.\right)\cdot \mathit{HI}\cdot {\mathit{d}}_{\mathit{ij}}\right)\right)}^{{}^{\left(\raisebox{1ex}{$\mathit{N}$}\!\left/ \!\raisebox{-1ex}{${\mathit{d}}_{\mathit{g}}$}\right.\right)}}}^{{\mathit{b}}_2}+{\mathit{e}}_{\mathit{ij}}$$

(14)

$${\widehat{\mathit{h}}}_{\mathit{ij}}=1.3+\mathit{H}{\mathit{o}}^{\left({\mathit{b}}_0+{\mathit{u}}_{\mathit{j}}\right)}{{\left(1-exp\left(\left(\raisebox{1ex}{${\mathit{b}}_1$}\!\left/ \!\raisebox{-1ex}{$10$}\right.\right)\cdot \mathit{HI}\cdot {\mathit{d}}_{\mathit{ij}}\right)\right)}^{{}^{\left(\raisebox{1ex}{$\mathit{N}$}\!\left/ \!\raisebox{-1ex}{${\mathit{d}}_{\mathit{g}}$}\right.\right)}}}^{{\mathit{b}}_2}+{\mathit{e}}_{\mathit{ij}}$$

(15)

where b ₀ – b ₄ are the fixed parameters of the model (common to all plots); (u _j , v _j ) ~ N(0, τ) are the random parameters (specific to each plot); and ${\widehat{\mathit{h}}}_{\mathit{ij}}$ and e _ij are respectively the height and error estimated by the model for the i-th observation (tree) in the j-th plot.

The values of the parameters and goodness-of-fit statistics for the local mixed models (Eqs 10 and 11) and for the generalized mixed models (Eqs 12 to 15) are shown in Tables 3 and 4, respectively. We compared the RMSE values obtained with the mixed effects equations with those obtained with fixed effects equations (fitted by ONLS); the values obtained with the local mixed model (Eq. 10) and the generalized mixed model (Eq. 12) for the Pinus grouping were 25.0% and 5.2% lower than those obtained with the local model (Eq. 7) and the generalized model (Eq. 8) without random parameters, respectively. For the group of Quercus species, the RMSE values obtained with the local mixed model (Eq. 11) and the generalized mixed model (Eq. 13) were 26.0% and 9.0% lower than those obtained with the local (Eq. 7) and the generalized models without random parameters (Eq. 8), respectively. For the other conifers, the RMSE values obtained with the local mixed model (Eq. 11) and the generalized mixed model (Eq. 14), were 9.8% and 14.3% lower than those obtained with Eqs (7) and (9), respectively. For the group comprising other broadleaf species, the RMSE values were 20.9% and 25.0% lower with the local mixed model (Eq. 11) and the generalized mixed model (Eq. 15) than with Eqs (7) and (9). The results obtained for BIC and R² were similar to those obtained for RMSE.

On inspecting the graphs of the residuals for the heights estimated by the models for each species grouping, we did not find any anomalies that would suggest non compliance of underlying hypothesis of independence of errors or homogeneity of variance. The magnitude of the bias in the residual values estimated by the two fitting methods (ONLS and NLMIXED) was consistent for all ranges and classes of heights observed by the defined species groupings.

Calibrated response

Calibration option CR1 for Eq. (12) in Pinus species and Eq. (13) in Quercus species was the most accurate when the total height of a subsample of 3 trees close (± 10%) to the mean breast height diameter for the plot was measured (Figure 2), as indicated by the slight decrease in the RMSE by 0.4% and 3.0% respectively for the two groups relative to the values estimated by the generalized model (Eq. 8) fitted without random parameters. For the other broadleaf species, calibration option CR2 for Eq. (14) was the most accurate (in terms of RMSE) when a subsample of 3 trees of mean, minimum and maximum breast height diameter were measured in each plot, as the RMSE was 21% lower than that estimated with the generalized model fitted without random parameters (Eq. 9). Finally, in the group of other conifer species, the RMSE value obtained with Eq. (15) and calibration option CR2 was 13% lower than that obtained with the generalized model fitted without random parameters (Eq. 9). Vargas-Larreta et al. (2009) found that the decrease in RMSE varied between 3.7 and 13.3% on calibrating the model of Sharma and Parton (2007) with data from 1 tree selected at random in the plot. However, Calama and Montero (2004) observed that use of a subsample of the 5 trees with the largest diameters significantly decreased the RMSE value; Castedo Dorado et al. (2006) observed that the bias estimated with Schnute’s generalized equation was lower when a subsample of the 3 trees with the smallest diameters was used for calibration than when the estimate was made without random parameters.

Figure 2

Value of RMSE for the calibrated equations (CR1 and CR2) fitted by NLMIXED and ONLS for the four groups of species shown.

The variation in the value of RMSE with respect to the number of trees used with the two calibration options for the four main groups of species studied is shown in Figure 2. This statistic was also compared with those values obtained when fitting the equations by the NLMIXED (minimum value of RMSE reached only using all trees as a calibration subsample) and ONLS (maximum value of RMSE using only fixed parameters) methods.

In the calibration process, the reduction in the RMSE value was particularly evident with the generalized mixed models for the other broadleaf species and other conifer species (21.0% and 13.0% respectively) compared with the generalized model fitted without random parameters; however, for the Pinus and Quercus groupings, the decrease in the value of this statistic was lower. Both calibration options resulted in an important reduction of RMSE for the local mixed model compared to the same model fitted without random parameters for all the groups analized. In accordance with Trincado et al. (2007), the use of a local mixed model in forest inventories with a subsample of trees to calibrate and then predict the total height of all trees not used in calibration allows retention of a simple model structure (i.e. without the need to include stand predictor variables) and may be an useful alternative to generalized mixed models when there is a lack of data to calculate stand variables.

Conclusions

Two generalized equations (Eqs 8 and 9) were derived from a local equation (Eq. 7) and used to estimate total tree height from breast height diameter and stand variables for the 25 species identified in the sample by using mixed models. The variability between plots is explained in terms of the random effect of each plot and from the stand variables included in the generalized models.

For species in the Pinus and Quercus groups, inclusion of the height measurements of 3 trees close (± 10%) of the mean breast height diameter from each plot improved the predictive capacity of the calibrated model. For the species included in other broadleaf species and other conifer species, the predictive capacity of the model was improved by including the total height measured in a subsample of 3 trees of minimum, mean and maximum breast height diameter. The possibility of using complementary data from the stands to calibrate the mixed models provides a clear advantage over models developed by other procedures, which require large amounts of data or are less accurate.