Can random components explain differences in the height–diameter relationship in mixed uneven-aged stands?

Abstract

• Context

Tree height prediction is an important issue in forest management since tree heights are usually measured only in a sample of trees. Although numerous model approaches have been used for this purpose, no agreement on which one is more appropriate has been achieved.

• Aims

To analyse the random effects of basic and generalised height–diameter (h–d) models fitted to multi-species uneven-aged forest stands, and to establish their ability to explain differences between ecoregions, plots and species.

• Methods

Height and diameter measurements for 29,084 trees from 187 sample plots located in the state of Durango (Mexico) were used. Basic and generalised h–d models were fitted in a mixed-models framework. The variability between ecoregions, plots and species was considered in the random effects definition. Model calibration for different height sampling designs and sampling sizes was also analysed.

• Results

Random components performed well in explaining the differences in the h–d relationship between the different plots and species; however, no significant variance for the random effects was found for the different ecoregions. A calibrated basic h–d model produced similar results to a fixed-effects generalised h–d model when a sufficiently large number of trees was used in the calibration process.

• Conclusion

From a practical point of view, if no calibration is carried out, different models should be used for the different species, so that at least the variation among species is captured.

1 Introduction

The forests of Durango (Mexico) cover an area of 4.9 million hectares, and timber resources account for about one quarter, in terms of volume, of the national forest resources in Mexico. Durango State is the top producer of pine growing stock (27.7 %), coniferous roundwood (35.3 %) and other roundwood (37.6 %) nationwide (SRNyMA 2006). The pine-oak forests of Durango have been ranked as the most important forests in Mexico because of their extent and economic value. These forests are of particular interest, not only because they represent a unique ecosystem, but also because they are owned and managed by local communities, known as Ejidos (Thoms and Betters 1998).

Most of the forests in Durango are irregular, and conifer species occur as mixtures with hardwood species. This irregularity refers to the spatial arrangement of trees (vertical and horizontal irregularity) and the variation in the age structure of trees and stands. This structure is the result of the management history, which has depended on land ownership, as well as the economic and social changes that have taken place in the state, and also on natural conditions (Wehenkel et al. 2011).

The management of such mixed uneven-aged forests is more complex than the management of even-aged ones, and one of the most pressing research problems is the development of growth models to define sustainable harvests. However, few studies have attempted to estimate tree growth and stand development in Durango (e.g. Corral-Rivas et al. 2004; Vargas-Larreta et al. 2009).

In this sense, height–diameter (h–d) models (i.e. models relating individual tree height with individual tree diameter) are very useful for stand growth dynamics, yield, site index and dominant height estimation (Curtis 1967), stand structural analysis (Morrison et al. 1992), damage appraisal and stand stability (Parresol 1992), and product recovery and carbon budgeting models (Newton and Amponsah 2007).

In general, most height–diameter (h–d) models have been applied to pure even-aged stands or plantations (e.g. Soares and Tomé 2002; López Sánchez et al. 2003). Despite the homogeneous characteristics of this type of forest, a single h–d model is not usually adequate for all possible situations within a stand because the height curve does not remain constant and increases in steepness with age. Height curves for good quality sites have steeper slopes than those for poor quality sites, and for a particular height, trees growing in high density stands will have smaller diameters than those growing in less dense stands (Curtis 1967; Bailey and Brooks 1994; Lappi 1997; Zhang et al. 1997). This is even more evident in mixed and uneven-aged stands, in which different species, ages, sizes, crown types and levels of shade tolerance coexist (Vargas-Larreta et al. 2009).

To account for this level of variance, basic h–d relationships can be improved by taking into account stand variables that introduce the dynamics of each stand into the model (Curtis 1967). Such models, called generalised h–d models, have been shown to be applicable to both even-aged (e.g. Soares and Tomé 2002; Castedo-Dorado et al. 2006) and uneven-aged stands (e.g. Sharma and Zhang 2004; Sharma and Parton 2007), and even to a mixture of both types of stands (e.g. Crecente-Campo et al. 2010). These models use diameter and stand-specific variables as regressors, accounting for the differences in the h–d relationship both within stands and over time (Curtis 1967; Soares and Tomé 2002; Sharma and Zhang 2004).

The hierarchical structure of the h–d data (i.e. trees within plots and plots within stands) usually results in a lack of independence among measurements, since observations from the same sampling unit are highly correlated (West et al. 1984). Development of the mixed-modelling methodology provided a statistical method capable of explicitly modelling this nested stochastic structure (Lappi 1997; Calama and Montero 2004; Castedo-Dorado et al. 2006). Mixed models are composed of a fixed functional part, common to the population, and random components acting at each sampling level. These models allow description of the variability in given phenomena among different sampling levels after defining a common fixed functional structure and a covariance structure for the random effects and residual terms (Lindstrom and Bates 1990; Calama and Montero 2005). Detailed information on nonlinear mixed-effects modelling for h–d relationships are provided by Calama and Montero (2004). General information and discussion on non-linear mixed models in the forestry context can be found in Hall and Bailey (2001). Finally, the multilevel case of mixed models has been discussed by several authors in a general context (e.g. Lindstrom and Bates 1990; Longford 1993; Goldstein 1995).

Mixed models may also improve the predictions obtained if it is possible to estimate the value of the random effects for an individual that has not been sampled. This approach is known as localisation or calibration and can be applied if supplementary observations of the dependent variable (total tree height in this case) are available (Lappi 1991; Jayaraman and Lappi 2001; Lynch et al. 2005; Calama and Montero 2005).

The h–d relationship is an important component of growth models. This model plus a compatible volume system (composed of a taper equation and a disaggregation equation) allows for stand volume classification by merchantable sizes, which are important tools for sustainable management of the species in the study area.

Generalised h–d models for uneven-aged stands were developed in a previous study for the region of El Salto (Durango, Mexico) (Vargas-Larreta et al. 2009). However, this previous study was very local, and results cannot be extrapolated to the state of Durango. Moreover, the authors did not consider the different hierarchical levels (ecoregions, plots and species) in the random effects definition, and comparisons with non-generalised models were not made.

The objective of the present study was to fit basic and generalised h–d models in a mixed-models framework and to analyse the random effects estimates of such models, in order to establish their ability to explain the differences in the h–d relationships between ecoregions, plots and species observed in mixed uneven-aged stands in Durango, México. Models for practical use are planned for development on the basis of these analyses; therefore, the impact on predictability of the model of having a subsample of complementary trees where both variables (diameter and height) are measured was also analysed for different height sampling designs and sampling sizes within each sampling level.

2 Materials and methods

2.1 Study area and data description

The forests of Durango State are located in the Sierra Madre Occidental, at altitudes ranging between 363 and 3,190 m above sea level. The climate in the area is temperate to tropical, with most rainfall occurring in the summer, annual average precipitation of 443–1,452 mm and annual average temperature of 8.2–26.2 °C (Silva-Flores et al., submitted). The forests are mainly managed by selective removals, with only a small number of shelterwood harvests (3 % of the productive forest area). Clearfellings, which require special management expertise, are rarely applied (Wehenkel et al. 2011).

The data used in the study were obtained from a network of permanent sample plots used to monitor the growth and yield of Durango’s forests. The plots used in this study were established between 2007 and 2010 and cover the main forest types and the current diameter distributions of commercial forests in Durango. The plots are 50 × 50 m in size and are distributed by systematic sampling (with some exceptions), with a variable grid ranging from 3 to 5 km, depending on the size of the Ejidos. The sampling plots are planned for re-measurement at 5-year intervals. Among other variables, tag number, species code, breast height diameter (d, centimetres, 1.3 m above ground level), total tree height (h, metres), height to the live crown (metres), azimuth (degrees Centigrade) and radius (metres) from the centre of the plot of all trees equal or larger than 7.5 cm in diameter were recorded. Currently, the database includes measurement data for 29,084 trees from 187 sample plots.

These sample plots were located between latitudes 23° 10′ N and 25° 20′ N, and between longitudes 104° 45′ W and 106° 40′ W (Fig. 1). They were located at altitudes ranging from 2,006 to 2,939 m above sea level. The soils were very variable, with sand (particles over 2 mm) percentages varying from 18 to 85 % (average value of 52 %), silt (particles between 0.2 and 2 mm) percentages varying from 1 to 87 % (average value of 25 %) and clay (particles below 2 mm) percentages varying from 11 to 42 % (average value of 23 %). The number of stems per plot varied from 16 to 498 (average value of 156), and the number of species per plot varied from 2 to 12 (average value of 7).

Fig. 1

Location, in the state of Durango (Mexico), of the research plots installed in the Ejidos in the four ecoregions analysed in this study: I (El Salto), II (San Dimas), III (Santiago y Topia) and IV (Tepehuanes)

Four different ecoregions, designated ecoregions I to IV, were included in the dataset: I (El Salto), II (San Dimas), III (Santiago y Topia) and IV (Tepehuanes) (Fig. 1). These ecoregions were defined on the basis of regional organisations previously designated by the National Forest Commission (CONAFOR), as Regional Forest Management Units (Unidades de Manejo Forestal Regional (UMAFORs)), which tend to focus on environmental services as a public policy goal. UMAFORs are defined as “representing” 30 % of the forest owners within a particular region, usually delineated along ecological zones, and thus include cross-community interdependence in a biological sense (Camille and García-López 2008).

Sixteen different species or groups of species (hereafter simply referred as species) were distinguished for posterior analysis as they presented similar growth patterns:

• 1: Pinus arizonica

• 2: Pinus ayacahuite

• 3: Pinus cooperi

• 4: Pinus durangensis

• 5: Pinus engelmannii

• 6: Pinus herrerae

• 7: Pinus lumholtzii

• 8: Pinus teocote

• 9: other pines: Pinus chihuahuana, Pinus douglasiana, Pinus michoacana, Pinus oocarpa, Pinus tenuifolia

• 10: Quercus sideroxyla

• 11: other Quercus: Quercus candicans, Quercus coccolobifolia, Quercus conzattii, Quercus cordifolia, Quercus crassifolia, Quercus durifolia, Quercus fulva, Quercus grisea, Quercus microphylla, Quercus obtusata, Quercus resinosa, Quercus rugosa, Quercus rysophylla, Quercus urbanii

• 12: Abies durangensis, Picea chihuahuana, Pseudotsuga menziesii

• 13: Cupresus lindleyi

• 14: Juniperus depeana

• 15: Alnus acuminata

• 16: Arbutus xalapensis

The following stand variables were calculated from the tree data recorded in each plot: stems per hectare (N, trees per hectare), stand basal area (G, square metres per hectare), quadratic mean diameter (dg, centimetres), dominant height (hdom, metres) and dominant diameter (ddom, centimetres). The last two variables were calculated from the proportion of the 100 trees with the largest diameter per hectare. Some of the trees with no species code, and those trees with broken or dead tops were then excluded from further analysis. Data from trees with broken or dead tops were not used to calculate hdom.

The scatter plots of total tree height against diameter at breast height for each species were visually examined to detect possible anomalies in the data. Extreme data points (i.e. outside the overall picture of the h–d pairs) were not observed. Some unusual data points (i.e. trees with quite large height and small diameters) were observed, but since we are attempting to develop a model for uneven-aged mixed forest (with a high variability), the unusual data points may be very valuable signals on h–d relationships rather than outliers. We therefore decided to keep these points in model fitting and to analyse them via the residuals.

The database (Table 1) was finally randomly divided into a fitting part (80 % of the plots) and an evaluation part (20 % of the plots).

Table 1 Summary statistics for the fitting and the evaluation datasets

2.2 Models analysed

In a first step, 27 basic h–d models (h as a function of only d) from Huang et al. (2000) were fitted to the dataset. In a second step, a total of 30 generalised h–d models selected from previous studies (López Sánchez et al. 2003; Sharma and Zhang 2004; Castedo-Dorado et al. 2006; Sharma and Parton 2007) were fitted to the dataset. Some modifications to the models were also tested (i.e. dg and mean height were replaced by ddom and hdom, respectively). Models that require age or site index were not used or were modified to exclude these variables, as the stands were uneven-aged.

2.3 Model fitting and comparison

To make an initial selection from the above-mentioned models, they were fitted for the different species, by the ordinary nonlinear least squares (ONLS) method with the MODEL procedure of SAS/ETS® (SAS Institute Inc. 2008). Statistical and graphical analyses were used to compare the performance of the models. Four statistical criteria obtained from the residuals were examined: the root mean square error (RMSE), the model efficiency (EF) (similar to the coefficient of determination for linear regression), the mean bias (E) and Akaike’s information criterion (AIC) (Akaike 1974), summarised as follows:

$$ \mathrm{RMSE}=\sqrt{\frac{{\displaystyle {\sum}_{l=1}^{l=n}{\left({y}_l-{\widehat{y}}_l\right)}^2}}{n-p}} $$

(1)

$$ EF=1-\frac{{\displaystyle {\sum}_{l=1}^{l=n}{\left({y}_l-{\widehat{y}}_l\right)}^2}}{{\displaystyle {\sum}_{l=1}^{l=n}{\left({y}_l-\overline{y}\right)}^2}} $$

(2)

$$ E=\frac{{\displaystyle {\sum}_{l=1}^{l=n}\left({y}_l-{\widehat{y}}_l\right)}}{n} $$

(3)

$$ \mathrm{AIC}=n \ln \left({\displaystyle {\sum}_{l=1}^{l=n}{\left({y}_l-{\widehat{y}}_l\right)}^2}/n\right)+2\left(p+1\right) $$

(4)

where y _l , ${\widehat{y}}_l$ and $\overline{y}$ are the measured, estimated and average values of the dependent variable for the lth tree, respectively; n is the total number of observations used; p is the number of model parameters; and ln is the natural logarithm. Subscript l was used to refer to tree for concordance with the mixed-models notation used later.

Model validation was not carried out because an independent dataset was not available, and other validation methods seldom provide any additional information compared with the fitting statistics (Kozak and Kozak 2003). However, the evaluation data (Table 1) were used in model calibration in order to provide informative results about the calibration process (see details later).

Once the best basic and generalised h–d models were selected, a nonlinear mixed-effects modelling framework was used (e.g. Lappi 1997; Calama and Montero 2004; Castedo-Dorado et al. 2006), because there was high variability between ecoregions, plots and species. Basically, the parameter vector of the nonlinear model can be defined (Pinheiro and Bates 2000) as:

$$ {\boldsymbol{\Phi}}_{ijk}={\mathbf{A}}_{ijk}\boldsymbol{\uplambda} +{\mathbf{B}}_{ijk}{\mathbf{b}}_{ijk} $$

(5)

where λ is the p × 1 vector of fixed population parameters (where p is the number of fixed parameters in the model), b _ijk is the q × 1 vector of random effects associated with the ith ecoregion, the jth plot and the kth species (where q is the number of random effects in the model), and A _ijk and B _ijk are design matrices of size r × p and r × q (where r is the total number of parameters in the model) for the fixed parameters and random effects specific to each ecoregion, plot and species, respectively.

The SAS macro NLINMIX (Littell et al. 2006) was used to fit the models. This macro implements a two-step process for estimating the parameters. First, it performs a standard nonlinear regression by least squares, and in a second step it performs successive calls to the procedure of estimating parameters of linear mixed models on a linearised approximation function until convergence is achieved. The macro can use two methods of expansion: expansion around zero or the best linear unbiased predictor (BLUP), which is the expected value for the random effects (Beal and Sheiner 1982), or expansion around the empirical best linear unbiased predictor (EBLUP) of the random effects (Littell et al. 2006). Both approaches produce reliable estimates (Davidian and Giltinan 1993; Pinheiro and Bates 2000). The EBLUP method sometimes generates slightly better results but requires more computation time, is less stable, and is very sensitive to model specification (Hartford and Davidian 2000). Moreover, Vonesh (1996) showed that consistent estimates with the EBLUP expansion method could only be obtained if the number of individuals and observations per individual are infinite. For our data, most of the models with various combinations of random effects failed to converge when the EBLUP expansion method was used. Therefore, all results presented in this study are based on the expansion around zero method.

Different combinations of parameters were assumed to be mixed (composed of a fixed part, and a random part, specific to each ecoregion, plot and species). A variance components structure was selected for the covariance structure of the random effects, since convergence was easily achieved and better fitting statistics were obtained than with other random effects variance structures. This matrix contains the variance components in a diagonal structure. The significance of the variance of the random effects was tested using the COVTEST option in the NLINMIX macro, which displays asymptotic standard errors and Wald tests for covariance parameters (SAS Institute Inc. 2009).

Finally, heterocedasticity was analysed via the residuals. If necessary, a special structure for the within-plot variance–covariance matrix could be used to include weighting factors to balance error variance (Calama and Montero 2004).

2.4 Model calibration

An advantage of mixed-effects models is that, if a subsample of m tree heights is available, such data can be used to predict the random effects vector b _ijk (i.e. to calibrate the model), with the following expression (Vonesh and Chinchilli 1997):

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

(6)

where $\widehat{\mathbf{D}}$ is a q × q variance–covariance matrix for the among-ecoregion, among-plot and among-species variability, common to all plots and estimated in the general fitting of the model; ${\widehat{\mathbf{R}}}_{ijk}$ is the m × m variance–covariance matrix for within-plot variability; ${\widehat{\mathbf{e}}}_{ijkl}$ is the residual vector m × 1, the components of which are given by the difference between the observed height of each tree included in the subsample, and the value predicted by the model including only fixed effects and ${\widehat{\mathbf{Z}}}_{ijk}$ is the m × q matrix of partial derivatives with respect to the random effects evaluated at ${\widehat{\mathbf{b}}}_{ijk}$.

For the best basic and generalised h-d models, the calibrated response was evaluated using the evaluation dataset for different height sampling designs and sampling sizes within each ecoregion, plot and species. This dataset was used to show the calibrated results for plots that were not included in the fitting process. The alternatives selected according to previous studies (Calama and Montero 2004; Castedo-Dorado et al. 2006; Crecente-Campo et al. 2010) were:

• 1. Total height of one to four randomly selected trees per species and plot.

• 2. Total height of the one to four largest trees per species and plot.

• 3. Total height of the one to four smallest trees per species and plot.

• 4. Total height of the one to four medium-size trees per species and plot.

• 5. Total height of three trees per species and plot: the smallest, the largest and the medium-size trees.

Some combinations of plot and species with less than four trees were, therefore, eliminated, reducing the evaluation dataset from 6,345 to 6,225 trees in 38 plots, and resulting in a modal value of five species per plot.

The five calibration alternatives were evaluated in terms of the previously defined statistics (RMSE, EF and E) and compared with the estimates obtained by ONLS in the individual fit of each equation to each of the ecoregion, plot and species groups. For the randomly selected trees, mean values of the statistics after 100 simulations were obtained, as this was considered a large enough value to obtain representative results and was in accordance with previous studies (e.g. Castedo-Dorado et al. 2006; Crecente-Campo et al. 2010).

Finally, the selected models were also fitted independently for each species, and random effects specific to each ecoregion and plot were added, to determine whether the consideration of species as a fixed effect led to a significant improvement in the accuracy of the predictions.

3 Results

3.1 ONLS fits

Statistics for the fitting of the basic h–d equations (not shown) showed that the models of Bertalanffy–Richards (Bertalanffy 1949; Richards 1959), Weibull (Yang et al. 1978) and Schnute (Schnute 1981) were the most accurate, with very similar values for the fitting statistics. When the models were fitted to each species, the model of Schnute (1981) (Eq. 7) produced slightly better results (lower average AIC and RMSE values) and was, therefore, selected as a basic h–d model for posterior analysis:

$$ {h}_l={\left({1.3}^{a_0}+\left({a}_1{}^{a_0}-{1.3}^{a_0}\right)\frac{1- \exp \left(-{a}_2\cdot {d}_l\right)}{1- \exp \left(-{a}_2\cdot 100\right)}\right)}^{\frac{1}{a_0}} $$

(7)

where h _l and d _l are, respectively, the total height and diameter at breast height of the lth tree and a _i are the parameters to be estimated.

From the total number of generalised h–d models analysed, a modification of the Bertalanffy–Richards (Bertalanffy 1949; Richards 1959) model (Eq. 8) showed the best results, as inferred from the fit statistics (Table 2). The same equation has previously been used by Sharma and Parton (2007) for several species in northern Ontario (Canada) and by Vargas-Larreta et al. (2009) for several species in El Salto (Durango, Mexico).

Table 2 Statistics for the fitting of the basic and generalised mixed h–d models

$$ {h}_l={\left(1.3+{a}_0\cdot ddo{m}_j{}^{a_1}\left(1- \exp \left(-{a}_2\cdot {\left(\frac{N_j}{G_j}\right)}^{a_3}{d}_l\right)\right)\right)}^{a_4} $$

(8)

where subscript j refers to plot and subscript l refers to tree; the other variables are as previously defined.

3.2 Basic h–d mixed model

Independent fitting of Eq. 7 to each ecoregion, plot and species, by the ONLS method, showed that most of the variability in parameters was derived from the differences between plots and, secondly, species. The a ₀ parameter showed the higher coefficient of variation for the different plots, while the a ₂ parameter showed the higher coefficient of variation among species. Variability among ecoregions was small, as inferred from the differences between the values of the parameters for the different ecoregions. Finally, all possible combinations of random effects for each ecoregion, plot and species were tried. However, no significant variance for the ecoregion-specific random effects was found when fitting the model. The best results in the fitting of Eq. 7 as a mixed model were obtained when all the parameters were expanded to include plot-specific random effects, and some parameters were expanded to include species-specific random effects, resulting in the following mixed model (once parameters were scaled so that they were all around the same order of magnitude, to avoid instabilities in the fitting algorithm; Littell et al. 2006):

$$ {h}_{jkl}={\left({1.3}^{a_0+{v}_{0j}+{w}_{0k}}+\left({\left[20\left({a}_1+{v}_{1j}+{w}_{1k}\right)\right]}^{a_0+{v}_{0j}+{w}_{0k}}-{1.3}^{a_0+{v}_{0j}+{w}_{0k}}\right)\frac{1- \exp \left(-\frac{a_2+{v}_{2j}}{50}\cdot {d}_{jkl}\right)}{1- \exp \left(-\frac{a_2+{v}_{2j}}{50}\cdot 100\right)}\right)}^{\frac{1}{a_0+{v}_{0j}+{w}_{0k}}} $$

(9)

where subscript j refers to plot, subscript k refers to species and subscript l refers to tree; v _ij and w _ik are the random effects specific to each plot and species, respectively, and the other variables are as previously defined.

Fit statistics for Eq. 9 (Table 2) showed an improvement over the basic ONLS model (Eq. 7). RMSE values were reduced by approximately 32 %, but E was greatly increased, although it kept at a reasonable level (approximately 0.5 m). This may suggest that differences have become biased either negative or positive. However, the plot of raw residuals against predicted values (Fig. 2a) showed no systematics trends and a homogeneous distribution. Trends in the residuals usually relate to missing covariates or heterogeneity of variance. Parameter estimates for Eq. 9 are shown in Table 3.

Fig. 2

Plots of raw residuals versus predicted values for Eqs. 9 (a) and 10 (b)

Table 3 Parameter estimates for the basic h–d models selected

3.3 Generalised h–d mixed model

The best results in the fitting of Eq. 8 as a mixed model were obtained when two parameters were expanded to include both plot and species-specific random effects. As previously, the ecoregion random effects variance was not significant. The resultant model (once parameters were scaled so that they were all around the same order of magnitude, to avoid instabilities in the fitting algorithm: Littell et al. 2006) was:

$$ {h}_{jkl}={\left(1.3+\left({a}_0+{v}_{0j}+{w}_{0k}\right)\cdot ddo{m}_j{}^{a_1}\left(1- \exp \left(-\frac{a_2}{50}\cdot {\left(\frac{N_j}{G_j}\right)}^{\frac{a_3}{4}}{d}_{jkl}\right)\right)\right)}^{a_4+{v}_{4j}+{w}_{4k}} $$

(10)

where all the variables are as previously defined.

Fit statistics for Eq. 10 (Table 2) showed a reduction in RMSE by approximately 29 % as compared with the generalised ONLS model (Eq. 8). However, an increase in E was obtained, although it kept at a reasonable level (approximately 0.14 m). The plot of raw residuals against predicted values (Fig. 2b) showed no systematics trends and a homogeneous distribution. Parameter estimates for Eq. 10 are shown in Table 4.

Table 4 Parameter estimates for the generalised h–d models selected (all parameters were significant at p < 0.05; standard errors in brackets)

While RMSE was only reduced by 11 % between Eqs. 9 and 10, E was reduced by 73 %. The differences between Eqs. 9 and 10 were also visually examined by species (Fig. 3). As confirmed, RMSEs were very similar, with only a slightly lower RMSE for Eq. 10. The greater RMSE differences were observed for species 9 and 12. E, however, was greater for all species for Eq. 9, particularly for species 6 and 12.

Fig. 3

Plots of the fitting statistics for Eqs. 9, 10, 11 and 12, for each species. Lines are used only for a better comparison of the differences between the models. Numbers on the X-axis refer to the species groups outlined in the text

3.4 Fits to individual species

In an attempt to reduce the sources of variation, the variation between species was treated as a fixed effect (i.e. by fitting the models independently for each species and adding random effects for each plot and ecoregion), producing Eqs. 11 and 12.

$$ {h}_{jkl}={\left({1.3}^{a_0{}_k+{v}_{0 jk}}+\left({\left[20\left({a}_{1k}+{v}_{1 jk}\right)\right]}^{a_0{}_k+{v}_{0 jk}}-{1.3}^{a_0{}_k+{v}_{0 jk}}\right)\frac{1- \exp \left(-\frac{a_{2k}+{v}_{2 jk}}{50}\cdot {d}_{jl}\right)}{1- \exp \left(-\frac{a_{2k}+{v}_{2 jk}}{50}\cdot 100\right)}\right)}^{\frac{1}{a_0{}_k+{v}_{0 jk}}} $$

(11)

$$ {h}_{jkl}={\left(1.3+\left({a}_{0\kern0.5em k}+{v}_{0 jk}\right)\cdot ddo{m}_j{}^{a_{1\kern0.5em k}+{v}_{1 jk}}\left(1- \exp \left(-\frac{a_{2\kern0.5em k}}{50}\cdot {\left(\frac{N_j}{G_j}\right)}^{\frac{a_{3\kern0.5em k}}{4}}{d}_{jl}\right)\right)\right)}^{a_{4\kern0.5em k}+{v}_{4 jk}} $$

(12)

where subscript j refers to plot, subscript k refers to species and subscript l refers to tree; a _{i k} are the fixed-effects parameters for the kth species, v _{i jk} are the random effects specific to each jth plot and kth species and the other variables are as previously defined.

In this case, the estimate of the variance of the ecoregion random effects was again not significant. Fitting statistics were clearly superior to those obtained with the previous models, with a large reduction in E and a small reduction in RMSE (Table 2). Results for Eq. 12 were again only slightly better than those for Eq. 11 (Fig. 3). Parameter estimates for Eq. 11 are shown in Table 3, and those for Eq. 12 are shown in Table 4. In Eq. 11, the a _{1 k} parameter (the model asymptote) showed clearly lower values for species 10, 11, 14, 15 and 16, which are the most shade tolerant species.

3.5 Model calibration

Results for the calibration process for all the mixed-models using the evaluation dataset showed small differences between Eqs. 9 and 10 and between Eqs. 11 and 12, using the random-trees selection method for calibration (Fig. 4). The results for all species together (Fig. 5) showed a similar trend, with quite different RMSEs for the fixed-effects model (i.e. the mixed model with 0 trees used in calibration) and quite similar RMSEs for the mixed-effects model using four trees for calibration. Similar results were found for EF and E (Fig. 5). The greatest reductions in RMSE between a fixed-effects model and a calibrated mixed-effects model were observed for the basic h–d models, i.e. for Eqs. 9 and 11. As the evaluation dataset has a modal value of five species per plot, the calibrated results were aggregated in Figs. 4 and 5. So calibration with one tree per species means, a measurement of five trees per plot; calibration with two trees per species means a measurement of 5 × 2 = 10 trees per plot, and so on.

Fig. 4

Plots of the statistics obtained in the calibration process including the random tree selection method, for Eq. 9 (first row), Eq. 10 (second row), Eq. 11 (third row) and Eq. 12 (fourth row), for each species. Zero trees for calibration represent the fixed-effects model. The X-axis shows aggregate values to an average number of trees per plot, because the evaluation dataset has a modal value of five species per plot. Numbers in the legend refer to the species groups outlined in the text

Fig. 5

Plots of the statistics obtained in the calibration process including the random tree selection method, for Eqs. 9, 10, 11 and 12, for all species together. Zero trees for calibration represent the fixed-effects model. The X-axis shows aggregate values to an average number of trees per plot, because the evaluation dataset has a modal value of five species per plot

In a practical sense, when using Eq. 9, at least two trees of each species (i.e. a modal value of 10 trees per plot in the evaluation dataset) must be measured to calibrate the model and to obtain the same results as the uncalibrated Eq. 12 (Fig. 5). For Eq. 11, the measurement of only one tree of each species per plot (a modal value of 5 trees per plot in the evaluation dataset) is enough to calibrate the model and obtain the same results as the uncalibrated Eq. 12. These general results are, however, influenced by the fact that Eq. 12 showed very poor results in the calibration process for species 12 and 13 (Fig. 4), maybe because there was a low number of observations included in the evaluation dataset (55 and 15 observations for each species, respectively) for these two species. For the rest of the species, Eq. 12 showed better calibrated results than Eq. 11 (Fig. 4).

The results obtained for the different calibration methods with Eq. 11 are shown in Fig. 6. This equation was selected as a reference for testing the different calibration methods because it produced comparable results to Eq. 12 for the total data set (Fig. 5) and needs fewer measurements to be applied. The medium-size tree selection method and the random tree selection method performed better than others in terms of RMSE and EF, followed very closely by the all-sizes tree selection method, which was superior in terms of E (Fig. 6). The other tree selection methods produced higher errors than the above-mentioned ones (Fig. 6).

Fig. 6

Plots of the fitting statistics obtained with the different tree selection options for calibration of Eq. 11. The X-axis shows aggregate values to an average number of trees per plot, because the evaluation dataset has a modal value of five species per plot. Lines and points in the legend refer to size of trees used in calibration

4 Discussion

Several basic h–d models were tested in this study for modelling the h–d relationship in mixed uneven-aged stands in Durango (Mexico). Confirming the findings of Zhang (1997), the best models in this study were the Bertalanffy–Richards (Bertalanffy 1949; Richards 1959), Weibull (Yang et al. 1978) and Schnute (Schnute 1981) models. According to Lei and Parresol (2001), the Schnute function together with the Bertalanffy–Richards function are probably the most flexible and versatile functions available for modelling height–diameter relationships. Finally, the Schnute’s equation was selected for fitting as a basic h-d mixed model because it showed superior fit statistics for almost all species groups.

From a total of 25 generalised h–d models tested in this study, the equation that provided the best fit to the dataset was derived from the Bertalanffy–Richards model and included diameter at breast height, dominant height, stand basal area and trees per hectare in its formulation. The relationship of trees per hectare and dominant height to height growth has already been discussed by several authors (e.g. Calama and Montero 2004; Saunders and Wagner 2008; Vanclay 2009). Stand density (measured for example by stand basal area or trees per hectare) is the most important factor that affects the h–d relationship in a stand (Zhang et al. 1997; Saunders and Wagner 2008); in dense stands, trees of the same diameter are usually taller than those in less dense stands. Dominant height has also previously been included in generalised h–d models by many authors (e.g. Castedo-Dorado et al. 2006; Vargas-Larreta et al. 2009; Amaral Paulo et al. 2011). In the present study, models that included dominant height provided more accurate results than those including mean height. This may be advantageous, since fewer trees need to be measured to estimate dominant height than to estimate mean height, and this greater sampling effort may limit the future use of models using mean height (López Sánchez et al. 2003). This generalised h–d model showed better fitting statistics than the basic h–d model (Table 2).

However, when the basic and generalised h–d models were fitted in a mixed-models framework, fit statistics differed little between them (Fig. 3, Table 2), confirming the expectation that random effects specific to each plot and species will capture most of the observed variability between these sampling units, and that the incorporation of additional predictor variables (i.e. stand variables) in a basic h–d model has a major effect on the ability of the fixed-effects model in explaining between-individual variability, but not on the mixed-effects model (Trincado et al. 2007). Another advantage usually assigned to mixed-models is that they take into account the correlations between observations that belong to the same sampling unit, thus dealing with the dependence of error terms which should be avoided in regression models (Calama and Montero 2004; Castedo-Dorado et al. 2006).

A result that seems somewhat counter intuitive is that bias increased for Eqs. 9 and 10 as compared with the respective ONLS Eqs. 7 and 8. This may be attributed to the fact that Eqs. 7 and 8 do not make any distinction between species and all the observations have the same importance in the calculation of fitting statistics. However, Eqs. 9 and 10 differentiate between species, assigning a random effect for each one, and in some cases, as with species 15 and 16, the explained variance was very low (Fig. 3) and the associated bias very high. The same happens for the bias obtained for species 6, 12 and 13 for Eq. 9, which showed large values, thereby increasing the overall bias (calculated over all species without normalising by the number of observations for each species).

Expansion of fixed parameters to account for differences between ecoregions led to non-significant variance of the random effects. This is not surprising, since only four ecoregions were sampled in this study. Several t tests were performed to compare the average values of the height-diameter ratios for the different ecoregions. These average values (i.e. 0.72, 0.68, 0.76 and 0.56 for ecoregions 1, 2, 3 and 4, respectively) were significantly different between the four ecoregions. However, group-level variance estimates of zero often arise when fitting multilevel or hierarchical linear models, especially when the number of groups is small (Chung et al. 2013). The same may be expected for non-linear models.

Equations 11 and 12, i.e. the species-specific ones, showed better precision (RMSE) and bias (E) than Eqs. 9 and 10, respectively (Table 2, Fig. 3), indicating that a general mixed-model for the whole population may not be sufficient to explain the observed variation among species. However, the differences between the models were greatly reduced when they were calibrated with only one tree by species. Only a marginal gain in accuracy was observed when two, three or four sample trees per species were used for calibration (Fig. 5). Other studies have reported similar behaviour using linear or nonlinear mixed-effects h–d equations for other tree species (e.g. Jayaraman and Lappi 2001; Calama and Montero 2004; Trincado et al. 2007).

For the whole dataset (fitting plus evaluation), there was a modal value of six species present in each plot (SD = 1.75). This means that if one tree per species was used for calibration, an average of six trees would have to be measured in each plot to calibrate Eq. 11. This number is substantially lower than the number of trees that should be measured to calculate dominant height (a variable included in Eq. 12) (25 dominant trees per plot in this study, as all the plots had an area of 2,500 m²). One problem is whether we should focus our efforts on measuring only the height of dominant trees, or measuring heights for all species. Foresters usually prefer to measure dominant height, since it is one of the key driving variables in growth and yield models and can, sometimes, be obtained as an output of growth simulators. However, if the purpose of the h–d relationship is only to know missing heights, the calibrated version of Eq. 11 showed better results in terms of prediction accuracy (as showed by RMSE and E values) and measurement effort. In any case, as the area of the sample plot decreases, the number of trees used to calculate dominant height also decreases, thus reducing the measurement effort required for the application of the fixed-effects of Eq. 12.

When fitting Eqs. 11 and 12, some parameters were not significant at the 0.05 level (Table 2). If these parameters were assigned a value of 0, then some variables would be eliminated from the model, resulting in a less accurate model. We therefore decided to keep these variables in the model and assign a value of 1 to the parameters (note that for parameters that are multiplicative or in the exponent of a power term, assigning a value of 1 means that this parameter is dropped from the model). Therefore, slightly modified models were used in these cases. Convergence was achieved in all cases with the simpler models, and the accuracy of the models (showed by the RMSE and E values) was higher than that obtained by simply eliminating the parameters and their associated variables.

Calibration of mixed-effects models has been shown to be a good alternative for height prediction even when using a basic h–d model and only one height measurement per sampling level is available (e.g. Trincado et al. 2007), as in the current study. Other studies report similar results for height prediction (e.g. Calama and Montero 2004; Castedo-Dorado et al. 2006; Sharma and Parton 2007) as well as prediction of other growth variables (e.g. Calama and Montero 2005). Different calibration methods were tested for Eq. 11 in this study, because it was the equation that showed the best trade-off between fit statistics and measurement effort. The medium-size tree selection method, the random tree selection method and the all-size tree selection method performed better than others. This can be attributed to the fact that most forests in Durango are irregular, and the smaller and larger diameters show more variation among stands, whereas the mean diameters are much more similar between stands. Some authors (e.g. Castedo-Dorado et al. 2006; Crecente-Campo et al. 2010) obtained the best results when the smaller trees in the plot were used for calibration and attributed this to the fact that their models included dominant height as a fixed effect, and therefore heights corresponding to the larger trees did not provide much additional information for calibration. In the current study, Eq. 11 does not contain stand variables as regressors, but the large-tree method gave the worst results in calibration (Fig. 6). However, the smaller trees were not the best possible sample in this study, perhaps because the selected model was not restricted to pass through any point. In the aforementioned studies, the models were restricted to pass through the point (ddom, hdom), implying that the models cannot change much in the upper part of the h–d relationship, and small trees provided much more information for calibration than larger trees (Crecente-Campo et al. 2010). The greater the number of measurements included in the subsample, the greater the decrease in RMSE and increase in EF (Fig. 5). However, a large sample is often not justifiable, because of the higher cost of sampling (Castedo-Dorado et al. 2006). This was particularly important in this study, since calibration should be carried out independently for each sampling unit.

For the species-specific Eqs. 11 and 12, calibration can be done for each plot independently, as the plot-specific random effect associated with each species is different. A macro program may be very useful in this case. However, when calibrating Eqs. 9 and 10, calibration should be done for the whole dataset for which we want to obtain calibrated predictions (see Appendix A). This is because each species-specific random effect is the same for all the plots in which one species is present, and each plot-specific random effect is the same for all the species included in a particular plot. Therefore, each plot–species combination has a different random effect, obtained by summing the corresponding plot-specific and species-specific random effects.

Height–diameter models have previously been developed for other multi-species, multi-height and multi-age forests (e.g. Sharma and Parton 2007; Vargas-Larreta et al. 2009). However, the developed models were always species-specific (as Eqs. 11 and 12 in this study), and the utility of developing a general model like Eqs. 9 and 10 had not been shown. The approach used in this study could be applied in other irregular stands, where a general model with species-specific and plot-specific random effects could be calibrated by measuring only one tree of each species per plot.

However, the specific models developed in this study can only be used for the forests in Durango, because the variance and covariance parameters were calculated with the variability observed in the research plots, which were all located in this state. Variability in other areas, such as the bordering states of Chihuahua, Nayarit and Jalisco, where the Sierra Madre Occidental is located, may be different and therefore any predictions in these areas using these models should be taken with caution until the developed models are tested.

Following these results, we therefore recommend measuring one tree height for each species in the stands in order to apply the species-specific Eq. 11. We recommend the measurement of medium-sized trees for this purpose, because the random tree selection method can give worse results than the ones showed in this study if only one repetition is carried out (note that the values showed for the statistics of the random tree selection method are mean values after 100 simulations). If heights are not measured, the uncalibrated Eq. 11 should be used. If dominant height is available, and assuming that trees per hectare and basal area are available because all diameters would have been measured, then the uncalibrated Eq. 12 should be used. Finally, Eq. 12 can also be calibrated if dominant height and additional height measurements are available.

5 Conclusions

Random components of a basic h–d mixed model performed very well in explaining the observed differences in the h–d relationship between different plots and species in mixed uneven-aged stands in Durango (Mexico). However, the variance of the random effects for the different ecoregions was non-significant. We likely do not have enough information to estimate such variance since only four ecoregions were included in the dataset. A generalised h–d model is not necessary if a calibrated basic h–d model is used, as similar results are obtained when only one tree per species is used in the calibration process. If the models are not being calibrated, species-specific models should be used because they at least capture the variation between species.

Appendix. Height estimation—an example

This Appendix shows an example of the height estimation using the calibrated Eqs. 9 and 11. Equations 10 and 12 should be calibrated in the same way as Eqs. 9 and 11, respectively.

We suppose that the following data belong to three new plots in which several species are present, but we are interested in calibrating the models for species 1, 2 and 3, which are the commercial ones. Thus, we measure one tree of these species in each one of the plots, resulting in:

Plot	Species	d (cm)	h (m)
1	1	15	12
1	2	17	12
1	3	14	12
2	1	20	13
2	2	19	13
2	3	16	13
3	1	17	11
3	2	15	11
3	3	18	11

1.1 Case 1: use of Eq. 9

In this case, calibration should be done for the whole dataset for which we want to calibrate Eq. 9, using the expression shown in Eq. 6.

Note that subscript i in Eq. 6 is finally not necessary because there were no differences between ecoregions, so hereafter it was omitted.

The estimated variances and the covariance of the random effects (Table 3) are the elements of the variance–covariance matrix $\widehat{\mathbf{D}}$, which should be arranged according to the dataset in the following way: three species, with two random parameters by species and three plots with three random parameters by plot, resulting in a 15 × 15 diagonal matrix:

$$ \widehat{\mathbf{D}}=\left[\begin{array}{ccccccccccccccc}\hfill 0.01090\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0.01090\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0.01090\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.07487\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.07487\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.07487\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.1206\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.1206\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.1206\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.06554\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.06554\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.06554\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.7793\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.7793\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.7793\hfill \end{array}\right] $$

The variance–covariance matrix for the random error term is determined by assuming that all estimations have constant variance (σ ², Table 3) and that the errors are not correlated:

$$ {\widehat{\mathbf{R}}}_{jk}={\sigma}^2\times {\mathbf{I}}_9=\left[\begin{array}{ccccccccc}\hfill 5.678\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 5.678\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 5.678\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 5.678\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 5.678\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 5.678\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 5.678\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 5.678\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 5.678\hfill \end{array}\right] $$

where I ₉ is the identity matrix with dimension (9 × 9) equal to the number of data used for calibration.

The partial derivatives matrix ${\widehat{\mathbf{Z}}}_{jk}$ with respect to the random effects is calculated with the partial derivatives of Eq. 9. Using $z=\left[1- \exp \left(-\frac{a_2}{50}\cdot {d}_{jkl}\right)\right]/\left[1- \exp \left(-\frac{a_2}{50}\cdot 100\right)\right]$ for simplification, the partial derivatives are:

$$ \begin{array}{l}\frac{\partial {h}_{jkl}}{\partial {w}_{0k}}={\left(\frac{1}{a_0}\right)}^2{\left[\left({\left(20{a}_1\right)}^{a_0}-{1.3}^{a_0}\right)z+{1.3}^{a_0}\right]}^{\frac{1}{a_0}}\cdot \\ {}\kern3em \left[\frac{a_0\left(\left( \ln 20+ \ln {a}_1\right){\left(20{a}_1\right)}^{a_0}z-{1.3}^{a_0}z \ln 1.3+{1.3}^{a_0} \ln 1.3\right)}{\left({\left(20{a}_1\right)}^{a_0}-{1.3}^{a_0}\right)z+{1.3}^{a_0}}- \ln \left(\left({\left(20{a}_1\right)}^{a_0}-{1.3}^{a_0}\right)z+{1.3}^{a_0}\right)\right]\end{array} $$

$$ \frac{\partial {h}_{jkl}}{\partial {w}_{1k}}={20}^{a_0}{a}_1{}^{a_0-1}z{\left[\left({\left(20{a}_1\right)}^{a_0}-{1.3}^{a_0}\right)z+{1.3}^{a_0}\right]}^{\left(\frac{1}{a_0}-1\right)} $$

$$ \frac{\partial {h}_{jkl}}{\partial {v}_{0j}}=\frac{\partial {h}_{jkl}}{\partial {w}_{0k}} $$

$$ \frac{\partial {h}_{jkl}}{\partial {v}_{1j}}=\frac{\partial {h}_{jkl}}{\partial {w}_{1k}} $$

$$ \begin{array}{l}\frac{\partial {h}_{jkl}}{\partial {v}_{2j}}=\frac{\left({1.3}^{a_0}-{\left(20{a}_1\right)}^{a_0}\right) \exp \left[\left(2-\frac{d_{jkl}}{50}\right){a}_2\right]}{\left(50{a}_0{\left( \exp \left(2{a}_2\right)-1\right)}^2\right)}.\\ {}\kern3em \left({d}_{jkl}\left(- \exp \left(2{a}_2\right)\right)+100\left( \exp \left(\frac{a_2}{50}{d}_{jkl}\right)-1\right)+{d}_{jkl}\right)\cdot {\left[\left({\left(20{a}_1\right)}^{a_0}-{1.3}^{a_0}\right)z+{1.3}^{a_0}\right]}^{\left(\frac{1}{a_0}-1\right)}\end{array} $$

The values for the partial derivatives are:

Plot	Species	d (cm)	h (m)	$\frac{\partial {h}_{jkl}}{\partial {w}_{0k}}$	$\frac{\partial {h}_{jkl}}{\partial {w}_{1k}}$	$\frac{\partial {h}_{jkl}}{\partial {v}_{0j}}$	$\frac{\partial {h}_{jkl}}{\partial {v}_{1j}}$	$\frac{\partial {h}_{jkl}}{\partial {v}_{2j}}$
1	1	15	12	6.0907	8.2505	6.0907	8.2505	3.3777
1	2	17	12	6.0354	9.0770	6.0354	9.0770	3.5535
1	3	14	12	6.0887	7.8158	6.0887	7.8158	3.2712
2	1	20	13	5.8419	10.2163	5.8419	10.2163	3.7346
2	2	19	13	5.9183	9.8494	5.9183	9.8494	3.6844
2	3	16	13	6.0720	8.6707	6.0720	8.6707	3.4716
3	1	17	11	6.0354	9.0770	6.0354	9.0770	3.5535
3	2	15	11	6.0907	8.2505	6.0907	8.2505	3.3777
3	3	18	11	5.9835	9.4697	5.9835	9.4697	3.6243

The partial derivatives should be arranged in the matrix into different columns, to reflect how the data is arranged in the dataset, to obtain a 9 × 15 matrix (9 observations used in calibration and 15 random parameters) as follows:

$$ {\widehat{\mathbf{Z}}}_{jk}=\left[\begin{array}{ccccccccccccccc}\hfill \frac{\partial {h}_{jk l}}{\partial {w}_{0k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{1k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{0j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{1j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{2j}}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{0k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{1k}}\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{0j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{1j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{2j}}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{0k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{1k}}\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{0j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{1j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{2j}}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill \frac{\partial {h}_{jk l}}{\partial {w}_{0k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{1k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{0j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{1j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{2j}}\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{0k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{1k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{0j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{1j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{2j}}\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{0k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{1k}}\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{0j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{1j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{2j}}\hfill & \hfill 0\hfill \\ {}\hfill \frac{\partial {h}_{jk l}}{\partial {w}_{0k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{1k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{0j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{1j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{2j}}\hfill \\ {}\hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{0k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{1k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{0j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{1j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{2j}}\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{0k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {w}_{1k}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{0j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{1j}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial {h}_{jk l}}{\partial {v}_{2j}}\hfill \end{array}\right] $$

resulting in:

$$ {\widehat{\mathbf{Z}}}_{jk}=\left[\begin{array}{ccccccccccccccc}\hfill 6.0907\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 8.2505\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 6.0907\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 8.2505\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 3.3777\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 6.0354\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 9.0770\hfill & \hfill 0\hfill & \hfill 6.0354\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 9.0770\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 3.5535\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 6.0887\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 7.8158\hfill & \hfill 6.0887\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 7.8158\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 3.2712\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 5.8419\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 10.2163\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 5.8419\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 10.2163\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 3.7346\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 5.9183\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 9.8494\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 5.9183\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 9.9494\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 3.6844\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 6.0720\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 8.6707\hfill & \hfill 0\hfill & \hfill 6.0720\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 8.6707\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 3.4716\hfill & \hfill 0\hfill \\ {}\hfill 6.0354\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 9.0770\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 6.0354\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 9.0770\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 3.5535\hfill \\ {}\hfill 0\hfill & \hfill 6.0907\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 8.2505\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 6.0907\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 8.2505\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 3.3777\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 5.9835\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 9.4697\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 5.9835\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 9.4697\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 3.6243\hfill \end{array}\right] $$

The matrix ${\widehat{\mathbf{e}}}_{jkl}$ containing the errors obtained with the fixed-effects model (Eq. 9 with w _0k , w _1k , v _0j , v _1j and v _2j equal 0) is:

$$ {\widehat{\mathbf{e}}}_{jkl}=\left[\begin{array}{c}\hfill 1.7949\hfill \\ {}\hfill 0.9066\hfill \\ {}\hfill 2.2621\hfill \\ {}\hfill 0.6825\hfill \\ {}\hfill 1.0767\hfill \\ {}\hfill 2.3432\hfill \\ {}\hfill -0.09339\hfill \\ {}\hfill 0.7949\hfill \\ {}\hfill -0.5154\hfill \end{array}\right] $$

Therefore, the random parameters calculated using Eq. 6 were:

$$ {\widehat{\mathbf{b}}}_{jk}=\left[\begin{array}{c}\hfill {w}_{01}\hfill \\ {}\hfill {w}_{02}\hfill \\ {}\hfill {w}_{03}\hfill \\ {}\hfill {w}_{11}\hfill \\ {}\hfill {w}_{12}\hfill \\ {}\hfill {w}_{13}\hfill \\ {}\hfill {v}_{01}\hfill \\ {}\hfill {v}_{02}\hfill \\ {}\hfill {v}_{03}\hfill \\ {}\hfill {v}_{11}\hfill \\ {}\hfill {v}_{12}\hfill \\ {}\hfill {v}_{13}\hfill \\ {}\hfill {v}_{21}\hfill \\ {}\hfill {v}_{22}\hfill \\ {}\hfill {v}_{23}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill 0.01561\hfill \\ {}\hfill 0.02045\hfill \\ {}\hfill 0.05229\hfill \\ {}\hfill 0.002747\hfill \\ {}\hfill 0.009959\hfill \\ {}\hfill 0.03455\hfill \\ {}\hfill 0.05549\hfill \\ {}\hfill 0.04252\hfill \\ {}\hfill -0.0002459\hfill \\ {}\hfill 0.03248\hfill \\ {}\hfill 0.02262\hfill \\ {}\hfill -0.01373\hfill \\ {}\hfill 0.1756\hfill \\ {}\hfill 0.1310\hfill \\ {}\hfill -0.03699\hfill \end{array}\right] $$

This estimation of the random effects can be used to calibrate the model to obtain a specific response for the plots and species of this example, using Eq. 9.

The estimated heights and associated errors using the fixed effects Eq. 9 and the calibrated Eq. 9 are, therefore:

Plot	Species	d (cm)	h (m)	h jkl fixed	e jkl fixed	h jkl cal.	e jkl cal.
1	1	15	12	10.21	1.79	11.54	0.46
1	2	17	12	11.09	0.91	12.58	−0.58
1	3	14	12	9.74	2.26	11.55	0.45
2	1	20	13	12.32	0.68	13.41	−0.41
2	2	19	13	11.92	1.08	13.11	−0.11
2	3	16	13	10.66	2.34	12.22	0.78
3	1	17	11	11.09	−0.09	10.95	0.05
3	2	15	11	10.21	0.79	10.17	0.83
3	3	18	11	11.52	−0.52	11.89	−0.89
				${\displaystyle \sum \left|{\widehat{\mathbf{e}}}_{jkl}\right|}$	10.47		4.54

It can be observed how the prediction errors are reduced from using a fixed-effects prediction to using a calibrated prediction.

1.2 Case 2: use of Eq. 11

This is the usual way of calibrating a mixed-effects model with only one sampling level (i.e. the plot level), as done in many previous studies (e.g. Calama and Montero 2004; Castedo-Dorado et al. 2006; Sharma and Parton 2007). In this case, the calibration process has to be done for each species independently, because each species has a different model (Eq. 11 with the parameters shown in Table 3).

Hereafter, we include an example for species 1 on plot 1. The other species and plots would be calibrated in the same way.

Not that for Eq. 11 and species 1, only the variance for the random effect v ₁ was significant, therefore only one random effect has to be calculated.

The estimated variances and the covariance of the random effects (Table 3) are the elements of the variance–covariance matrix $\widehat{\mathbf{D}}$, conforming a simple (1 × 1) diagonal matrix (1 random parameter):

$$ \widehat{\mathbf{D}}=\left[0.05110\right] $$

The variance–covariance matrix for the random error term is determined by assuming that all estimations have constant variance (σ ², Table 3) and that the errors are not correlated. In this case, it only contains one observation:

$$ {\widehat{\mathbf{R}}}_{jk}={\sigma}^2\times {\mathbf{I}}_1=\left[5.308\right] $$

where I ₁ is the identity matrix with dimension (1 × 1) equal to the number of data points used in calibration.

The partial derivatives matrix ${\widehat{\mathbf{Z}}}_{jk}$ with respect to the random effects is calculated with the partial derivatives of Eq. 11, which have the same expression as those from Eq. 9.

The values for the partial derivatives are:

Plot	Species	d (cm)	h (m)	$\frac{\partial {h}_{jkl}}{\partial {v}_{1j}}$
1	1	15	12	7.7013

Therefore the corresponding matrix is:

$$ {\widehat{\mathbf{Z}}}_{jk}=\left[\frac{\partial {h}_{jk l}}{\partial {v}_{1j}}\right]=\left[7.7013\right] $$

The matrix ${\widehat{\mathbf{e}}}_{jkl}$ containing the errors obtained with the fixed-effects model (Eq. 11 with v _0jk , v _1jk and v _2jk equal 0) is:

$$ {\widehat{\mathbf{e}}}_{jkl}=\left[1.2840\right] $$

Therefore, the random effect estimated using Eq. 6 was:

$$ {\widehat{\mathbf{b}}}_{jk}=\left[{v}_{11}\right]=\left[0.06060\right] $$

This estimation of the random effect can be used to calibrate the model to obtain a specific response for the plot and species of this example, using Eq. 11.

The estimated heights and associated errors using the fixed effects of Eq. 11 and the calibrated Eq. 11 are, therefore:

Plot	Species	d (cm)	h (m)	h jkl fixed	e jkl fixed	h jkl cal.	e jkl cal.
1	1	15	12	10.71	1.28	11.18	0.82