Height-Diameter Models for Scots Pine and Birch in Finland

Height-Diameter (H-D) models for two shade-intolerant tree species were estimated from longitudinal data. The longitudinal character of the data was taken into account by estimating the models as random effects models using two nested levels: stand and measurement occasion level. The results show that the parameters of the H-D equation develop over time but the development rate varies between stands. Therefore the development of the parameters is not linked to the stand age but to the median diameter of the basal-area weighted diameter distribution (DGM). Models were estimated with different predictor combinations in order to produce appropriate models for different situations. The estimated models can be localized for a new stand using measured heights and diameters, presumably from different points in time, and the H-D curves can be projected into the future.


Introduction
Many studies have presented models for the prediction of the height-diameter (H-D) relationship of a stand.Most of these models use a representative sample of height sample trees from the target stand (Curtis 1967, Arabatzis and Burkhart 1992, Huang et al. 1992, Lynch and Murphy 1995, Fang and Bailey 1998) and the main focus in these studies lies in finding the best functional form of the model.However, in many situations, the sample size needed is too large for practical purposes.This is because height measurements are timeconsuming and in many situations, for example in a stand wise inventory for forest management planning, only one or a few sample trees from a stand can be measured.Therefore, in recent studies, models that can predict the H-D relationship of a stand using few or even no sample trees have been developed (Lappi 1997, Eerikäinen 2003, Mehtätalo 2004).In these models, the accuracy of the prediction can be improved by enlarging the number of height sample trees.
The H-D relationship of a stand is not stable but develops over time (Curtis 1967, Flewelling and de Jong 1994, Lappi 1997).However, Mehtätalo (2004) observed that with a shade-tolerant tree species (Norway spruce, Picea abies) the age at which maturity is reached varies from stand to stand and he gave two reasons for this.First, after remaining as undergrowth for a long time, shade-tolerant trees may suddenly begin to grow as rapidly as young saplings and secondly, the development of a stand from a sapling stage to a mature stand takes longer on poor sites than on rich sites, i.e. the development rate varies between stands.On the other hand, stands that are reaching maturity seem to have almost equal mean tree size.Thus, instead of linking the development of H-D relationship with stand age, Mehtätalo (2004) linked it with basal area weighted median diameter of the stand (DGM).The latter of the reasons given above for the variation in the age at which maturity is reached might hold also with shade-intolerant tree species.Therefore, this study presents an analysis similar to that in Mehtätalo (2004) but with two shade-intolerant tree species.
The aim of this study is to model the H-D relationship of Scots pine (Pinus sylvestris) and birch (Betula pendula and Betula pubescens) in Finland.The methodology is the same as was used in Mehtätalo (2004), but here it is applied to shade-intolerant tree species.Methodologically, the aim is to study if the development of the H-D curve of shade-intolerant tree species should also be linked with mean tree size rather than with stand age.Furthermore, the aim of this study is to show that the model formulation used in Lappi (1997) and Mehtätalo (2004) can be successfully applied with several tree species with only minor changes in the model formulation.

Data
The modeling data are a subset of a larger dataset collected by the Finnish Forest Research Institute (Gustafsen et al. 1988).The sample stands of the data were selected randomly from those sample plots of the 7th National Forest Inventory which are situated on mineral soils and forest land.Three fixed-radius sample plots were established in all stands.Each sample plot was measured 3 times with 5 year intervals.In addition, when establishing the plots, the growth of the trees over the previous 5 years was recorded in some stands.Thus, the number of measurement occasions for each stand varied from 1 to 4.
In this study, the three plots were combined to obtain the data of a stand.Only stands with on average more than 10 Scots pines / measurement occasion were selected to the Scots pine data and stands with more than 8 birches / measurement occasion to the birch data.All trees of other tree species than the one in question were removed from the data of any given tree species.However, before doing this, DGM and basal area were calculated from all trees belonging to the dominant story of the stand.The Scots pine data included 46338 observations from 497 stands (1774 measurement occasions) and the birch data 2979 observations from 61 stands (190 measurement occasions).Table 1 summarizes the modeling datasets of this study.

Model Development
The model is the linearized form of the exponential function (Korf function), which in many studies has been found good for the description of the H-D relationship of a stand (Zakrzewski and Bella 1988, Huang et al. 1992, Arabatzis and Burkhart 1992, Fang and Bailey 1998).The model formula is where H is tree height and D diameter at breast height and A, B and C and λ are parameters.The parameter λ is interpreted as the expected difference between diameter at ground level and that at breast height.It is an alternative for the more commonly used subtraction of the breast height from the height; for more discussion on the matter see Lappi (1991Lappi ( , 1997)).
Trying to fit model (1) to the data showed that the model is clearly overparametrized.Thus, the first step in the analysis was to reduce the number of parameters to be estimated.To do this, model (1) was fitted separately for each stand and measurement occasion using different combinations of parameters C and λ.The value of λ giving the lowest mean error variance over all stands and measurement occasions was selected as the fixed value of λ (see Mehtätalo 2004).The selected values were 7 for Scots pine and 6 for birch.For each stand and measurement occasion, the value of C giving the lowest error variation was selected and these values were modeled as a function of DGM to fit a heuristic trend function for parameter C (see Lappi 1997).The trend function of C for Scots pine was C = 0.9823 + 0.05753 × DGM but for birch no trend was found and thus the constant C = 1.809 was used.
The model for tree i in stand k at time t is where parameters A kt and B kt need to be estimated.The next step in the analysis was to study whether DGM is a better descriptor of the stage of development than the stand age.If it is, then in the subsequent longitudinal analysis (Diggle et al. 2002) rather than using stand age, we use DGM as the variable describing the development of the stand.Model (2) was fitted separately for each stand and measurement occasion.The obtained parameter estimates were plotted against stand age and DGM (Fig. 1).Only the estimates of parameter A are shown here; it is interpreted as the asymptote of the H-D curve of a stand, i.e. it is the maximum tree height of the stand.In any one stand, parameter A seems first to develop rather rapidly and later level out at some level, which can be interpreted as the maximum height of a tree of that tree species growing at that site (Figs.1a  and 1b).This level varies considerably between stands.Another, more interesting feature is that in those stands where the maximum height is low, the overall development of the stand takes longer than in those stands where the maximum height is high.In other words, the development rate of parameter A varies between stands, being higher in stands where trees get higher.This implies that two stands of equal age are not at the same stage of development.Plotting the parameters against DGM (Figs. 1c and 1d) shows that the form of the trend as a function of DGM is quite similar in all stands except for a vertical shift in the level of the curve.This implies that stands with equal DGM are at the same stage of development, but, because of site properties, location, etc., the asymptote of the H-D curve is not equal in all stands with any given DGM.Thus, a good strategy in modeling is to link the H-D relationship with DGM in the model and take the vertical shift in the parameters into account with stand-specific predictors and random parameters.The estimates of A and B and their estimation errors are strongly correlated.To reduce these correlations, the diameter was reparametrized as The model using this parametrization is ln ( ) Parametrization (3) provides interpretations for the parameters of (4): A is the expected logarith-mic height of a tree with a diameter of DGM + 10 and B is the expected difference in the logarithmic height between diameters of 30 and 10 cm.The next steps in the analysis were as follows (for more in-depth description, see Mehtätalo 2004 and Lappi 1997): 1) Model 4 was fitted with OLS for each stand and measurement occasion.2) Appropriate trend functions for the estimates of parameters A and B were searched for and fitted as multilevel random effects models using stand and measurement occasion levels.
3) The trend functions were written into Model 4. The nonlinear parameters were used as fixed constants and the linear parameters were re-estimated.An appropriate model for the residual variation as a function of tree diameter was defined and fitted concurrently.4) Step 1 was repeated using WLS, where the weights were calculated as the inverse of the variance function obtained in step 3. Furthermore, steps 2 and The analysis was carried out with the R-implementation of the S-language (Chambers 1998, Venables and Ripley 2002, see http://www.r-project.org),where package nlme (Pinheiro and Bates 2000) was used for the random effects models.

Fitted Models
For both tree species, an appropriate trend function for parameter A was the Chapman-Richards function (Richards 1959) and for parameter B, a linear function of the form was used.In ( 5) and ( 6 where z kt is the nonlinear part of ( 5) and the parameters are as explained before.It is assumed that the random effects are normally distributed with a mean of 0 and constant variance.Covariances cov(α k , β k ) and cov(α kt , β kt ) may be nonzero but all other covariances between the random effects and the error term are zero.The error term ε kti is assumed to be normally distributed with a variance that depends on tree diameter according to the formula var max , ( ) where Δ = 4.5 for Scots pine and Δ = 9.5 for birch.These values were determined by visually exam-ining a figure where the means and standard deviations of the OLS residuals were plotted in one cm diameter classes (see Lappi 1997) The estimates of the fixed parameters and variances of the random parameters are in Tables 2  and 3.For both tree species, 5 models with different predictor combinations were estimated.Model I uses only DGM as the predictor.In model II, variables describing the geographical location are included.Models III and IV include, in addition to the predictors of model II, stand characteristics measured from the whole growing stock and furthermore, model V includes those measured by tree species.The difference between models III and IV is that model IV includes stand age.Only predictors with statistically significant coefficients were included.The significance level used was 1% in Scots pine models and 5% in the models for birch.The significant predictors were searched stepwise, refitting the model and eliminating the least significant parameter until all remaining coefficients were significant.
One can see that the variances of the random parameters decrease when the number of predictors increases.This happens because the variation explained by the additional fixed predictors belongs to the random part in the more sparse models.where measured heights are in vector y, the prediction of the fixed part in vector µ, random parameters in vector b and residuals in vector e.Vector b includes stand effects and time effects for those points in time from which we have measurements.
The number of time effects depends on how many points in time the measurements have been taken from.It is convenient to write the time effects after

Prediction of DGM in the Future
In the prediction of the H-D curve of a stand at a given point in time, the DGM of the stand at that point in time needs to be known.With current and past points in time this is not a problem, since in Finnish inventories DGM is a very commonly measured mean stand characteristic.However, in order to be able to predict the H-D curve in the future, a model for the DGM was estimated.Since the DGM used in the models is the DGM of all trees species belonging to the dominant canopy layer, one and the same model can be used with Scots pine, birch and Norway spruce.To construct such a model the modeling data comprised both the full data used in this study and that used in Mehtätalo (2004).However, if a cutting had been carried out in the stand during the last 5 years prior to the measurement, the observation was deleted from the data and either the time series before or that after the cutting, whichever was longer, was used.
As seen in Fig. 2a, the DGM of a stand with a given age varies considerably, especially with older stands.However, the change in DGM is stable with respect to stand age.Thus, if we have an observation of stand DGM at some age, we can predict the DGM at some other age rather well.The nonlinear trend in Fig. 2a where T kt is age of stand k at time t, u and v are fixed parameters, u k is a stand level random parameter and e kt is the residual.The random parameters and residual error were assumed to be normally distributed with a constant variance.The parameter estimates obtained were u = 0.3260 (s.e.0.03816) and v = 0.6321 (s.e.0.009543), var(u k ) = 0.3212 2 and var(e kt ) = 0.06856 2 .The linearized model seems to fit the data well (Fig. 2b).When utilizing the model to predict the DGM of a stand with a given age, the stand effect u k is predicted using the best linear unbiased predictor (Eq.12).When the predicted logarithmic DGM is trans-formed back to the arithmetic scale, the bias correction needs to be applied (Eq.16).

Application Example
To demonstrate the use of the estimated models, H-D curves were predicted with each of the models for one stand selected from the modeling data.The stand was a 35-year-old mixed-species forest with a DGM of 12 cm and a basal area of 22m 2 /ha, of which 17.7m 2 /ha was Scots pine and 3.4m 2 /ha birch.The correction for bias was applied to all predictions (Equation 16).Fig. 3 shows the fixed part predictions obtained with each of the models in Tables 2 and 3.The predictions obtained with the various models differ slightly from each other.However, for both  tree species all models give too low heights in this stand.Fig. 4a demonstrates the effect of localization on the predictions.Models II and V for Scots pine were localized for the sample stand by predicting the stand and time effects of the H-D models using one measured height sample tree from the stand (Eq.12).The localized models give much better height predictions than the fixed parts only.Note that in Fig. 4a the localized models are very close to each other even if the fixed part predictions are not.This is because the information of one measured sample overrides the information of the additional fixed predictors of model V.This demonstrates the somewhat self-evident fact that it is more efficient to improve the prediction of the H-D curve of a stand by measuring heights and diameters than by measuring the covariates of model V.
Fig. 4b demonstrates the projection of the H-D curve of a stand into the future.The measurements were made at the age of 35 years and the H-D curve is projected to the age of 45 years.This required knowledge about the DGM at the age of 45 years, which was predicted using model ( 17).The random effects of model ( 17) were predicted with BLUP using the known DGM at the age of 35 years to obtain a stand-specific DGM-age -curve.It was used to predict the DGM of the stand at the age of 45 years.The projected H-D curves were calculated using the predicted DGM.The projections obtained with the localized model are again clearly better than the predictions of the fixed part only but they seem to be underestimates of the height.In fact, all localized models in Fig. 4 seem to give slightly underestimated heights.This is because the expected H-D curve is so far from the observations that one measurement does not move it far enough.Using more than one sample tree would reduce the bias of the predictions.

Discussion
In this study H-D models for Scots pine and birch were estimated for practical use in Finland.The models can be used to predict the H-D relationship of a stand with known DGM.If other stand characteristics than DGM are measured, they can be used through selecting from models I,…,V the model that best suits the situation.Measured height-diameter pairs can be used in localizing the model for a target stand and, due to the longitudinal character of the model, information from any point in time, i.e. any stage of development, can be utilized.Furthermore, the H-D curve of a stand in the future can be predicted, but this requires the prediction of DGM at that point in time.These properties of the model are discussed in Lappi (1997) and Mehtätalo (2004).
The number of parameters in the model is quite high.Some of the parameters are so called second level parameters, which are included in order to improve the properties of parameter estimates and to provide interpretations for the actual parameters of the model.The statistical significance of these parameters was not tested and a more parsimonious model might have lead to an equally accurate prediction.However, because the aim of this study was predicting, not studying the effects of different factors on the height-diameter relationship, and because the second level parameters did not cause any harm in the estimation phase and the models also worked well in the prediction, there was no need to decrease the number of parameters.Mehtätalo (2004) observed that the development of the H-D curve of a shade-tolerant tree species (Norway spruce) depends on mean tree size in the stand rather than on stand age.The present study continued this work and showed that the observation made with shade-tolerant tree species is true also with shade-intolerant tree species.The reason for this is that the site properties affect the development rate of a forest stand, so that stands on poor sites develop more slowly and for longer than stands on rich sites.Because of this, the models predicting the development of the H-D curve of a stand perform better when mean diameter of the stand is used as the variable describing the maturity of the stand instead of stand age.This effect on the performance is probably also true of models predicting other things than H-D curves.Hence, when modeling the development of any stand characteristics, for example, diameter distribution and stand growth, the use of stand age as the only variable determining the stage of development of the stand should be viewed critically.
However, the method for taking into account the effect of DGM on the development of H-D curves which has been presented here is not the only correct one; other approaches may also lead to equally good results.For example, the development of parameters A and B can be linked with stand age and the DGM can be used in the prediction of the random effects as Lappi (1997) did in his models 6 and 7. Thus, if both DGM and age are known, it is not obvious that the development of H-D curves should be linked with DGM and linking it with age may work as well, if the DGM is taken into account in the models.However, using age only will not lead to as good models as will the use of DGM only and thus, if one of them must be chosen, it would be preferable to use the DGM.

Fig. 1 .
Fig. 1.Estimates of parameter A of model (2) against stand age (a and b) and stand DGM (c and d) in Scots pine data (a and c) and birch data (b and d).The estimates of the same stand at different points in time are connected by lines.
In step 5, the fixed parameters p 1a and p 1b were written as linear combinations of the predictors used in the model, i.e. p 1a = b 0a + b 1a x 1a + b 2a x 2a + … and p 1b = b 0b + b 1b x 1b + b 2b x 2b + …, where x 1a , x 1b , x 2a , x 2b ,… are the additional predictors and b 0a , b 0b , b 1a , b 1b ,… are their coefficients.

Fig. 2 .
Fig. 2. Stand DGM against stand age in the combined data.Subfigure a shows the relation on the arithmetic scale and subfigure b shows the relation between logarithmic DGM and age.The solid line in subfigure b is the expected development obtained with the fixed part of model (17).

Fig. 3 .
Fig. 3. Predicted H-D curves of Scots pine (a) and birch (b) in a 35-year-old mixed pine-birch stand selected from the modeling data.The predictions are calculated with each of the models I-V using only the fixed part of the model.The marks show the observed heights and diameters.

Fig. 4 .
Fig. 4. Predicted H-D curves of the Scots pines in the stand of Fig. 3. Dashed lines are the predictions of the models obtained using the fixed part only and solid lines are localized using one measured height-diameter pair at the age of 35 years (the big solid ball in the figures).Subfigure a shows the predictions obtained using models II and V at the age of 35 years.In subfigure b, the H-D curves after 10 years growth are predicted with model II.Small symbols are the observations (• at the age of 35 years and ▲ at the age of 45 years).The DGM of the stand at the age of 35 years was 9.6 cm and the predicted DGM at the age of 45 years was 11.4 cm.

Table 1 .
Some characteristics of the modeling data.

Table 2 .
The parameter estimates for the model of Scots pine.
Note: The predictors are as follows: yk, the north coordinate in the Finnish uniform coordinate system, 1000 km; alt,

Table 3 .
The estimated models for birch.