Nonlinear fixed and random generalized height–diameter models for Portuguese cork oak stands

Abstract

• Objective

The objective of the research was to develop a generalized height–diameter model for Quercus suber L. in Portugal, which can be applied both to undebarked and debarked trees, with diameter at breast height over cork larger than 2.5 cm.

• Methods

A nonlinear fixed effects model (NLFEM) and a nonlinear mixed effects model (NLMEM) approaches were used. Parameters estimates were obtained using the SAS macro NLINMIX, which uses a linear approximation to the marginal likelihood function by expanding it with a first-order Taylor series on the random effects. The option of expanding on the random effects at their current empirical best linear unbiased predictors (EBLUP) was used. The fitted models were evaluated using an independent data set, together with an existing model specific for undebarked trees. To obtain subject specific predictions with the NLMEM, a conventional and an improved calibration procedures were applied, considering four different tree sub-sampling designs. Both proposed models included dominant height and stand density as covariates to explain plot variability.

• Conclusions

Validation indicated that, even in the situations where the NLMEM calibration is not possible, this model should be preferred. The differences between the validated models, which were more evident for young stands, were considered. No large differences in predictive accuracy were found between the calibrated NLMEM using the conventional or the improved calibration procedures, for all the considered sub-sampling designs.

1 Introduction

An accurate estimate of total tree height (h) is often required in forest management and in growth and yield studies. However, height is expensive to measure, and therefore tree height is often measured only for a sample of trees. Height of the trees not measured in this sample can be estimated using a height–diameter model. These models can also be used to estimate height growth from diameter growth estimates in forest growth and yield models.

The height–diameter relationship varies from stand to stand, and in the same stand the relationship varies over time. Therefore, a simple curve with the diameter as the sole regressor has a local and temporal applicability to the stand where the fit data was taken from. A generalized height–diameter model, with a wider range of geographical applicability, can be built by taking into account stand variables (e.g., Sharma and Parton, 2007; Temesgen and von Gadow, 2004), using models with random coefficients for each sampling unit (e.g., Lappi, 1997) or including fixed dummy variables for regional effects (e.g., Huang et al., 2000). All or some of these approaches can be included in the same model (e.g., Adame et al., 2008; Calama and Montero, 2004; Saunders and Wagner, 2008).

Cork oak (Q. suber) is one of the most valuable tree species in Portugal, generally oriented to the production of cork. The cork is an outer layer in the bark periderm built up by the phellogen. The removal of the cork layer is carried out at periodic intervals (every 9 years or more after the perimeter at breast height over bark achieves 70 cm), and due to the species capacity of regenerating a traumatic phellogen, the cork layer is sustainably renewed during the tree’s lifetime.

Few generalized and local height-diameter models for Q. suber can be found in the literature. Sánchez-González et al. (2007), using data from the second Spanish Forest Inventory, has built a generalized height–diameter model for cork oak trees in Spain. The Spanish Forest Inventory uses the value of 7.5 cm as a threshold for diameter measurement, and additionally the data set did not include trees larger than 75 cm in diameter under bark at breast height, neither trees in stands with dominant height values greater than 14 m that are common in Portugal (Cañellas et al., 2008). Although the data used to fit the Sánchez-González et al. (2007) height–diameter model only included trees in which cork thickness had been measured (three or four trees randomly selected in each plot), stand variables such as dominant or maximum diameter under cork for each plot were computed using the diameter at breast height under cork, estimated by subtracting the mean cork thickness of the three or four trees where cork thickness was measured, assuming the same cork age for all trees in each plot (Sánchez-González et al. 2007). Paulo and Tomé (2009) developed a height–diameter model applicable only to juvenile trees (trees before the first debarking) in even-aged stands, using diameter over cork as regressor and for the computation of the diameter-dependent stand level regressor variables.

In the present research a generalized height–diameter model was developed for cork oak in Portugal that can be applied to trees with diameter over cork larger than 2.5 cm, located in stands where mature (debarked at least once) and juvenile trees may coexist. This is a frequent situation since it is common for trees in the same stand and with the same age to have the first debarking occurring in different years.

Recognizing the lack of independence of the data that are collected on trees sampled in plots from multiple stands, the use of mixed effects model was explored. The inclusion of random parameters allows modeling the variability among different plots. The inclusion of random effects provides consistent estimates of the fixed parameters and their standard errors, and may improve the predictions obtained if it is possible to estimate the random effects for a plot that has not been sampled. This is accomplished by an approach known as calibration or localization, which demands supplementary observations of the dependent variable (total height in this case).

In order to determine the importance of the mixed model approach on the improvement of the height predictions, the resulting mixed model was compared for bias and precision with the fixed effects model, also developed during the research. The prediction performances of the developed models, when applied to juvenile stands, were compared with the prediction accuracy of the height–diameter model of Paulo and Tomé (2009), which has been specifically built for juvenile stands.

All different alternative models were compared using an independent data set. The subject-specific responses from the mixed effects models were estimated considering two distinct approaches proposed in the literature and four different tree sampling designs, in order to find out which design should be preferred when performing the mixed effect model calibration.

2 Material and methods

2.1 Data

The data set used to develop the model consists of 5,982 height–diameter pairs of measurements taken from 112 plots well-distributed over the cork oak distribution area in Portugal, with a wide range of ages of trees and stand characteristics.

Plots in non-debarked stands (juvenile stands), in a total of 33 plots, were established between 1993 and 2006, with plot area varying between 2,025 and 2,827 m². Thirteen of these plots were re-measured once at a 3 or 6 year interval.

Data from adult stands, in a total of 70 plots, came from 53 inventory plots distributed over four stands and measured once between 2000 and 2004, and 17 permanent plots distributed over the main regions of cork production and measured during 2001. The plot areas ranged from 2,827 m² to 5,184 m².

An independent data set was used for calibration and validation of the chosen height–diameter models. It includes data from a collection of 15 permanent plots installed during 2007 (not yet remeasured) in four different stands where no trees had yet been debarked, with stand ages going from 9 to 15 years. It also includes data from 18 inventory plots measured during 2004 and eight permanent plots measured in 2007, all located in mature stands with unknown age.

Diameters at breast height were measured for every tree in the two datasets, using calipers or diameter tapes for bigger trees. In mature trees, diameter under bark was directly measured after debarking or, when this measurement was not possible, calculated from diameter over bark and bark thickness evaluated as the average of two measurements with a probe. In juvenile trees, due to the difficulty of measuring virgin-cork thickness (due to its elasticity and its uneven surface with several cracks due to the quick initial diameter growth), cork thickness was estimated using the model proposed by Tomé (2004), which was developed with data taken from stem disks. Crown radii along the four cardinal directions were measured for every tree. Total height of the trees was measured with an ultrasonic hypsometer.

The following stand variables were computed for each plot: number of trees per ha (N), basal area under bark (Gu), minimum diameter under bark (dumin), maximum diameter under bark (dumax), dominant diameter under bark (dudom), dominant height (hdom), quadratic mean diameter under bark (dug), spacing coefficient (Sc), and several distance independent competition indices defined as the ratios between tree diameter under bark and hdom, dudom, dug, dumin or dumax. Dominant diameter under bark and dominant height were defined as the mean diameter and the mean height of the 25 thickest trees per hectare, and not of the 100 thickest trees per hectare, as is usual, due to the very common low density values of the cork oak stands. Spacing coefficient is defined as

$$ {\text{S}}{{\text{c}}_{\text{ij}}} = \frac{{100}}{{\sqrt {{{{\text{N}}_{\text{ij}}}}} \,\overline {{\text{c}}{{\text{w}}_{\text{ij}}}} }} $$

(1)

where N_ij is the number of trees per hectare and $\overline {{\text{c}}{{\text{w}}_{\text{ij}}}}$ is the mean crown diameter of the trees, in the i^th plot at the j^th measurement for both variables. This is a commonly used index to evaluate cork oak stand density (Natividade, 1950), as it represents the ratio between the average distance between trees if they were regularly spaced $\left( {\frac{{100}}{{\sqrt {{{{\text{N}}_{\text{ij}}}}} }}} \right)$ and the mean crown diameter of the trees $\overline {{\text{c}}{{\text{w}}_{\text{ij}}}}$.

It is important to notice that although the developed height–diameter model will be used to estimate height from trees with diameter over cork larger than 2.5 cm, all stand variables were computed considering only trees with diameter over cork larger than 7.5 cm, since this is the definition used in the Portuguese National Forest Inventory and in most forest inventories carried out for management purposes.

After a graphical analysis at the plot level, a small number of extreme observations, attributed to measurements errors, were excluded from the data sets. Trees that were dead, damaged or forked below breast height or presenting a broken top were also excluded. Summary of tree and stand-level descriptive attributes for both fit and validation data sets are presented in Table 1.

Table 1 Summary statistics for fitting and validation data sets

2.2 Model selection

Several functions have been used to model the height–diameter relationship (e.g., Fang and Bailey, 1998; Huang et al., 1992; Temesgen and von Gadow, 2004).

According to Lei and Parresol (2001), the function form used to model the height-diameter relationship should: (i) increase monotonically, (ii) have an upper asymptote, and (iii) have an inflection point. We questioned the third requirement. Both diameter and height growth curves have an inflection point, but this may not be necessarily so for the relationship between height and diameter. Furthermore, although the fitting data set includes trees with large variability of diameter values (corresponding to young and old individuals), no evidence of an inflection point was found when plotting height against diameter under cork with our fitting data (Fig. 1).

Fig. 1

Plot of total tree height (h) against diameter at breast height under cork (du) for model fitting data (n = 5,982)

The functions included in the comprehensive list presented by Huang et al. (2000) were considered as possible models to be assessed in this study. Functions without an asymptote (functions 1, 7 and 8 in Huang et al., 2000) or functions with more than three parameters were discarded due to frequent convergence problems. Functions that were particular cases of others were also discarded (e.g., functions 23 and 24 in Huang et al., 2000) and only the more general functions were considered. The Schnute function (Schnute, 1981), in its original form, directly restricts the function to pass through two distinct points. We chose the points (0, 1.30) and (dudom, hdom) as did Dorado et al. (2006). A total of eleven models were considered for model selection, since they presented the required characteristics.

The candidate functions were modified, if necessary, to guarantee that they pass through the point (d = 0, h = 1.30 m) to prevent negative height estimates for small trees and/or to guarantee a good estimation for small trees. Table 2 shows the original and the reparameterized form of the considered functions, as well as the expression of the corresponding asymptotes.

Table 2 Formulation of the local height-diameter functions considered

The candidate functions were fitted independently to the data from each of the measurement dates of each plot by nonlinear least squares (NLS) method, using the MODEL procedure of SAS (SAS Institute Inc. 2004). For a selection of the functions that better describe the relationship between height and diameter at each measurement, the mean square error of prediction (MSE_ij) was calculated for each model for plot i at the measurement date j as

$$ {\text{MS}}{{\text{E}}_{\text{ij}}} = \frac{{\sum\nolimits_{{{\text{k}} = 1}}^{{{{\text{n}}_{\text{ij}}}}} {{{\left( {{{{{\hat{e}}}}_{\text{ijk}}}} \right)}^2}} }}{{{{\text{n}}_{\text{ij}}} - 1}}, $$

(2)

where ê_ijk is the residual from the k^th tree in the i plot at the measurement date j and n_ij is the number of trees in plot i at the measurement date j. Since each plot has only one or two measurements, j takes values 1 or 2. For each function, the sum of MSE obtained for all plots and measurement occasions, denoted as MSE_ij, was computed and used as criterion to compare the candidate models and select the model with the best value of this statistic. Additionally, adjusted R-square (adj-R²) obtained by fitting each equation to the complete data set was also used as a criterion for model selection.

Both statistics were only considered for plots and measurement occasion with more than 30 observations (a total of 77) to assure meaningful estimates.

2.3 Fitting strategies

2.3.1 Nonlinear fixed effects model

The selected model was fitted with nonlinear least squares techniques (maximum likelihood estimation methods) using the SAS macro NLINMIX (SAS Institute Inc. 2004) and the complete fit data set, without specifying any random effects in the model.

The parameters of the base equation were expressed as linear functions of stand-level variables and competition indices defined previously in the data section. The selection of the stand-level variables and/or competition indices was done with the help of plots of the individual coefficient estimates of the base equation for each plot and measurement occasion versus the corresponding values of the stand-level variables and competition indices.

In order to avoid multicollinearity problems due to over-parameterization of the model, only variables largely correlated with estimated coefficients were included. After including these covariates in the model, the significance of their coefficients was assessed by asymptotic t-tests.

If the plot of the studentized residuals versus the fitted values showed that the within-individual error variance was heterogeneous, variance functions were included. Two frequently used and flexile variance functions for growth modeling were tested (Davidian and Giltinan, 1995): the power function model

$$ {\text{g}}\left( {{{{\mu }}_{\text{ij}}}{{,\alpha }}} \right) = {{\mu }}_{\text{ij}}^{{{\alpha }}} $$

(3)

and the exponential function model

$$ {\text{g}}\left( {{{{\mu }}_{\text{ij}}}{{,\alpha }}} \right) = \exp \left( {{{\alpha }}{{{\mu }}_{\text{ij}}}} \right) $$

(4)

where α is a fixed effect parameter to be estimated, and ${{{\mu }}_{\text{ij}}}$ is the mean function, which is defined by ${{{\mu }}_{\text{ij}}} = {\text{E}}\left[ {{{\text{y}}_{\text{ij}}}\left| {\mathbf{\beta }} \right.} \right]$, where β is the (p x 1) vector of fixed effects parameters common to the population. Thus, the expression for the intra-individual variance-covariance matrix takes the form

$$ {{\mathbf{R}}_{\text{ij}}} = {{{\sigma }}^2}{{\mathbf{G}}_{\text{i}}} $$

(5)

where σ² is the value of the error variance of the model, and G_i is a (n_ij x n_ij) diagonal matrix whose elements that describe the non-constant variance of the errors are given by one of the variance function models. During the iterative estimating process, σ² and μ_ij are substituted by their current estimates, ${{{\hat{\sigma }}}^2}$ and ${{{\hat{\mu }}}_{\text{ij}}}$. The choice between the two functions was made with the graphical observation of the studentized residuals and the value of the Akaike information criterion (AIC).

2.3.2 Nonlinear mixed effects models

Using the selected base model, and to account for the correlation of the data due to their nested structure (trees inside plots), specific subject models considering both fixed and random parameters (each measurement occasion for each plot) were fitted. The large number of trees in comparison to the number of re-measurements for the same individual (measured twice maximum) makes it reasonable to assume an absence of serial correlation (Dorado et al., 2006; Vanclay et al., 1995).

A general expression for a nonlinear mixed effects model can be written as (e.g., Lindstrom and Bates, 1990; Pinheiro and Bates, 2000; Vonesh and Chinchilli 1997)

$$ {{\mathbf{y}}_{{{\mathbf{ij}}}}} = {\text{f}}\left( {{{{\Phi }}_{\text{ij}}}\,,\,{{\mathbf{x}}_{\text{ij}}}} \right) + {{\mathbf{e}}_{{{\mathbf{ij}}}}} $$

(6)

where y_ij is the (n_ij x 1) vector with the n_ij observations of variable y from the i^th plot at the j^th measurement occasion, x _ij is the (n_ij x 1) vector from the n_ij observations from the predictor variable taken from the i^th plot at the j^th measurement occasion, f is a nonlinear function of the predictor variable and the parameter vector, and e _ij is a (n_ij x 1) vector from the error terms. The regression coefficients can be defined as

$$ {{\mathbf{\Phi }}_{{{\mathbf{ij}}}}} = {{\mathbf{A}}_{{{\mathbf{ij}}}}}{\mathbf{\beta }} + {{\mathbf{B}}_{{{\mathbf{ij}}}}}{{\mathbf{b}}_{{{\mathbf{ij}}}}} $$

(7)

where β is the (p x 1) vector of fixed effects parameters common to the population, and b_ij is the (q x 1) vector of random effects associated with the i^th plot at the j^th measurement occasion. A_ij and B_ij are design matrices (r x p) and (r x q), for the fixed and random effects specific to each plot at each measurement occasion, respectively, where r is the number of the parameters in the model, p is the number of fixed parameters, and q the number of random effects. Elements of the design matrices are 0, 1, or the values of covariates associated with the fixed or random effects. Basic assumptions for nonlinear mixed model fitting include the multivariate normal distribution for the random-effects and the error term vector:

$$ {{\mathbf{b}}_{{{\mathbf{ij}}}}}\sim {{\text{N}}_{\text{q}}}\left( {{\mathbf{0}},{\mathbf{D}}} \right) $$

(8)

$$ {{\mathbf{e}}_{{{\mathbf{ij}}}}}\sim {\text{N}}\left( {0,{{\mathbf{R}}_{{{\mathbf{ij}}}}}} \right) $$

(9)

where D is a positive-definite variance-covariance matrix (q x q), for the random effects assumed to be common to every subject, and R_ij is the (n_ij x n_ij) intra-individual variance–covariance matrix.

Parameter estimates and analysis were performed using SAS macro NLINMIX (SAS Institute Inc. 2004), which uses a linear approximation to the marginal likelihood function by expanding it with a first-order Taylor series on the random effects. The option of expanding on the random effects at their current empirical best linear unbiased predictors (EBLUP) was used.

First, the nature of the parameters of the selected equation was specified as fixed and random effects or purely fixed. Separate fits for each plot and measurement occasion were obtained, and the variability across subjects was assessed by considering the individual confidence intervals for the parameters, using only the plots and measurement occasion with more than 30 observations to assure meaningful estimates by separate fitting. The parameters with less overlap in confidence intervals across subjects were considered as needing the inclusion of random effects. With the complete fitting data set, likelihood ratio tests (LRT) and, for non-nested models, the AIC statistic values were used to compare models with alternative sets of fixed and random effects, and to confirm the choices suggested by the graphical screening. As two random effects were specified for the model, an unstructured variance–covariance matrix D for random effects (Schabenberger and Pierce, 2002) was first assumed.

If graphical diagnostics of studentized residuals showed that the within-individual error variance was heterogeneous, variance functions were included as described in the previous section.

At this stage, the need for an unstructured D matrix was tested, using a likelihood ratio test where the nested model had a diagonal matrix D.

As a next step, to explain variability of parameters among subjects, appropriate covariates were studied. To investigate possible relationships between mixed effects and individual subject attributes, estimates of the random effects were plotted against the potential covariates mentioned above for nonlinear fixed effects models. The fixed effects parameters of the equations were expressed as linear functions of the chosen stand-level variables and likelihood ratio tests (LRT) or, for non-nested models, the AIC statistics were used to compare models with alternative sets of covariates. Only variables strongly correlated with estimated random effects were included.

After the reformulation of the model with the inclusion of covariates, the needs for an unstructured variance–covariance matrix D and for the inclusion of a variance function were re-assessed.

2.4 Evaluation of the predictive performance of the models

2.4.1 Population-specific responses

Mean responses are obtained with the fitted nonlinear fixed effects model to predict the total heights of the trees in the validation data set, and with the nonlinear mixed effects models by setting the random effects as zero, as they appear in a linear form in the expression for the nonlinear mixed effects model [see equation (7) from section 2.3.2].

For each model and every k^th observation in this data set, the difference between observed and predicted height was computed (ê_k). Then, the mean residual (M_r) and the mean of the absolute value of the residuals (M_|r|), were calculated to evaluate bias and precision respectively, as follows:

$$ {{\text{M}}_{\text{r}}} = \frac{{\sum\nolimits_{{{\text{k = 1}}}}^{\text{n}} {{{{\widehat e }}_{\text{k}}}} }}{\text{n}} $$

(10)

$$ {{\text{M}}_{{\left| {\text{r}} \right|}}} = \frac{{\sum\nolimits_{{{\text{k}} = 1}}^{\text{n}} {\left| {{{{\widehat e }}_{\text{k}}}} \right|} }}{\text{n}} $$

(11)

where n is the number of observations in the validation data set. Modelling efficiency (ef) or the proportion of variation explained by the model was also computed as:

$$ {\text{ef}} = {1} - \frac{{\sum\nolimits_{{{\text{k}} = 1}}^{\text{n}} {{{\left( {{{{{\hat{e}}}}_{\text{k}}}} \right)}^2}} }}{{\sum\nolimits_{{{\text{k}} = 1}}^{\text{n}} {{{\left( {{{\text{y}}_{\text{k}}}{ - }\overline {\text{y}} } \right)}^{{2}}}} }}, $$

(12)

where y_k is height value for observation k, $\overline {\text{y}}$ is the mean value of the observed heights in the validating data set, and the other variables are as defined above.

The validation data set was also divided into ten total height and ten diameter classes according to their deciles, and the values of M_r and ${{\text{M}}_{{\left| {\text{r}} \right|}}}$ were computed for each class, to study how the predictive ability of the models varied with the dimensions of the trees.

When evaluating the precision of the predictions for the different height and diameter classes, the absolute value of the residuals was expressed as a percentage as follows:

$$ {{\text{M}}_{{{\text{p }}\left| {\text{r}} \right|}}}\% = {100}\,\frac{{\frac{{\sum\limits_{{1 = 1}}^{{^{\text{n}}{\text{p}}}} {\left| {{{{\widehat e }}_{\text{k}}}} \right|} }}{{{{\text{n}}_{\text{p}}}}}}}{{{{\text{c}}_{\text{p}}}}}, $$

(13)

where n_p is the number of observations in class p, and c_p is the mean total height of the observations in class p.

2.4.2 Subject-specific responses

It is possible to tailor the nonlinear mixed effects models to individual subjects (plots in our case), by measuring heights in a subsample of trees in the plot, together with the predictor variable and covariates of the model. In previous applications (e.g., Fang and Bailey, 2001; Calama and Montero, 2004) predictions of the random effects for nonlinear mixed effects models, that have been fitted by using the first-order Taylor series expansion at the EBLUP, took the following form (Vonesh and Chinchilli, 1997):

$$ {{\mathbf{\hat{b}}}_{{\mathbf{s}}}} \approx {\mathbf{\hat{D}\hat{Z}}}_{{\mathbf{s}}}^{{\mathbf{T}}}{\left( {{{{\mathbf{\hat{R}}}}_{{\mathbf{s}}}} + {{{\mathbf{\hat{Z}}}}_{{\mathbf{s}}}}{\mathbf{\hat{D}\hat{Z}}}_{{\mathbf{s}}}^{{\mathbf{T}}}} \right)^{{ - 1}}}\left( {{{\mathbf{y}}_{{\mathbf{s}}}} - {\text{f}}\left( {{{\mathbf{A}}_{{\mathbf{s}}}}{\mathbf{\hat{\beta }}},{{\mathbf{x}}_{{\mathbf{s}}}}} \right)} \right) $$

(14)

where ${\mathbf{\hat{D}}}$ is the (q x q) estimated variance-covariance matrix for the among plot variability (common for all plots), ${{\mathbf{\hat{R}}}_{{\mathbf{s}}}}$ is the estimated (n_s x n_s) variance–covariance matrix for the within-plot variability, y_s is the (n_s x 1) response variable vector, x_s is the (n_s x 1) predictor variable vector of dimension n_s, taken from plot s, ${\mathbf{\hat{\beta }}}$ is the (p x 1) estimated vector of β, and A_s is the (r x p) design matrix for the fixed effects. y_sand x_s are both measured in the subsample taken in plot s. The ${{\mathbf{\hat{Z}}}_{{\mathbf{s}}}}$ matrix of partial derivatives of dimension (n_s x q) is defined as

$$ {{\mathbf{\hat{Z}}}_{{\mathbf{s}}}} = {\left. {\frac{{\partial {\text{f}}\left( {{{\mathbf{A}}_{{\mathbf{s}}}}{\mathbf{\beta }} + {{\mathbf{B}}_{{\mathbf{s}}}}{{\mathbf{b}}_{{\mathbf{s}}}},{{\mathbf{x}}_{{\mathbf{s}}}}} \right)}}{{\partial {\text{b}}_{\text{s}}^{\text{T}}}}} \right|_{{{{\beta = \hat{\beta },}}{{\text{b}}_{\text{s}}}{ = 0}}}} $$

(15)

where B_s is the (r x q) design matrix for the random effects.

Once the random effects are predicted, predictions of the response variable can be obtained by

$$ {{\mathbf{\hat{y}}}_{{\mathbf{s}}}} = {\text{f}}\left( {{{\mathbf{A}}_{{\mathbf{s}}}}{\mathbf{\hat{\beta }}} + {{\mathbf{B}}_{{\mathbf{s}}}}{{{\mathbf{\hat{b}}}}_{{\mathbf{s}}}},{{\mathbf{x}}_{{\mathbf{s}}}}} \right). $$

(16)

The values obtained by this calibration procedure are referred as conventional predicted values.

Meng and Huang (2009) proposed the following equation for prediction of the random effects for nonlinear mixed effects models expanded at EBLUP (Lindstrom and Bates, 1990; Vonesh and Chinchilli, 1997)

$$ {\mathbf{\hat{b}}}_{{\mathbf{s}}}^{{\mathbf{k}}} \approx {\mathbf{\hat{D}}}{\left( {{\mathbf{\hat{Z}}}_{{\mathbf{s}}}^{{\left( {{\text{k}} - 1} \right)}}} \right)^{\text{T}}}{\left[ {{{{\mathbf{\hat{R}}}}_{{\mathbf{s}}}} + {{\mathbf{Z}}_{{\mathbf{s}}}}{\mathbf{\hat{D}}}{{\left( {{\mathbf{\hat{Z}}}_{{\mathbf{s}}}^{{\left( {k - 1} \right)}}} \right)}^T}} \right]^{{ - 1}}}\left[ {{{\mathbf{y}}_{{\mathbf{s}}}} - {\text{f}}\left( {{{\mathbf{A}}_{{\mathbf{s}}}}{\mathbf{\hat{\beta }}} + {{\mathbf{B}}_{{\mathbf{s}}}}{\mathbf{\hat{b}}}_{{\mathbf{s}}}^{{\left( {{\mathbf{k}} - {\mathbf{1}}} \right)}},{{\mathbf{x}}_{{\mathbf{s}}}}} \right)} \right] + {\mathbf{\hat{Z}}}_s^{{\left( {{\text{k}} - 1} \right)}}{\mathbf{\hat{b}}}_{{\mathbf{s}}}^{{\left( {{\mathbf{k}} - {\mathbf{1}}} \right)}} $$

(17)

where ${\mathbf{\hat{b}}}_{{\mathbf{s}}}^{{\left( {{\mathbf{k}} - {\mathbf{1}}} \right)}}$ and ${\mathbf{\hat{b}}}_{\text{s}}^{\text{k}}$ are the successive predictions of the random effects obtained iteratively and

$$ \widehat{{\mathbf{Z}}}^{{{\left( {{\text{k - 1}}} \right)}}}_{{\text{s}}} = \left. {\frac{{\partial {\text{f}}{\left( {{\mathbf{A}}_{{\mathbf{s}}} \beta + {\mathbf{B}}_{{\mathbf{s}}} {\mathbf{b}}_{{\mathbf{s}}} ,{\mathbf{x}}_{{\mathbf{s}}} } \right)}}} {{\partial {\text{b}}^{{\text{T}}}_{{\text{s}}} }}} \right|_{{\beta {\text{ = }}\widehat{\beta }{\text{,b}}_{{\text{s}}} {\text{ = }}\widehat{{\text{b}}}^{{{\left( {{\text{k - 1}}} \right)}}}_{{\text{s}}} }} . $$

(18)

The initial predictions of the random effects for the first iteration are the ones obtained by a conventional calibration procedure. The iteration procedure stops when the desired precision of the predictors is obtained, and the predictions of the response variable are then obtained, as mentioned above, for the conventional method. Meng and Huang (2009) found that when the improved procedure was applied to a nonlinear mixed effects tree height growth model, the bias of the local predictions and the variance of their errors were reduced.

The values obtained by this calibration procedure are referred as improved predicted values.

In order to compare both methods of prediction of the random effects, once the random effects for each plot in the validation data set were predicted using both the conventional and the improved methods, the heights of the other trees in each plot, which were not included in the subsample, were predicted. These two different calibration procedures were evaluated and compared using the values of the above-mentioned statistics (M_r, and M_|r|).

Using these same statistics, the proposed model was also compared with the model of Paulo and Tomé (2009), considering only the plots where no tree has yet been debarked, since this model was fitted to data from juvenile cork oak stands. This comparison is important in order to test if the models developed during this research, based on diameter under bark estimates for juvenile trees (trees with virgin cork), create predictions for juvenile stands comparable to those of a model that was developed specifically for these stands and based in diameter over cork measurements.

For both calibration methods, subsamples of heights from each plot in the validation data set were selected with the following designs:

• A. Total height of the three trees with the smallest diameter in each plot

• B. Total height of the dominant trees in each plot

• C. Total height of the first tree in each group of five trees consecutively measured in each 5 cm diameter class of the plot (1st, 6th, 11th …)

• D. Total height of the first tree in each group of five trees consecutively measured in the plot (1st, 6th, 11th …)

Designs A and B are non-random selection methods, and designs C and D represent two alternatives for tree systematic selection during forest inventory. If the trees are randomly distributed in the plots, the systematic selection is close to a random sample. Design A was considered since it has shown good results in other researches (e.g., Dorado et al., 2006). It is important to consider design B, since the height of dominant trees is usually measured for the estimation of dominant height and site index (when age is known), and therefore there are no additional costs to measure additional trees during the forest inventory. Design C corresponds to a systematic selection of trees based in the Draudt method (Patrone, 1963), very common in forest inventories in Portugal. Finally, design D represents one possible method for systematic selection of the trees that can be easily implemented in forest inventories. Due to the large variability found in cork oak stands, design D was thought to give better results than measuring, for instance, the 20% largest trees, as Calama and Montero (2004) did for stone pine plots, although both alternatives result in measuring a similar number of trees per plot.

3 Results

3.1 Model selection

When fitting the models in Table 2 to the fitting data set independently for each of the measurement dates of each plot, convergence for all of the plots was only achieved for four of the considered models, even when different initial values for the parameters were provided. Only these models were considered as candidate models for the next selection.

Table 3 shows the computed values of the sum of the mean square error (MSE_ij) and of the adjusted R-square (when the model was fitted to the complete data set) for each of the four height–diameter candidate models. Functions (9) and (11) have similar results. Function (11) was selected to model the height–diameter relationship, since it presented a slightly better value of adjusted R-square (adj-R²).

Table 3 Goodness of fit for the candidate local height-diameter functions used

3.2 Nonlinear fixed effects model

By plotting the individual coefficient estimates of the chosen base equation for each plot and measurement occasion against their stand-level variables and competition indices values, the covariates number of trees per hectare (N) and dominant height (hdom) were selected and the following final fixed effects model was obtained:

$$ {{\text{h}}_{\text{ijk}}} = 1.3 + \frac{{{\text{d}}{{\text{u}}_{\text{ijk}}}}}{{\left( {{\alpha_0} + \frac{{{\alpha_1}}}{{{{\text{N}}_{\text{ij}}}}}} \right) + \left( {{{{\beta }}_{{0}}}{ + }{{{\beta }}_{{1}}}\,{\text{hdo}}{{\text{m}}_{\text{ij}}}} \right)\,{\text{d}}{{\text{u}}_{\text{ijk}}}}} + {{\text{e}}_{\text{ijk}}}, $$

(19)

where variables are defined as Table 1.

None of the parameters seems to be correlated with any of the competition indices tested.

The graphical diagnostic of the studentized residuals versus the fitted values of the model showed heterogeneous error variance. For this reason, the two considered variance functions were evaluated. The computed AIC values were used to compare the two alternatives. The model including the exponential function resulted in a smaller AIC value (18,047) than the one with the power function (18,347), and the first was chosen.

Figure 2 shows the plots of the studentized residuals versus the fitted values resulting from the fitted model, before and after the inclusion of the variance function in the model. The Q-Q plot of the studentized residuals (not shown) showed that the residuals follow approximately the normal distribution. Table 4 presents the estimates of the fixed parameters of the final model.

Fig. 2

Plots of studentized residuals against predicted height values for the nonlinear fixed effects model, before (a) and after (b) the inclusion of the exponential variance function in the model

Table 4 Parameter estimates and variance components for the selected height-diameter models

3.3 Nonlinear mixed effects model

Plots of the ninety-five percent confidence intervals of the two coefficients of the chosen fixed base model (function 11 from Table 2) for each subject (plot in this case) gave a clear indication that random effects were needed to account for subject-to-subject variability in both coefficients.

By plotting the estimates of the random effects against the available covariates, a relationship between u₁ with the number of trees per hectare (N) and u₂ with dominant height (hdom) was evident (Fig. 3). Similarly to the results obtained for the nonlinear fixed effects model, no relationship was found between the random effects and the tree competition indices.

Fig. 3

Relationship between the two random effects and number of tree per hectare (N) and stand dominant height (hdom), before the introduction of covariates in the model parameters

Likelihood ratio tests and the AIC criterion were used to check if all the random effects and the stand level covariates were warranted. All strongly favor the more general model.

After the inclusion of these two variables in the model, the plot of the new estimates of the random effects against the other covariates showed no relationship, and the final model was defined according to the following expression, where variables are defined as in section 3.2.1 above:

$$ {{\text{h}}_{\text{ijk}}} = 1.3 + \frac{{{\text{d}}{{\text{u}}_{\text{ijk}}}}}{{\left( {{\alpha_{{0}}}{ + }\frac{{{\alpha_{{1}}}}}{{{{\text{N}}_{\text{ij}}}}}{ + }{{\text{u}}_{{1}}}} \right) + \left( {{{{\beta }}_{{0}}}{ + }{{{\beta }}_{{1}}}\,{\text{hdo}}{{\text{m}}_{\text{ij}}}{ + }{{\text{u}}_{{2}}}} \right)\,{\text{d}}{{\text{u}}_{\text{ijk}}}}} + {{\text{e}}_{\text{ijk}}} $$

(20)

As in the nonlinear fixed effects model, and despite the inclusion of random effects in the model, the graphical diagnostic of the studentized residuals showed heterogeneous error variance. For this reason, the two considered variance functions were evaluated and the computed AIC values used to compare the two alternatives. The model including the exponential function resulted in a smaller AIC value (17,061) than the one with the power function (17,237). Again, the exponential function was chosen.

The model was then re-assessed for the inclusion of all the random effects and the need for an unstructured variance–covariance matrix D, both of which proved to be significant in the final model.

Figure 4 shows the plot of the studentized residuals versus the fitted values resulting from the final model, before and after the inclusion of the variance function in the model. The Q–Q plot of the studentized residuals (not shown) showed that the residuals follow the normal distribution, supporting the normality of the model error as required by the basic assumptions for nonlinear mixed model fitting. The estimates of the fixed parameters for both models and variance parameters for the final mixed effects model are presented in Table 4.

Fig. 4

Plots of studentized residuals against predicted height values for the nonlinear mixed effects model, before (a) and after (b) the inclusion of the exponential variance function in the model

3.4 Evaluation of the predictive performance of the models

Table 5 presents the values of the mean residual, the mean of the absolute value of the residuals in percentage, and the model efficiency for the developed fixed and mixed effects models, computed from the complete validation data set. For the mixed effects model, the values are shown for the population-specific and for the subject-specific responses.

Table 5 Comparison of the models' predictive performances. Statistics were computed considering the complete validation data set. Tree sub-sampling designs are defined in section 2.4.2

The nonlinear mixed effect model (population- or subject-specific responses) is always more precise than the fixed effects model. Larger differences were found when the fixed effects model was compared with the subject-specific response obtained by the improved calibration procedure and tree sub-sampling designs B or D. In this case, the mean of the absolute value of residuals decreased from 0.8287 m to 0.6594 m, corresponding to 13.4% and 10.8% respectively of the mean value of the tree total heights (6.2 m) in the validation data set (Table 1), an improvement of 20.4% in the model precision.

However, the nonlinear fixed effects model predictions were less biased than those of the nonlinear mixed effects model, even when subject-specific responses are considered. The population-specific response resulted in a mean value of the residuals very similar to the one obtained when tree sub-sampling design A was used to calibrate the mixed model, indicating that measuring the height of the three trees with the smallest diameter in each plot had few impact on the bias of the model predictions.

To compare the models developed during this research with the Paulo and Tomé (2009) model for juvenile trees, in order to appraise if the estimation of diameter under cork for undebarked trees (where virgin cork is not usually measured) influenced the prediction capability of the developed models, the statistics presented in Table 5 were again computed considering only the 15 plots from the validation data set where no debarked trees exist, and consequently where the Paulo and Tomé (2009) model could be used for tree height estimation.

The Paulo and Tomé (2009) model, although developed specifically for cork oak juvenile stands, did not result in better estimates for juvenile stands (Table 6). This model only showed smaller values of the mean of the residuals and absolute residuals when compared with the fixed effects model. If the nonlinear mixed effects model is not calibrated, the use of the population-specific response would be recommended, since the results are similar to the Paulo and Tomé (2009) model and have no constraint in its application. This result supports the generalized use of the fitted mixed effects model in Portugal, avoiding the use of different models for juvenile and adult stands.

Table 6 Comparison of the models predictive performances. Statistics were computed considering only juvenile plots from the validation data set (plots with no debarked trees)

For the entire tree sub-sampling designs considered, the differences between the predicted tree heights using the conventional or the improved calibration procedures (Table 5 and 6) were small. In both cases, the sub-sampling designs B and D perform better.

The mean value of residuals and the mean of absolute value of residuals by height deciles classes, for the predictions of the fixed effects model, population-specific estimates of the mixed effect model, and subject-specific estimates for the improved calibration procedure considering tree sub-sampling design D are presented in Fig. 5. Differences in the predictive ability of the models are greater for trees less than 5 m in total height, mainly in terms of precision but also in terms of bias. Subject-specific estimation with the nonlinear mixed effects model has precision values less than 15% of the total tree height, a value that is only attained by the other two models (fixed effects model and population-specific model) for values of total height equal or greater to 6 m.

Fig. 5

Bias (mean value of the residuals) and precision (mean of absolute value of the residuals in percentage) for the considered models, by tree total height deciles classes. The horizontal axis values refer to the mean height of the observations in each total height decile class. Smoothed curves were superimposed on the data points

The limits of the 95% confidence intervals for the mean of the residuals in each decile class are presented in Table 7.

Table 7 Confidence intervals (α = 0.05) for the mean value of the residuals of each model, by each decile class

Analyzing the predictive ability by diameter classes (Fig. 6) again showed a similar trend, where the calibrated response is less biased and more precise in smaller trees with diameters less than 10 cm. For larger trees, the predictive ability of the population-specific estimates and subject-specific estimates of the mixed effect model present similar results.

Fig. 6

Bias (mean value of the residuals) and precision (mean of absolute value of the residuals in percentage) for the considered models, by diameter at breast height deciles classes (under cork). The horizontal axis values refer to the mean diameter of the observations in each diameter decile class. Smoothed curves were superimposed on the data points

The fixed effect, population-specific and subject-specific responses for four plots are presented, as an example, in Fig. 7. The large range of diameters under cork values among these four plots illustrates how the calibration of the nonlinear mixed effect model increases the predictive ability of the model, especially in plots with small trees (for example plot 3001001).

Fig. 7

Height–diameter prediction of the fixed effect, population-specific and subject-specific mixed effects models in four plots from the validation data set. Values of the horizontal axis refer to the average diameter at breast height from each of the decile classes

4 Discussion

One of the desirable characteristics of a height–diameter model is to pass through the point (d = 0, h = 1.30) to prevent negative or otherwise poor height estimates for small trees. This was considered important for the developed models, since they are expected to predict total height of trees as soon as the diameter at breast height reaches 2.5 cm. Despite the fact that many of the frequently used functions, in their original form, do not present this characteristic, they can be restricted to guarantee that they pass through this point. Many times this is not equal to simply adding a constant of 1.3 to the initial model like some authors did (e.g., Adame et al., 2008; Huang et al., 2000). All the tested candidate models were formulated by restricting them to passing through the point (d = 0, h = 1.30). Some of the functions presented in Huang et al. (2000) were discarded because they did not present an upper asymptote.

Sánchez-González et al. (2007) developed a generalized height–diameter model for Q. suber, but the selected function has no upper asymptote. Although the fitting data set did not include trees larger than 75 cm, from a biological point of view this is an important characteristic that a height–diameter model should present.

The function selected to model the height–diameter relationship for Q. suber, one of the functions proposed by Prodan (1968), is presented in Table 2 as function (11). This function includes two parameters (α related to the growth rate and β related to the upper asymptote), is increasingly monotonous, has an upper asymptote, and has a concave shape. Plotting height against diameter under cork on our fit data set did not show any evidence of an inflection point. Therefore, the function was thought to be appropriate for modelling the height–diameter relationship for Q. suber. A similar result was presented for Quercus pyrenaica Willd. by Adame et al. (2008).

The fitting characteristics of the fixed and mixed effects developed models improved when the parameters were expressed as a function of stand variables (stand density and dominant height). The relation of these two variables to height growth has already been referred by other authors (e.g., Vanclay, 2009; Saunders and Wagner, 2008; Calama and Montero, 2004), and the sign of the associated coefficients confirms that they are positively related with tree height. Generally, higher values of tree height for a given diameter value are expected to be found in stands with higher dominant height (associated to site index differences) and in denser stands (associated with the height/diameter ratio). The considered competition indices were never significant when simultaneously tested in the models with the included stand variables. This indicates that tree height is more dependent on stand characteristics than on tree position in the stand, since the fitting data set included plots with a large range of stand variable values.

Using an independent data set for validation, the predictive ability of the fixed effects, population-specific mixed effects and subject-specific mixed effects models was evaluated. A conventional and an improved methodology to obtain predictions of random effects for the nonlinear mixed effects model, expanded at EBLUP, were used.

Meng and Huang (2009) showed that the results obtained from these two methods are significantly different, and that the improved calibration procedure reduced bias and variance of the estimates. Using the computing program developed by Meng and Huang (2009) the estimates obtained by both procedures were compared, but the differences were not evident, despite the tree sub-sampling designs considered.

Possibly, none of the tested designs allowed an accurate prediction of the random effects due to great tree variability in the plots. This may result from the fact that Q. suber is among the Quercus species with the highest inter-tree variability, due to the high genetic diversity (Toumi and Lumaret, 1998). Furthermore, cork oak trees are traditionally and regularly pruned and the branch patterns and branch dominance in the trees are altered, possibly inducing a greater variability in the height–diameter relationship, except in the younger trees where pruning aims to remove the lower branches and to avoid the dominance of lateral branches. This could be an explanation for the greater impact of the calibration on the predictive ability of the nonlinear mixed effects model for smaller trees.

Tree sub-sampling designs D (measuring total height of 20% of the trees in the plot) and B (measuring height of the dominant trees in the plot) were found to be the designs that resulted in less biased predictions, although did not differ from the other two sub-sampling designs (A and C) for precision. Since dominant height is a regressor in the model that have to be estimated, the calibration procedure does not imply additional costs if design B is applied.

As dominant height is a regressor of the chosen mixed effects model, it was expected that although the sample is biased, with tree sub-sampling design A (measuring the total height of the three trees with the smaller diameter in each plot) the predictive ability of the calibrated model would have better results, since heights corresponding to the smallest trees would carry more additional information to the calibration procedure, as verified by Dorado et al. (2006).

However, the greater differences in the predictive ability of the models, as evaluated by the validation process, are found between the nonlinear fixed effects model and the calibrated responses of the mixed effects models. In some applications of height–diameter functions, such as their use as a module of an empirical growth model, calibration is not possible. For that situation, in accordance with the results obtained by other authors (e.g., Dorado et al., 2006; Temesgen et al., 2008), it was found that the population-specific response of the nonlinear mixed effects model is more accurate then the predictions of the nonlinear fixed effects model, and should be chosen.

When looking for the prediction ability along the total height classes, it was evident that the differences were larger for trees with total height under 5–6 m. This reinforced the usefulness of the mixed model approach, since in addition to improving the resulting estimates, the proposed mixed effects model can be generally used, since diameter at breast height reaches a value of 2.5 cm (over cork). Therefore, there is no need to use a specific model for juvenile trees. Furthermore, it was also shown that estimating diameter under cork for juvenile trees, using the Tomé (2004) model, does not decrease the prediction ability of the proposed model.