Generalized height-diameter and crown diameter prediction models for cork oak forests in Spain

A generalized height-diameter equation, along with a crown diameter prediction equation for cork oak forests in Spain were developed based on data from the Second Spanish Forest Inventory. Nine generalised height-diameter equations were selected as candidate functions to model the height-diameter under cork relationship, while for the crown diameter prediction model five linear and non-linear equations were tested. The equations were fitted using the mixedeffects model approach. The Stoffels & Van Soest power equation, constrained to pass through the point of dominant diameter and dominant height, was selected as the generalised height-diameter model. Regarding the crown diameter prediction model, the parable function without the intercept and with quadratic mean diameter incorporated as a fixed effect into the b parameter, proved to be the model with best prediction capabilities. The models were validated by characterising the model error using the PRESS (Prediction Sum of Squares) statistic. Both equations will be submodels of the ALCORNOQUE v1.0, a management oriented growth and yield model for cork oak forests in Spain.


Introduction
Tree height and crown dimensions are important tree characteristics used in many growth and yield models (Soares and Tomé, 2001).Height-diameter curves for tree species have been long used in forest inventories and growth models for predicting missing total height measurements (Curtis, 1967;Wykoff et al., 1982;Huang et al., 1992).In forest inventories, height is usually measured only for a subsample of trees, while diameter is measured for all the sampled trees.Height-diameter equations can either be used for local application or they can have a more generalised use (Krumland and Wensel, 1988;Tomé, 1989;Soares and Tomé, 2002).The former (local application) is normally only dependent on tree diameter and is only applicable to the stand where the height-diameter data were gathered, whereas generalized height-diameter equations are a function of tree diameter and stand variables and can be applied at the regional level.Height-diameter models are principally applied in height estimations in forest inventories
and as one of the main modules in management-oriented growth models.Crown width is used in tree and crown level growthmodelling systems, where simple competition indices are not available to adequately predict recovery from competition when a competitor is removed (Vanclay, 1994) and as predictor variable in diameter and height growth equations (e.g.Monserud and Sterba, 1996;Wykoff et al., 1982).Crown width is also used in calculating competition indices based on crown overlap (Biging and Dobbertin, 1992;Daniels et al., 1986).Crown width models can be formulated from opengrown trees or from stand-grown trees.Equations for predicting the dimensions of crowns in open locations consider maximum biological potential, and so are known as «maximum crown width» (MCW) equations, while those for stand-grown trees which generally have a smaller crown due to competition, are called «largest crown width» (LCW) equations (Hann, 1997).MCW models predict potential crown size and are primarily used in computing the crown competition factor (Krajiceck et al., 1961).LCW models predict the actual size of tree crowns in forest stands, and have many applications including estimations of crown surface area and volume in order to asses forest health (Zarnoch et al., 2004), tree-crown profiles and canopy architecture (Hann, 1999;Marshall et al., 2003), forest canopy cover (Gill et al., 2000) and the arrangement of trees in forest visualization programs (Hanus and Hann, 1998).
When modelling crown diameter, a simple linear model between crown width and diameter at breast height is often adequate (e.g.Cañadas, 2000;Paulo et al., 2002;Benítez et al., 2003).Other studies suggest that these linear models can be improved with quadratic expressions of diameter (Bechtold, 2003).Non-linear models have also been used, such as the power function and the monomolecular function (Bragg, 2001;Tomé et al., 2001).
The height-diameter and crown diameter prediction equations developed for other tree species used diameter over bark as a predictor variable (Curtis, 1967;Wykoff et al., 1982;Cañadas, 2000;Bechtold, 2003, Lizarralde et al., 2004).In cork oak stands, the main product is cork, which is periodically removed in harvest intervals of 9-14 years, depending on the ecological characteristics of the area.After harvesting, the cork cambium (phelogen) adds a new layer of cork to the outer bark (phellem) of the tree every year (Caritat et al., 2000).The diameter growth is thus the sum of wood and cork growth.Cork growth depends on ecological and silvi-cultural conditions and on the growth rate which varies considerably among trees (Montero and Cañellas, 1999).As a consequence of this variability, diameter over cork is not considered suitable to use as a predictor variable in growth models.Therefore, for the purposes of this study, the predictor variable used was diameter at breast height under cork.
Cork oak stands in Spain can be differentiated into open cork oak woodlands (low tree density, «dehesas») and cork oak forests (higher tree density) (San Miguel et al., 1992;Montero and Cañellas, 1999) according to ecological, silvicultural and productive characteristics.Although the main activity in open cork oak woodlands is cork extraction, they also provide grazing for domestic and wild livestock.The compatibility of these two uses is achieved by reducing the number of trees per hectare.Open cork oak woodlands occupy more than 300,000 ha in the west and southwest of Spain; they have an open structure with 20-100 trees per hectare, 10-50% canopy cover and a well developed understory of annual grasses (Montero and Torres, 1993;Montero et al., 1994).
Cork oak forests, covering a total 170,000 ha, are mainly located in Catalonia and the south of Andalusia.These forests have a higher density and a substantial understory of shrubs such as Arbutus unedo, Phyllirea latifolia, Cistus sp., Erica sp., etc. (Montero and Torres, 1993;Montero et al., 1994).
The main objective of this study is to develop a generalized height-diameter model and a crown diameter prediction model for Quercus suber L. grown in cork oak forests in Spain.Bearing in mind that the data used in this work are from different stand and regional scenarios, both equations might have a regional application.These models may prove useful not only in numerous forest management applications but also as two of the main modules of management oriented growth models, such as that developed by the authors in the CIFOR-INIA for cork oak forests in Spain.

Data
Data for developing the models were provided by the Second National Forest Inventory (2NFI) (ICONA, 1990).The 2NFI plots are systematically distributed using a grid of one square kilometre.Each plot consists of four concentric subplots with radii of 5, 10, 15 and 25 m.For each of these subplots, the minimum tree diameter recorded is 7.5, 12.5, 22.5 and 42.5 cm, respectively.In order to expand the data to the whole hectare for each minimum diameter, the following expansion factors were used: 127.32, 31.83, 14.15 and 5.09, respectively.At plot establishment, the following data were recorded for every sample tree: species, diameter over bark at 1.30 m to the nearest millimetre and total height to the nearest quarter meter.Diameters were measured with callipers in two perpendicular directions.In a second stage, three or four trees per plot were randomly selected, recording cork thickness and crown diameter among others variables for each tree.Cork thickness was obtained by averaging two perpendicular measurements taken at 1.30 m using a bark gauge.Two crown diameters were measured per tree, one being the horizontal diameter of the axis of the crown which passes through the centre of the plot and the second being perpendicular to the first.The arithmetic mean crown diameter calculated from these two field measurements is the crown diameter considered as a dependent variable.
431 plots mainly located in Catalonia and the south of Andalusia were chosen from the 2NFI database using BASIFOR software (Río et al., 2001).The plots selected met the following criteria: (1) at least 75% of the basal area was cork oak, (2) at least 50% of the number of trees per hectare were cork oak, (3) basal area above 10 squared meters per hectare, and (4) number of trees per hectare above 100 (Montero and Cañellas, 1999).
For cork oak, diameter at breast height under cork was the main predictor variable used to predict other variables at tree level.Diameter at breast height under cork can be calculated as the difference between diameter at breast height over cork and cork thickness.The last variable was measured only on three or four trees per plot, so the fitting data set is composed of 1,660 observations in the 431 plots.Table 1 shows a characterization of the data set.
In order to estimate stand variables such as dominant diameter under cork for each plot, the diameter at breast height under cork was calculated by subtracting the mean cork thickness of the three or four fullsampled trees from the diameter over cork, assuming the same cork age for all trees in each plot.This is a normal assumption in Spanish cork oak forests (Montero and Cañellas, 1999;Montes et al., 2005).Since cork age is unknown in NFI data, it is not possible to fit a function relating diameter under cork to diameter over cork.

Candidate functions
Nine generalised height-diameter equations (Table 2) were selected as candidate functions to model the height-diameter under cork relationship (Krumland and Wensel, 1988;Tomé, 1989;Soares and Tomé, 2002).All functions tested are non-linear and constrain the height-diameter relationship to pass through the point (1.30, 0) and also through the point of dominant height-dominant diameter (H 0 , D 0 ).The first constriction prevent negative height estimates for small trees, and the second ensures good predictions for larger diameters (Krumland and Wensel, 1988;Tomé, 1989;Cañadas, 2000).
The equations analysed for the crown diameter predictor model are displayed in Table 3: linear, parable, power, monomolecular and Hossfeld I.

Model fitting and evaluation
The available fitting data set consists of measurements taken from trees located within different plots.This hierarchical nested structure leads to lack of independence, since a greater than average correlation is seen detected among observations coming from the same plot (Gregoire, 1987;Fox et al., 2001).

2004), including both fixed and random component.
A general expression for a linear or nonlinear mixed effects model can be defined as (Lindstrom and Bates, 1990;Vonesh and Chinchilli, 1997;Pinheiro and Bates, 1998): where y ij is the jth observation (tree) of the response variable taken from the ith sampling unit (plot) [j=1,…n i ]; x ij is the jth measurement of a predictor variable taken from the ith plot; Φ i is a parameter vector, r × 1(where r is the number of parameters in the model), specific for each sampling unit; f is a linear or Height-diameter and crown diameter models for Spanish cork oak forests 79  et al. (1999) [h7] Mønness (1982) [h8] Michailoff (1943) modified by Tomé (1989) [h9] Prodan (1965) modified by Tomé (1989) h: total tree height (m).du: diameter at breast height under cork (cm).H 0 : dominant height (m).D 0 : dominant diameter under cork (cm).a, b: fitting parameters.
nonlinear function of the predictor variables and the parameter vector; and e ij is the residual noise term.In vector form: where y i is the (n i × 1) vector including complete observations from the ith plot [y i1 , y i2 ,...y ij ,...y inj ] T ; x i is the n i × 1 predictor vector for the n i observations of the predictor variable x taken from the ith plot [x i1 , x i2 ,...x ij ,...x inj ] T ; and e i is a n i × 1 vector for the residual terms [e i1 , e i2 ,...e ij ,...e inj ] T .
The main features of mixed-effects models are that they allow parameter vectors to vary randomly from plot to plot; regression coefficients are broken down into a fixed part, common to the population, and random components, specific to each plot.The parameter vector Φ i , can then be defined as (Pinheiro and Bates, 1998): where λ is a p × 1 vector of fixed effects, b i is a q × 1 vector of random effects associated with the ith plot with mean zero and variance σ 2 b , and A i and B i are design matrices of size r × p and r × q, for the fixed and random effects specific to each plot, respectively.In the basic assumptions, the residual within-plot errors are independently distributed with mean zero and variance σ 2 e and are independent of the random effects.The approach used in modelling variance and correlation structures is basically the same for linear mixed-effects models as for nonlinear ones.Details can be found in Lindstrom and Bates (1990), Pinheiro and Bates (1998) and Vonesh and Chinchilli, (1997).The linear mixed-effects models were fitted using the restricted maximum likelihood method implemented in the PROC MIXED procedure of the SAS/ETS software (SAS Institute Inc., 2004), while the SAS macro NLINMIX was used to fit the nonlinear models.
In those equations with more than one parameter, a determination was made as to which of the parameters in the model would be considered as parameter of a mixed effect, composed of a fixed part (common to all data in the sample) and of a random part (specific for every sampling plot), and which would be considered as parameters of a purely fixed effect (Fang and Bailey, 2001).
The evaluation of the models was based on Akaike's information criterion (AIC), the Schwarz's Bayesian information criterion (BIC), the -2 × logarithm of likelihood function (-2LL) and on numerical and graphical analyses of the residuals.Three statistics were examined: the bias, which reflects the deviation of the model with respect to observed values; the root mean square error (RMSE), which analyses the precision of the estimates; and the coefficient of determination (R 2 ).The expressions may be summarized as follows: where y i , y ˆand y ¯are the measured, estimated and mean values of the dependent variable, respectively, n the total number of observations used to fit the model, and p the number of model fixed parameters.
Another important step in evaluating the models was to perform a graphical analysis of the residuals and assess the appearance of the fitted curves overlaid on the data set.
Once the best crown diameter equation had been selected, several variables characterizing the stand were included in the mixed model as fixed effects (Hökkä, 1997;Pinheiro and Bates, 1998;Singer, 1998).Stand variables tested were basal area (G m 2 ha -1 ), number of trees per hectare N, dominant diameter (D 0 cm), quadratic mean diameter (Dg cm), dominant height (H 0 m), mean basal area of the trees larger than i tree where du i > du j (BAL m 2 /ha), and crown ratio.Criteria for including explanatory variables were the level of significance for the parameters, reduction in the values of the components of the variance-covariance matrices, significant decreases for Akaike's information criterion (AIC), the Schwarz's Bayesian information criterion (BIC) and the -2 × logarithm of likelihood function (-2LL), as well as the rate of explained variability.
The validation of the selected function for both models was done through characterisation of the model error.Since an independent validation data set was not available, the PRESS (Prediction Sum of Squares) statistics were used (Myers, 1990): [7] where y i is the observed value of observation i, y ˆi,i is the estimated value for observation i in a model fitted without this observation and n is number of observations.The bias and precision of the estimations obtained with the selected models were analysed by computing the mean of the press residuals (MPRESS) and the mean of the absolute values of the press residuals (MAPRESS), using a SAS macro.The selected models were also evaluated by examining the magnitude and distribution of the press residuals across the different predictor variables.

Height-diameter equation
The results obtained by fitting the candidate equations are shown on Tables 4 and 5. Models [h2], [h5] and [h9] did not meet the convergence criterion.As this circumstance persists when the convergence criteria are decreased or the initial parameter values are changed, these models will not be considered further.
All the parameters were found to be significant at the 5% level except parameter b in model [h1] (Table 4).Results from the comparative analysis (Table 5) suggest that the Stoffels and Van Soest model [h3] with the a parameter varying randomly between plots performed best, therefore this model was finally selected: [8] where h ij is total height of the jth tree in the ith plot (m); H 0 is dominant height of the ith plot (m), du is diameter at breast height under cork (cm), D 0 is dominant diameter under cork of the ith plot (cm), u i is the random effect associated with the ith plot with mean zero and variance 0.064 and e ij is the residual error term of the jth observation in the ith plot with mean zero and variance 1.4447.Figure 1a shows the plot of the residuals versus height estimated from the selected model.No trends were detected that suggest the presence of heteroscedasticity.
Height-diameter and crown diameter models for Spanish cork oak forests 81 the validation procedure, the mean (MPRESS) and the mean of the absolute values (MAPRESS) of the press residuals were computed for the Stoffels and Van Soest model.The values obtained, although different from zero, were small: 0.0568 m for MPRESS and 0.9698 m for MAPRESS.Plots of the mean and the absolute mean of the press residuals across the different predictor variables (Fig. 2) showed that the selected model is accurate although it tends to slightly overestimate height predictions.

Crown diameter equation
Tables 4 and 5 also show the results obtained by fitting the crown diameter models tested.The monomolecular function [cw4] did not converge so it was discarded.All the parameters were found to be signif icant at the 5% level except parameter a in model [cw2], so this model was refitted without parameter a (Table 4).Results from the comparative analysis (Table 5) indicated that the model which performed best was the parable model [cw2] without the intercept and with the b parameter divided into a fixed part and a random between-plot component.Consequently, this model was selected.Figure 1b shows the plot of the residuals versus crown width estimated from the selected model.No trends were detected that suggest the presence of heteroscedasticity.
In order to give a regional character to the selected model, several variables characterizing the stand were tested for inclusion in the mixed model as fixed effects.The best predictive capabilities were found by incorporating quadratic mean diameter as a fixed effect in the  b parameter.All parameters were significant the 0.05 level.Table 6 shows the comparison between local and generalized models.After the inclusion of quadratic mean diameter as a fixed effect, the between-plot variability decreases and the model performs more adequately, so the following model was finally selected: where cw ij is the crown diameter of the jth tree in the ith plot (m); du is diameter at breast height under cork (cm), Dg is the quadratic mean diameter in the ith plot (cm), u i is the random effect associated with the ith plot with mean zero and variance 0.0007 and e ij is the residual error term of the jth observation in the ith plot with mean zero and variance 1.0577.
For the validation procedure, the mean (MPRESS) and the mean of the absolute values (MAPRESS) of the press residuals were computed for the selected model.The values obtained were 0.0306 and 0.9689 m respectively.As can be seen from Figure 3, the precision shows no notable trend over the predictor variables.

Discussion
This study presents height-diameter and crown width prediction equations for Spanish cork oak forests based on data from the Second National Forest Inventory (2NFI) (ICONA, 1990).Due to the hierarchical nested structure of the data set, with measurements being taken Height-diameter and crown diameter models for Spanish cork oak forests 83 trees located in different plots, the multilevel mixed model approach was applied.Modelling the height-diameter relationship as a stochastic process considering random variability has been proposed by many authors since was first applied by Lappi (1997).
More recently, in Spain, has been proposed by Calama and Montero (2004) for Pinus pinea L. and by Castedo Dorado et al. (2006) for Pinus radiata D. Don.Regarding crown width, as far as we know, the mixed-model methodology has not previously been applied.The height-diameter equations tested in this are constrained to pass through the point (H 0 , D 0 ).This formulation has already been proposed in heightdiameter models by Krumland and Wensel (1988), Tomé (1989), Cañadas (2000), Calama and Montero (2004) and Diéguez-Aranda et al. (2005), and guarantees that the asymptote is near to the dominant height and that the height growth rate is smaller for the greatest dominant heights (Soares and Tomé, 2002).In addition, conditioning the model in terms of dominant trees makes it sensitive to variations in stand characteristics, being appropriate for most of the cork oak forests in Spain.
In the validation stage, the greatest prediction errors were obtained for larger diameter classes.It was also found that the prediction error increases slightly with dominant heights, in particular, for values higher than 13 m.This can be due to the very few observations of such heights in our data set.The highest site index value defined by Sánchez-González et al. (2005) was 14 m, so trees larger than 14 m are not very common in cork oak forests.For this reason, the model should not be applied in stands with dominant height values above 14 m.
The data set does not include trees smaller than 5 cm.Therefore, in order to avoid unrealistic height predictions, the developed height-diameter model must not be applied outside the diameter range 7.5 to 75 cm.The development of a specific height-diameter model would be required for cork stands in the juvenile stage (du < 7.5 cm) because periodical cork stripping is inexistent at that stage.
In order to illustrate the behaviour of the heightdiameter model, height vs. diameter at breast height under cork was plotted (Fig. 4).Height was calculated using the selected model ([h3]) for f ive dominant heights corresponding to the site index classes defined for cork oak forests in Spain (Sánchez-González et al., 2005) and considering a dominant diameter of 40 cm for all of them.As it can be observed, the curves assume biologically reasonable shapes.
The crown diameter model provides adequate crown diameter predictions for Spanish cork oak forests.The selected model, like the rest of the functions tested, uses diameter at breast height under cork as predictor variable because it is by far the most common variable used in crown diameter prediction models (Bechtold, 2003).The parable function, without the intercept and with the quadratic mean diameter incorporated as a fixed effect into the b parameter, proved to be the model with the best prediction capabilities.The signs of all parameters were consistent and biologically reasonable.Diameter at breast height under cork (du) is the strongest predictor of crown diameter for cork oak.Beyond du, moderate improvements are gained by including the mean square diameter which reduces the root mean square error from 0.9511 to 0.9332.Tome et al. (2001), in a crown diameter model for juvenile stands using the monomolecular function, expressed the shape parameter (b) as a function of stem diameter, number of trees per hectare and the ratio between the diameter at breast height of subject tree and the quadratic mean diameter of the stand, whereas for adult stands density was not included.Paulo et al. (2002) reported an improvement in a model to relate crown diameter to diameter at breast height in open cork oak woodlands by including a crown shape parameter and distance to the nearest tree.
If the present crown diameter model is compared with that developed by Tome et al. (2001) in the SUBER model for adult open woodlands, it can be appreciated that in the latter, the asymptotic value is attained at 29.93 m while in our model, considering 57 cm to be the maximum quadratic mean diameter (which is the maximum value considered for our data set), it will be attained at 16.13 m.The first value can be considered a maximum biological potential while our asymptotic value is the largest crown diameter attained for trees growing in cork oak forests with higher densities and a substantial understory of shrubs.In addition, formation and fructification pruning treatments, which are normal practice in open woodlands, have an influence on crown shape and lead to flatter crowns.These silvicultural treatments are not common in cork oak forests, where the cork growth is not modified through pruning.
Both developed models, generalized height-diameter and crown diameter prediction models, were described as a stochastic process, where a fixed model explains mean value, while unexplained residual variability is described and modelled by including random parameters acting at plot and residual levels.This approach would allow us to calibrate developed models for new locations using complementary observations of the dependent variable if available (Calama and Montero, 2004).

Conclusions
The height-diameter model developed in this study gave reasonably precise estimates of tree heights and is recommended for use in cork oak forests within the range of conditions defined above.The crown diameter model provides reliable estimations of crown width and is sensitive to quadratic mean diameter variations.Therefore, it could be used to characterize cork oak forest structure, which in turn is used to simulate stand development.Mixed-model techniques were used to estimate fixed and random-effect parameters for height-diameter and crown diameter models.The inclusion of random-effects specific to each plot allow us to deal with the lack of independence among observations derived from the special hierarchical structure of the data (trees within plots).Both models may contribute signif icantly to the integrated management model developed by the authors which can be used as an aid to define the optimum silvicultural practices for cork oak forests in Spain.

Figure 1 .
Figure 1.Plots of residuals versus predicted values for the selected models: (A) generalised height-diameter equation; (B) crown width equation.

Figure 2 .
Figure 2. Mean and absolute mean of PRESS residuals (MPRESS and MAPRESS) by diameter, dominant diameter and dominant height classes for the selected height-diameter model.The number of observations in each class is given in brackets.Dotted lines indicate standard error for the mean and dashed lines indicate standard deviation.

Figure 3 .
Figure 3. Mean and absolute mean of PRESS residuals (MPRESS and MAPRESS) by diameter and quadratic mean diameter classes for the selected crown diameter model.The number of observations is given in brackets.Dotted lines indicate standard error for the mean and dashed lines indicate standard deviation.

Figure 4 .
Figure 4. Height-diameter relationship considering a dominant diameter of 40 cm for each site quality for five dominant heights that correspond to site index classes defined for cork oak forests in Spain(Sánchez-González et al., 2005).

Table 1 .
Characterisation of data set, tree and stand variables

Table 3 .
Crown diameter functions analysed cw: crown width (m).du: diameter at breast under cork (cm).a, b and c: fitting parameters.

Table 4 .
Parameter estimates, corresponding standard errors and P-values for the models analysed

Table 5 .
Values of the goodness-of-fit statistics for fitting and cross-validation phases for the models analysed a See Tables2 and 3for forms of the functions.* Not significant (p > 0.05).

Table 6 .
Comparison of fitting statistics and estimated variance components (approximated standard errors in brackets) of the local and generalized approaches for the crown diameter model selected : diameter at breast height under cork (cm).Dg: quadratic mean diameter (cm).σ 2 : variance terms.-2LL: -2 × logarithm of likelihood function.AIC: Akaike's information criterion.BIC: Schwarz's Bayesian information criterion.RMSE: root mean squared error.R 2 : coefficient of determination. du