Comparing Regression Estimation Techniques when Predicting Diameter Distributions of Scots Pine on Drained Peatlands

We compared different statistical methods for fitting linear regression models to a longitudinal data of breast height diameter (dbh) distributions of Scots pine dominated stands on drained peatlands. The parameter prediction methods for two parameters of Johnson’s SB distribution, fitted to basal-area dbh distributions, were: 1) a linear model estimated by ordinary least squares (OLS), 2) a multivariate linear model estimated using the seemingly unrelated regression approach (SUR), 3) a linear mixed-effects model with random intercept (MIX), and 4) a multivariate mixed-effects model (MSUR). The aim was to clarify the effect of taking into account the hierarchy of the data, as well as simultaneous estimation of the correlated dependent variables on the model fit and predictions. Instead of the reliability of the predicted parameters, we focused on the reliability of the models in predicting stand conditions. Predicted distributions were validated in terms of bias, RMSE, and error deviation in the generated quantities of the growing stock. The study material consisted of 112 successively measured stands from 12 experimental areas covering the whole of Finland (total of 608 observations). Two independent test data sets were used for model validation. All the advanced regression techniques were superior to OLS, when exactly the same independent stand variables were included. SUR and MSUR were ranked the overall best and second best, respectively. Their ranking was the same in the modeling data, whereas MSUR was superior in the peatland test data and SUR in the mineral soil test data. The ranking of the models was logical, but may not be widely generalized. The SUR and MSUR models were considered to be relevant tools for practical forest management planning purposes over a variety of site types and stand structures.


Introduction
The structure of a forest stand in terms of its breast height diameter (dbh) distribution is of great importance.In practical forestry, the dbh distribution is useful for estimating e.g.stand total volume, the quantity of timber assortments, and the developmental stage of the stand.Furthermore, it enables prediction and simulation of the future yields of wood and the target stand states for management objectives, e.g.cutting regimes (Hyink andMoser 1983, Franklin et al. 2002).The use of tree-specific models in growth simulators requires that the diameters are known or can be predicted using standard stand characteristics (Bailey andDell 1973, Päivinen 1980).During the last couple of decades considerable modeling efforts have been made in forest research on predicting stand structure.In Finland, for example, models for predicting the basal area-dbh distribution of Scots pine, Norway spruce and silver birch dominated stands (the commercially most important tree species in northern Europe) have been presented by e.g.Maltamo et al. (1995), Haara et al. (1997), Maltamo (1997), Siipilehto (1999), and Kangas andMaltamo (2000a, 2000b).Both parametric methods (e.g.Beta, Weibull and Johnson's S B functions) and non-parametric methods (e.g.k-nearest neighbour, kernel function or percentile based approach) have been applied in numerically describing the dbh distributions.
The dbh distribution models have mainly been constructed for stands growing on mineral soil sites.However, stands growing on drained peatlands are also very important natural resources in Finland, where about 4.7 million hectares of peatland (about 52% of the total peatland area) have been drained for forestry purposes in order to increase wood production (Hökkä et al. 2002).Scots pine is one of the most common tree species on both drained and on pristine peatlands.It is the dominant tree species on about 3.2 M. ha, i.e. it covers 68% of the total drained forested peatland area in Finland.According to recent scenarios based on data from the 8th National Forest Inventory (NFI), the annual cutting possibilities (mainly thinnings) may increase up to 15-20 M. m 3 in the course of the next 20 years (Nuutinen et al. 2000).
On pristine peatland, the stands are often sparse and they have a heterogeneous age, size and spatial structure (Gustavsen and Päivänen 1986, Norokorpi et al. 1997, Macdonald and Yin 1999).After drainage, this structural inequality may remain or even increase during the first few decades owing to the increase in stand stem number resulting from the improved regeneration and growing conditions for the trees (Hökkä and Laine 1988, Sarkkola et al. 2003, Sarkkola et al. 2005).However, during the course of post-drainage succession the impact of increasing inter-tree competition decreases the stem number even if no cuttings were carried out.Furthermore, the initial pre-drainage stand properties and the ecohydrological conditions of the site affect the stand structure for a considerable period of time after drainage (Sarkkola et al. 2005).Consequently, the initial decreasing diameter distributions gradually become less of a skewed, bell-shaped form.However, high variation among the tree size dimensions may still remain (Sarkkola et al. 2005).
Despite the increasing significance of peatland forests as an important wood resource, relatively little attention in research has been paid to stand structure.Some Weibull distribution models for predicting the stand structure on drained peatlands have earlier been presented: Hökkä et al. (1991) for Scots pine or pubescent birch dominated stands in northern Finland, Sarkkola et al. (2003) for Norway spruce, and Sarkkola et al. (2005) for Scots pine dominated stands in southern Finland.Furthermore, some other models have been constructed from data representing both mineral soils and peatland stands.For example, Mykkänen (1986) and Kilkki et al. (1989) presented Weibull distribution models based on NFI data, and Kangas and Maltamo (2000b) constructed distribution-free models for local upland and peatland forests located in central and eastern Finland.However, there are still no valid specific models for drained peatlands that would be simple to apply in practical use, but simultaneously also applicable for describing the wide range of variation in stand structures and accurately predicting the dbh distribution.
When using a parametric approach in describing the stand dbh distribution, there are two main ways of predicting the diameter distribution of a stand by using measured mean and sum stand characteristics only.In the parameter prediction method (PPM), a priori estimated regression models are applied for prediction of the distribution of the target stand (e.g.Schreuder et al. 1979).The other possibility would be to use the parameter recovery method (PRM) in which the relationships between distribution parameters and stand variables, moments, or percentiles are derived in a closed form, and the parameters estimated for the target stand are solved from the resulting system of equations (Hyink and Moser 1983).The variables used in recovery may be based on direct measurements or regression relationships.It should be noted that unless regression relationships are utilized, PRM is possible only for as many parameters as there are known distribution-related stand variables, and only stand variables that are mathematically related to the diameter distribution can be used.
In most of the previously published studies on PPM, the models for single parameters have been estimated separately, even though some recent papers have utilized some more advanced approaches (e.g.Kangas and Maltamo 2000b, Cao 2004, Liu et al. 2004, Robinson 2004, Newton et al. 2005).When fitting parameter prediction models to a longitudinal data, there are two points that should be considered in the selection of the fitting method and modeling approach.First, the model can be multivariate, i.e. several models should be fitted simultaneously to the same dataset, because the fitting of the individual models separately might be not an efficient estimation approach due to the correlation between different models (Zellner 1962).Secondly, there may be autocorrelation between successive observations.Longitudinal data are a special case of hierarchical data.Since the path-breaking paper of Laird and Ware (1982), the mixed-effects modeling approach has become the standard approach with longitudinal data.The recent development of statistical packages has provided possibilities to take into account the autocorrelation between models and between explanatory variables.
The aim of this study is i) to construct linear regression models for two S B distribution parameters estimated for the dbh distributions of Scots pine dominated stands on drained peatlands, ii) to compare the reliability of different statistical methods for fitting these models to a longitudinal data taking either 1) the hierarchy of the data or 2) the cross-model correlation, or 3) both of these properties into account on model fitting and predictions.The regression approaches were validated in terms of the accuracy of the predicted basal area-dbh distributions in generating the quantity of the growing stock.Independent test data were used for model validation.Furthermore, estimates of the most relevant earlier models constructed for Scots pine stands were compared to our ones.The best performing models will be recommended as a tool for practical forest management planning.

Study Material
The material of this study consisted of 112 successively measured sample plots (stands) from 12 experimental areas of the Finnish Forest Research Institute (Metla) (91 plots) and the University of Helsinki (21 plots) that cover the whole Finland.They were established on drained peatlands during 1925-1973.The original aim of these sample plots have been to monitor the effects of drainage and to compare the growth and yield of thinned and unmanaged stands (Lukkala 1929) or educational purposes in order to monitor postdrainage stand growth and yield, as well as to demonstrate the effects of silvicultural operations on stand development (Sarkkola and Päivänen 2001).The sample plots used in this study were subjectively selected from the whole data set (about 500 plots) choosing long-term developmental series of stands dominated by Scots pine growing on drained peatlands.The stands represent the first post-drainage tree generation with varying densities and site fertility.Even though our stands do not represent a nationwide probability sample, we assumed that they represent well the structural variation of Scots pine stands in drained peatlands in Finland.For details of the selection, see Sarkkola et al. (2005).
The stands were located in Finland within an area delimited by 60°01´-67°10'N and 23°07´-26°40´E.The mean annual temperature sum (the sum of daily average temperatures exceeding +5°C, based on hourly measurements) varies between 1400 and 800 degree days.The annual precipitation varies between 750 and 600 mm, of which 200-350 mm fall as snow.
The material covered all the peatland site types that are naturally favoured by Scots pine and primarily those drained for forestry purposes in Finland (classification according to Laine 1989): -Vaccinium myrtillus (MT II) site type characterized by V. myrtillus L., V. vitis-idaea L., and herbs of mesic sites, e.g.These site types form a nutrient regime gradient ranging from minerotrophic to ombrotrophic sites.Type I sites represent forested minerotrophic (VT I) or ombrotrophic (DsT I) pine peatlands consisting of Sphagnum spp.residue-dominated peat.Type II sites represent minerotrophic originally treeless sites or sparsely forested composite pine peatland sites, which have the properties of both treeless and fully-stocked sites typically consisting of Carex dominated peat.The main differences between type I and type II sites in terms of the development of Scots pine stands after drainage have been dealt with in detail by Hökkä and Ojansuu (2004) and Sarkkola et al. (2005).
The stands had been managed by silvicultural (light) thinnings, with the exception of 18 stands that had been left unmanaged.The initial drainage (ditching) of the sites had been performed between 1909 and 1955.The oldest sample plots were established in 1928, and had 11 repeated measurement rounds up until year 1996.The average number of repeated observations was six.The period between adjacent measurements was typically either 5 or 10 years but varied from 3 years to 25 years.The longest monitoring period of stand development on a given sample plot was 66 years.
All the trees growing on the sample plots, which varied in size from 300 m 2 to 2500 m 2 , had been measured for dbh, the smallest dbh being 1 cm.The undergrowth, primarily consisting of Norway spruce or pubescent birch that appeared particularly on type II sites, was not included in the studied dbh distributions.Some average standard stand characteristics are presented in Table 1.The total number of observations, i.e. the total number of successive stand measurements (dbh distributions) was 608.The average number of tallied trees was 382 (range 22-5406) per dbh distribution.

Test Material
Two independent data sets were used in testing the models constructed for dbh distributions.One of the test materials consisted of 52 Scots pine dominated, permanent sample plots located on drained peatlands in 14 different peatland areas in central and northern Finland (HARKO-data).The stands belong to the specific experimental designs objected to study the effects of thinning on stand growth and yield.Each thinning experiment contains a control and three different thinning treatments of varying intensity.For the testing dataset, we chose one sample plot (stand) representing each of the treatments from each of the experimental areas.The stands were measured in 1992-1996, and the smallest measured dbh value was 5 cm.The size of the selected sample plots within a test site unit varied between 600 m 2 and 1200 m 2 .Some standard site and stand characteristics of the test site units are presented in Table 2.The HARKO data covered the same peatland site types as the modelling data.However, these stands were mainly located in different geographical and climatic regions than the modeling data.In contrast to the modeling data, the HARKO data included a relatively large range of thinning intensities.In addition to the unthinned controls, relatively heavy thinnings were also included (thinning removal was at its greatest about 80% of the total basal area).Thus, this data provided important additional information about the validity of the models.The total number of test stands (sample plots) was 531.
The second test material consisted of Scots pine dominated permanent sample plots established on mineral soil sites located throughout Finland (socalled INKA data).The stands are based on a subsample of the NFI data, and they were originally established to study tree growth (Penttilä and Honkanen 1986).The INKA sample plots consisted of a cluster of three circular plots within a stand.The total number of tallied trees was about 100-120 per stand (see Gustavsen et al. 1988).For the purpose of the present study, the cluster of three plots was combined in order to obtain reliable dbh distributions.The second round of INKA measurements was used as cross-sectional test data.The stands were measured in 1976-1983.The average number of trees (only pines were included) in the INKA stands was 79 stems/plot, ranging from 11 to 135 stems per sample plot.The total number of stands was 505.
The reason for including mineral soils was that previous dbh distribution models have been applied to both mineral soils and peatlands, and it would therefore be possible to compare the behaviour of the dimension variation between peatland and mineral soil stands.Johnson (1949) presented a family of parametric distributions, which are based on transforming the distribution of the original variable into a standard normally distributed variable by applying different transformation functions.In our case, the dbh distribution for Scots pine was described using Johnson's S B function which is, according to Hafley and Schreuder (1977), the most flexible parametric distribution together with the beta function.Johnson's S B probability density function (1) applies transformation function (2), which are

Characterizing the Diameter Distribution
in which γ and δ are shape parameters, x and λ are location and scale parameters, d is the diameter observed in a stand plot.The parameters were solved using the method of conditional (x = 0) maximum likelihood (ML), as in the study of Schreuder and Hafley (1977) with the exception of basal-area-weighting (Siipilehto 1999).The average parameters of the estimated S B distribution of these data are shown in Table 1.Prediction models were formulated for parameters λ and δ.For the third parameter, a recovery approach was applied by solving it from the Eq. 3 by setting the observed basal area-median diameter (d gM ) equal with the median of the predicted basal area-dbh distribution.(3)

Regression Methods
We applied the direct parameter prediction approach in modeling the S B scale parameter λ and the shape parameter δ (Siipilehto 1999, Liu et al. 2004, Newton et al. 2005).Thus, the distribution was first fitted to the data using conditional Maximum Likelihood (ML), and the obtained estimates, treated as true values for the stand, were then modeled against stand variables.The linear model can be written as where vector y includes the responses, matrix X the predictors, vector b the coefficients to be estimated and vector e the residuals (McGulloch and Searle 2000).If the errors are independent and have equal variances, then the efficient estimator of the coefficients is the Ordinary Least Squares If the OLS assumptions of errors do not hold, then the OLS estimator is no longer efficient.Allowing heteroscedasticity and correlation between observations (i.e.assuming that var e = σ 2 V, where V is any positive definite matrix), the efficient estimator is the Generalized Least Squared (GLS) estimator When several distribution parameters are modeled in the PPM approach, one has a multivariate model (Kangas andMaltamo 2000b, Robinson 2004).We have two models of form (4) to be estimated from the same data with responses, design matrices, parameters and residual vectors of y 1 and y 2 , X 1 and X 2 , b 1 and b 2 , and e 1 and e 2 , respectively.The multivariate model can be written as or equivalently which shows that it is a special case of linear model ( 4).This model has a Seemingly unrelated regression (SUR) model structure.In this model, estimating each of the component models separately would lead to efficient estimation if the residuals of the separate models were not correlated, or if the design matrices X 1 and X 2 are the same in both models (Zellner 1962).In our case, the residuals are correlated and design matrices are unequal.Thus, the approach of Zellner (1962) is used, where the dispersion matrix of the data, var e*, is first estimated from the residuals of separate OLS fits, and b* is then estimated by GLS.
Another feature causing violation of the OLS assumptions is the hierarchy of the data.In our case, the hierarchy arises from that we have repeated measurements from several stands, which are a sample from a population of stands.A model accounting for the hierarchy is the mixedeffects model (Laird and Ware 1982) where Z is the design matrix of the random part, c the vector of the stand specific random parameters and u the vector of measurement occasion specific residuals.By defining e = Zc + u, we see that model ( 6) is a special case of model ( 4 In a repeated measurement (longitudinal) data with several responses, an efficient estimation method should be able to take into account both the hierarchy (autocorrelation) of the data and the correlations between the models.Thus, a multivariate mixed-effects model, combining the mixed-effects modeling and SUR approaches, would be appropriate.A model of this kind can be treated as a special case of a hierarchical model, where an additional level of hierarchy is added to the model of the longitudinal data (Snijders and Bosker 1999) and the assumptions about between-model covariances are parameterized in the variance-covariance-matrix of the observations.For example, the MLwiN software uses this approach in fitting a multivariate mixed-effects model (Rasbasch et al. 2004) and it leads to efficient estimation that accounts for both the data hierarchy and the between-model correlations (Goldstein 1995).
In this study, four alternative regression methods were applied in predicting the S B distribution parameters λ and δ.The regression approaches were as follows: 1) Linear model (model 4) estimated by ordinary least squares (OLS) 2) Multivariate linear model (model 5) estimated using the seemingly unrelated regression approach (SUR) 3) Linear mixed-effects model with random intercept (model 6) (MIX).4) Multivariate mixed-effects model estimated as a mixed-effects model with an additional level of hierarchy to allow simultaneous estimation (MSUR).

Model Construction
The stand characteristics d gM , G, and N were linked together by performing a transformation called shape index (see Siipilehto 1999).This shape index is given by the symbol ψ in formula (7).
The shape index (ψ) can be determined also as a squared proportion of the quadratic mean (D q ) and basal area-median diameter (d gM ).Typically, this proportion is less than 1.0, which means that d gM is greater than D q (see Table 1 and Fig. 1).The lowest shape index values are found from bimodal or descending dbh-frequency distributions, while the greatest values result in from peaked dbh distributions with a relatively high mean dbh (see Siipilehto 1999).With the shape index we can describe, for example, the variation in N conditional to fixed d gM and G. Especially large variation is found in stands with a small median diameter (d gM < 15 cm), as shown in Fig. 1.For example, considering a stand with d gM = 5 cm and G = 5 m 2 ha -1 , the shape index values of between 0.2 and 0.9 represent a difference of about 10 000 trees per hectare.The average shape index values of these data are shown in Table 1.
The shape index (ψ), stand basal area (G), basal area median diameter (d gM ), stem number (N) and annual temperature sum (DDY) were used as the explanatory factors in the models of the S B parameters.Furthermore, variables describing site type (e.g.site type dummies formed by main site type groups, i.e. for type I-and type II sites) and performed cuttings (thinning dummy and thinning intensity as a continuous variable) were tested, but they were not significant in the models.The final model formulations were chosen from many alternative formulations through intensive study of residuals.No significant heteroscedasticity was observed in the residuals, but the distribution of the residuals of ln(λ) was slightly skewed to the right, resulting from the nature of the maximum diameter.As taking the skewness into account would have required a generalized linear mixed modeling approach, it was not accounted for in this study.In all the approaches, the fixed parts of the models ln(λ) and ln(δ) were as follows: Using the notation of Eq. 6, this implies that The MIX and MSUR models were estimated using MLwiN package (Rasbash et al. 2004) applying stand-level random effect on the constant.The experimental area was also tested as an additional level of hierarchy, but it was found insignificant.In MLwiN, the restricted iterative generalised least square (RIGLS) method was used to produce unbiased restricted maximum likelihood (REML) estimates for the parameters (Rasbash et al. 2004).OLS and SUR models were estimated using the SAS package by applying the MODEL procedure (SAS/ETS... 1993).

Model Evaluation
The accuracies of the constructed models were validated in terms of bias (10) and RMSE (11).
The model precision (s b , Eq. 12) was expressed as the standard deviation of the relative error in the estimated Y i excluding the effect of bias (e.g.Siipilehto 1999).Relative bias and RMSE were calculated by dividing RMSE by the mean value of the observed Y i , and they were expressed in percentages.The estimates for the stem number and the third and fourth powers of the cumulative frequencies of the dbh's (Sd 3 , describing the stand volume, and Sd 4 describing the stand value) were compared with the values derived from the original dbh measurements.(NOTE: the cumulative frequency of the second power of the dbh 's can be derived exactly from the basal area).A numerical solution by 1 cm-dbh classes was applied when transforming the basal area-dbh distribution into the dbh-frequency distribution.The advantage of Sd 3 and Sd 4 is that they do not require height information, while they can still provide reasonable estimates of the accuracy in the 'volume' and in the 'value' of the growing stock, respectively (see, Kilkki andPäivinen 1986, Maltamo et al. 1995).
in which Y i is the observed and Ŷi is the predicted stand characteristic, and e i is the relative prediction error (%) in stand i.
The specific ranking method was used in comparing the "goodness" of the models.The lowest test values given by formulae (10-12) were given the rank score of one, the second lowest score of two, and the highest value score of four.Thus, the lowest rank score sum belonged to the best performing model.The minimum possible score was 27 while the maximum was 108, if the issued model was always the "best" or "worst" in each test criterion and data set, respectively.(I.e. 4 models (score 1-4) × 9 criterions × 3 data sets).
Two earlier models were included here in order to compare the results in exactly the same way.The models were the Weibull distribution for Scots pine by Mykkänen (1986) and the percentile-based basal area-dbh distribution including N as a predictor by Kangas and Maltamo (2000b).These models were chosen because the modeling data included stands on both mineral soils and peatlands.The models of the present study are theoretically applicable to mineral soil stands because no peatland-specific independent variable was included in the final models (Eq.8 and 9).

Models
In the present paper the Johnson's S B distribution was used to describe the basal area-dbh distributions of Scots pine dominated stands growing on drained peatlands.The shape index (1) was curvilinear in relation to the parameters of the S B distribution (Fig. 2).In this case the shape index proved to be extremely useful, while its linearized (and normalized) form ln(ψ + 1.6) 8 alone explained 58% and 74% of the variation in the logarithmic S B parameters lnλ and lnδ, respectively.The shape of the dbh distributions varied from inverse J-shaped distributions (ψ ≤ 0.3) to peaked and positively skewed distributions (0.9 ≤ y < 1.0).The final estimated models are presented in Appendix 1.An example of the dynamics of the dbh distribution of a certain stand, including the observed and predicted distributions using the MSUR model, is given in Fig. 3.The total stem numbers varied between 4700 and 1010 ha -1 between the first and last measurements.

Model Performance in the Study Data
When generating the modeling data with predicted dbh distributions, each model provided good results for the generated number of stems as well as for the 'volume' (Sd 3 ) and 'value' (Sd 4 ) of the stock.The differences between the models were relatively marginal.In general, more sophisticated regression models improved the model behaviour compared with the OLS model (Table 3).The clearest improvement was in a decreased error variation, particularly in the model precision, s b .MSUR provided the most precise alternative according to the error deviation in Sd 3 and Sd 4 , but the bias in stem number (N) was the highest, although still less than 1%.The MIX model resulted in only slightly biased stand characteristics, but the precision was among the lowest.In the study data, the performance of both MIX and MSUR improved when the estimated random stand effects were utilized.The SUR model generally provided good results.

Validation of the Models in the Independent Data Sets
According to the results in the HARKO test data set, all the models seemed to predict the stem number, Sd 3 ('volume') and Sd 4 (the 'value' of the growing stock) relatively similarly as regards the biases and their standard deviations.Depending on the model, the average bias in N varied from 2.7% to 5.5% overestimation, the bias in Sd 4 from 0.6% to 1.7% overestimation, and that of Sd 3 was practically unbiased (Table 4).Even though each model provided good results, the SUR and MSUR models provided the best estimates and some of the smallest error deviations.
The MIX model provided precise estimates for stand 'volume', but overall the highest biases and clearly the worst precision for stem number.Finally, when the errors of the compared models were tested against the thinning intensity, no visible trends were found.For example, with MSUR the bias in N varied between only 3.3% to 3.9% overestimation.This is also evident from the small The biases in N were -8.5, -3.6, -2.9, and -2.8 %, respectively.
error deviation (s b of 1.7%).The MIX model was an exception in that showed a trend from 6.7% overestimation in unthinned stands to 3.8% overestimation in heavily thinned stands, thus also resulting in a greater error deviation, s b of 2.6%.However, in spite of using the temperature sum as an independent variable, the models showed a slight trend in the bias of N with respect to the annual temperature sum (DDY).For example, if DDY was less than 800 degree days (i.e.below the values of the modeling data) then the bias in N using MSUR was -5.0%, while with DDY greater than 1000 degree days it was -1.3%.The Weibull distribution of Mykkänen (1986) performed well in 'volume' and 'value' of the growing stock, while N was imprecise and underestimated by 9%.In contrast, the percentile based model of Kangas and Maltamo (2000b) provided excellent N estimates while the volume and the value estimates were more biased (0.8-1.7% overestimates) and slightly less precise.Both test models showed slight trends in the bias of N with respect to thinning intensity; unthinned stands were underestimated by 2.1% and 14.6% of N when using the models of Kangas and Maltamo (2000b) and Mykkänen (1986), respectively, while the estimates of N in heavily thinned stands were almost unbiased (-0.1% and 1.0%, respectively).The Weibull distribution model of Mykkänen (1986) resulted in a slightly increasing underestimation in N with respect to increasing DDY.The percentile based model (Kangas and Maltamo 2000b) did not show such a trend.
In the test data of the upland sites (INKA) the models did not behave in the same way as in the HARKO data, which were collected from drained peatlands.First of all, the higher variation in stand characteristics resulted in a significantly higher RMSE and almost doubled the relative error deviation (even four-fold in N) compared with the HARKO data (Table 4).One reason for the higher variation in the INKA data may be that the HARKO data originate from thinning experiments, which have been established in stands where silvicultural thinnings were appropriate to carry out, whereas INKA includes a larger range of stand densities and stockings.In this test data, however, SUR and MSUR again performed well.Typically, the lowest test values were found with SUR and the second best with MSUR.An exception to the HARKO data was the precision in N, in which MSUR was superior (4.8%) to the second best SUR (7.1%).
The earlier models that were included gave relatively similar results as before.The model of the Weibull distribution (Mykkänen 1986) constructed for mineral soil sites did not predict the stem number accurately (19% underestimate, 29% RMSE), but still resulted in satisfactory accuracy in the 'volume' and 'value' characteristics, i.e.Sd 3 and Sd 4 (Table 4).The percentile based model (Kangas and Maltamo 2000b) gave excellent accuracy in stem number (only 3% underestimate and 5% RMSE) and also good predictions for 'volume' and 'value' of the growing stock.
When the earlier test values were ranked and the sums of the rank scores calculated, we found that SUR and MSUR gave the lowest rank in modeling data (20), while OLS received the highest ranking, namely a score of 27 (Table 5).In the HARKO data set the order of model performance from first to last was MSUR, SUR, MIX, and OLS, while in the INKA data set the order was SUR, MSUR, MIX, and OLS.When the overall ranking was calculated across all the data sets, the result in rank order was 1. SUR, 2. MSUR, 3. MIX, and 4. OLS.It was notable that both of the test data sets gave a relatively similar ranking to that of the modeling data.Nevertheless, SUR clearly performed better in INKA than in the modeling data.In conclusion, all the advanced models received smaller rank score (50-77) than the OLS model ( 85).Note that the mixed-effects models (MIX, MSUR) were applied as a fixed model when predicting the S B parameters.Thus, improvement of the model predictions on the basis of the estimated random effects for each stand in the modeling data was not used here.

Discussion
The parameter prediction method has been criticized on the basis of the typically low correlations between distribution parameters and stand characteristics (see Knoebel and Burkhart 1991).In the Nordic countries, Mønnes (1982), Tham (1988) and Holte (1993) have used the S B distribution, and predicted the dbh-frequency distributions by means of the 'biologically sounder' percentile method.However, Siipilehto (1999Siipilehto ( , 2000) ) applied direct parameter prediction in constructing S B distribution models for Norway spruce, Scots pine and silver birch on mineral soil sites.The typically low degree of determination was overcome by using N as an additional predictor with the variables d gM and G in the form of the shape index.Different transformations of the shape index as well as other proportions (e.g.ln(N/G) or ln(G/g M )) were tested as explanatory variables.Some combinations gave equal or even slightly better fit but the presented models were found the most stable for prediction purposes.It is noteworthy that there were systematic differences in distribution parameters between site types (see also Sarkkola et al. 2005).The effect of site type was most obvious on parameter δ resulting in more peaked distributions on genuine forested sites than on sparsely forested composite sites.Similar differences could be seen with respect to nutrient regime.However, the effect of site type did not become statistically significant in the models.The most probable explanation is that much of the variation in stand structure due to site properties was implicitly accounted for by the other explanatory variables in the models; particularly the shape index explained most of the variation in SB-parameters, and it was firmly correlated with site properties.
The shape index showed a slightly wider variation in the modeling data of peatland sites (0.17-0.99) than that of the mineral soil sites (INKA), namely 0.20-0.99.Note that the shape index should not exceed one (1.05 was found critical threshold).If so, the shape of the distribution became extremely peaked and the number of stems became simultaneously more and more overestimated.In the present study, the shape index proved to be extremely useful, but also the general stand variables (d gM , G and N) could explain much of the variation in S B parameters, especially in λ.A recently developed alternative approach (Cao 2004), in which the stages of fitting and modeling the distribution are combined, was not utilized here.
The advantage of using the basal area-dbh distribution instead of the dbh-frequency distribution is primarily due to the common inventory practice (basal area median dbh and stand basal area are the main characteristics collected in the field).Secondly, using the basal area distributions provides more reliable total volume predictions (e.g.Gobakken and Naesset 2004).For example, when Kangas and Maltamo (2000a) applied dbh-frequency distributions for pine, they found a RMSE of 4.4% and 6.8% in the basal area, and a slightly higher RMSE of 5.2% and 8.5% in the volume when using the distribution free percentile method and the Weibull distribution, respectively.When they calibrated these distributions against the real basal area, the accuracy was highly improved.On the other hand, the percentile models that were constructed directly for the basal area-dbh distribution resulted in a RMSE in the volume of either 1.5% or 2.3%, irrespective of whether N was used as a predictor or not (Kangas and Maltamo 2000b).
When Kangas and Maltamo (2000c) investigated the model performance in different data sets, they found slightly better results in the INKA data set than we found in the present study.For example, the bias in N was only 2.1% compared with 3.15%, in our study.This might be due to certain restrictions that they applied for the ranges in stand variables in the INKA test data.We noticed relatively high overestimations in N in the HARKO data (2.7-5.5%) when using the S B distribution.This was partly due to the truncated dbh distributions (only dbh ≥ 5 cm were measured).If the predicted dbh distribution is truncated similarly, then the overestimation in N diminishes to 1-3 percent.The percentile based models gave better accuracy in N, probably because the minimum dbh (d 0 ) was modeled (Kangas and Maltamo 2000b).
Due to the large variation in the modeling data with respect to the stand developmental stages, and therefore large variation in the individual stand characteristics, the parameter estimates given by the different estimation methods were close to each other.Each model provided good results when the alternative regression approaches were compared with each other in the different data sets.However, some differences could still be found.Taking into account the hierarchical structure of the modeling data (MIX), cross-model correlations (SUR) or both (MSUR) improved the model predictability compared with the OLS model.In the INKA test data set this improvement proved to be the most obvious.SUR and MSUR provided the overall best approach when the stand structure was focused in terms of stem number, volume (Sd 3 ) and value (Sd 4 ) of the growing stock.
The advantage of SUR and MSUR was related to the high cross-model error correlation (0.74), which could be utilized in re-estimating the models.The variance components in MIX only slightly improved the accuracy of the model compared to OLS.One explanation for this may be the relatively low correlation between repeated observations within the stands.In fact, the values of the S B parameters may vary considerably between observations, even without any management operations between measurement occasions.This may be due to tree mortality and ingrowth, but also to variations in tree growth rate among tree size classes.The between-stand variation also differed only slightly from the residual variation (see appendix).Finally, applying mixed-effects models as fixed models when predicting the S B parameters in a FORTRAN application, means that the estimated random effects for each stand in the modeling data were not utilized.The random effect is unknown when predicting the distribution for new stands, and therefore we decided not to use it in this study.
In conclusion, SUR and MSUR provided the best estimates for stand characteristics on different site types.MSUR was statistically more relevant because the hierarchy of the modeling data was taken into account.The ranking of the models was logical, but may not be widely generalized.It is highly dependent on the structure of the data and that of the models.If the data were unbalanced, e.g. with respect to the age of the stands (or years elapsed since drainage), then the variance component (MIX, MSUR) model and the OLS model would have resulted in considerably different estimated parameters.In such a case, the variance component model could obviously provide a more reliable result than OLS (i.e.relevant within stand development).In any case, the significance of the estimated parameters is more reliable in mixed-effects models because the dependencies among repeated measurements are better accounted for.
The models that performed the best (SUR, MSUR) in this study are the most reliable parametric dbh distribution models so far constructed for Scots pine stands in Finnish conditions, at least as regards the variation in stem number.The models showed reliable performance even in the independent test data sets, which included a wide range of thinning intensity or site properties.In Norway, Gobakken and Naesset (2005) presented models for lidar-based estimation of basal area dbh distribution.They found percentile prediction of a two parameter Weibull function to perform better than a system of ten percentiles.In Finland, S B model for birch (Siipilehto 1999) has performed well but the distribution free percentile based models using N as an additional predictor (Kangas and Maltamo, 2000b) have proved to be generally the best of the numerous earlier distribution models (see Kangas andMaltamo 2000c, Maltamo et al. 2002).In the present study, the S B distribution gave at least comparable reliability.The differences in generated stand volumes were relatively negligible, but they were more obvious in the stem number.Flexible changes in the dbh distribution enabled accurate estimation of the stem number.This was especially evident in the predicted proportion of trees with a small diameter, which was not so evident for the stand volume estimation, but that can be regarded as important for the simulation of the stand dynamics.
), where the hierarchy causes correlations between observations (var(e) = Z var(c)Z' + var(u)) and the fixed parameters should be estimated using GLS instead of OLS(Gregoire et al. 1995, McCulloch andSearle 2001).The Maximum Likelihood (ML) and Restricted Maximum Likelihood (REML) estimation methods lead to GLS estimation of the fixed parameters after first estimating the variance components(McGulloch and Searle  2001, p. 308).

Fig 1 .
Fig 1.The variation in shape index (ψ) with respect to the basal-area median diameter (d gM ) in pine-dominated drained peatland stands.

Fig 2 .
Fig 2. The relationship between the shape index (ψ) and the estimated S B parameter λ and δ (left), as well as the relationships of the corresponding linear logarithmic transformations used in the models (right).

Table 2 .
Mean stand characteristics and derived shape index (ψ) in the HARKO test data set on peatland sites and the INKA test data set on mineral soils.

Table 3 .
Relative bias, RMSE and standard deviation of the relative error, s b (%) in total stem number N, in Sd 3 and in Sd 4 using the modeling data.The lowest values given by the models are indicated in bold and the highest in italics.

Table 4 .
Mykkänen (1986)stics in the independent HARKO and INKA data sets.Test model K & M was the percentile based dbh distribution model ofKangas and Maltamo (2000b);Mykkänen (1986)applied the Weibull distribution.The lowest values given by the models are indicated in bold and the highest in italics.

Table 5 .
Ranking of the compared models.Rank sums were calculated from the ordered test values (bias, RMSE and s b ) by giving the lowest test values rank one, the second lowest rank two and so on.