Local Prediction of Stand Structure Using Linear Prediction Theory in Scots Pine-Dominated Stands in Finland

This study produced a family of models for eight standard stand characteristics, frequency and basal area-based diameter distributions, and a height curve for stands in Finland dominated by Scots pine (Pinus sylvestris L.). The data consisted of 752 National Forest Inventorybased sample plots, measured three times between 1976 and 2001. Of the data, 75% were randomly selected for modelling and 25% left out for model evaluation. Base prediction models were constructed as functions of stand age, location and site providing strongly average expectations. These expectations were then calibrated with the known stand variables using linear prediction theory when estimating the best linear unbiased predictor (BLUP). Three stand variables, typically assessed in Finnish forest management planning fieldwork, were quite effective for calibrating the expectation for the unknown variable. In the case of optional distributions, it was essential to choose the weighting of the diameter distribution model such that the available input variables and the model applied were based on the same scale (e.g. arithmetic stand variables for frequency distribution). Additional input variables generally improved the accuracy of the validated characteristics, but the improvements in the predicted distributions were most noteworthy when the arithmetic mean and basal area-weighted median were simultaneously included in the BLUP estimation. The BLUP method provided a flexible approach for characterising relationships among stand variables, alternative size distributions and the height–diameter curve. Models are intended for practical use in the MOTTI simulator.


Introduction
In Finland, the description of a stand is commonly simplified to visually assessed mean and sum characteristics in forest management planning (FMP) fieldwork.Stand characteristics are assessed by tree species (see Solmun… 1997).The number of stems is typically assessed in young stands, and the basal area in advanced stands.However, the stem number is sometimes required additionally for basal area and basal area-weighted variables (Kuvioittainen… 1998).Alternative choices related to the stand characteristics assessed may cause problems in the use of FMP inventory data as input variables in simulation systems.Finnish simulators such as MELA (see DemoMELA, Siitonen et al. 1996), MOTTI (see MOTTI software, Hynynen et al. 2005), and MONSU (see MONSU, Pukkala 2004) are based on tree-level data while the SIMO (Kalliovirta 2006) simulator incorporates both stand-level and tree-level simulation options.In any case, diameter distributions are needed for the calculation of assortment volumes (Kalliovirta 2006).Thus, tree size distribution models are needed for converting stand-level information into tree-level information.
Most of the diameter distribution models in Finland are based on basal area and basal areaweighted medians (e.g.Päivinen 1980, Mykkänen 1986, Maltamo 1997, Kilkki et al. 1989).More recently developed Johnson's SB and percentilebased distribution models attempt to enhance the accuracy in predicted distributions by using the stem number as an additional independent variable (Siipilehto 1999, Kangas and Maltamo 2000a, Siipilehto et al. 2007).Stem frequencybased Weibull models have been presented for drained peatlands by Sarkkola et al. (2003Sarkkola et al. ( , 2005)).Models for young stands based on Weibull height distribution (Valkonen 1997, Siipilehto 2009) and the models for tree height in Finland (Veltheim 1987, Siipilehto 1999, Mehtätalo 2004, 2005, Eerikäinen 2009) require different sets of input variables.Therefore, we may still have an incomplete set of input variables for some present models and/or unsuitable models for the data in use.
Thus, the need for modelling individual stand characteristics or relationships between stand characteristics arises from the difference between known and needed input variables for applied models, but also because of the changes and alternatives in FMP or National Forest Inventory practices (e.g.Nissinen 2002, Nuutinen 1986, Eid 2001).At present, the utilisation of satellite images and laser scanning data may create pressure for changes -for example, by including the number of stems in the variables assessed, because of the increased accuracy (Packalén and Maltamo 2007, Maltamo et al. 2007, Vohland et al. 2007) in comparison with traditional FMP fieldwork (e.g.Haara andKorhonen 2004, Kangas et al. 2004).
To overcome the problem of missing input variables, Siipilehto (2006) presented a linear prediction application for eight stand characteristics of Norway spruce stands.Basic models for spruce were function of stand age, location, site and origin.The expectations for unknown stand characteristics were then calibrated with the known stand variables using linear prediction theory (see Lappi 1993Lappi , 2006)).For example, the original RMSE of 34% for basal area was reduced to 18% when the mean diameter and stem number were used for calibration (Siipilehto 2006).This work is a direct follow-up and extension to that study.Extension consisted of simultaneous prediction models for three optional diameter distributions and a height curve in order to provide all the necessary tools for stand structure in the same package.
The purpose of this study is to develop a family of models providing tools to predict the structure of Scots pine stands in Finland.Simultaneously to stand characteristics, models for unknown parameters of the diameter distributions and heightdiameter relationship are estimated.The models are intended to provide 1) availability of prediction with a limited number of known characteristics as a basic model; 2) calibration of the prediction with an arbitrary set of known stand variables based on linear prediction theory; 3) a theoretically consistent framework for comparing alternative size distribution models; and 4) a broad range of applicability ranging from sapling stands to mature stands.Tree diameter is defined as diameter-at-breast-height (dbh) and its measurement is assumed as being error-free.In this study, we focused on the extension to Siipilehto (2006), i.e. a prediction of stand structure in terms of dbh distribution and the height-dbh curve.
The distributions studied are Weibull for a dbhfrequency distribution (W N ) and basal area-dbh distribution (W G ) and Johnson's SB distribution for basal area-dbh distribution (SB G ). Height was modelled with a function by Näslund (1936).Models are intended for practical use in a MOTTI simulator.

Modelling and Test Data
The data was obtained from permanent sample plots, established in seedling stands between 1984-1989 (TINKA: H dom < 5 m when established) and more advanced stands (INKA: H dom ≥ 5 m) between 1976-1983.The stands are based on a subsample of the Seventh National Forest Inventory.Each standwise sample plot consisted of a cluster of three circular plots within a stand.The radius was selected such that the total number of trees tallied was about 100-120 per stand.A smaller radius has been applied within each circular plot to select one-third of the sampled area for height and other more detailed measurements (Gustavsen et al. 1988).For the purpose of the present study, the cluster of three largerradius sample plots were combined for reliable stand characteristics and especially for reliable distributions (see Shiver, 1988) representing the whole stand.Expected heights were given from Näslund's height curve for trees with missing tree height in INKA data set, whereas in sapling stand TINKA data, all crop trees have been measured for height and dbh.In the INKA data set, the smallest trees (dbh < 5 cm) were not measured if they were considered unsuitable for growing (e.g.having inadequate growing space).
The number of Scots pine-dominated INKA and TINKA stands was 562 and 190 respectively.In this study, data was restricted to a mean diameter (D) greater than 1.5 cm.This was done in order to avoid the anomalies in stand characteristics resulting from the low proportion of trees above breast height in the youngest state of stands.The remeasurements were carried out twice between 1981-2001, typically five and ten years after establishment.Because of the repeated measurements, the total number of observations in the study was 1858.Data was randomly divided into a modelling data set (75%) for estimation of the models and a validation data set (25%) for evaluating the estimated models.The final number of observations was 1392 for modelling (see Table 1) and 467 for validation (see Table 2).The modelling data covered dominant heights of 2.5 to 28.6 m (as seen in Table 1).
In current forest management planning fieldwork, the typical set of variables characterising stand structure (i.e.dbh and height distributions) is either 1) the number of trees per hectare, the mean diameter, and the mean height for young stands; or 2) the basal area per hectare, the basal area median dbh, and the corresponding height for advanced stands (see Solmun… 1997, Kuvioit-tainen… 1998).Hence, the randomly selected test data was divided into young and advanced stands due to the assumption of different combi-  1) Dbh distribution and height-dbh models are required for the tree-level information of stand structure.Selected functions and fitting methods for the parameters are given in section 2.3.2) The basic prediction models for the obtained parameters are constructed simultaneously with the models for stand characteristics in section 2.4.The estimated parameters of the basic models are given in the results, Tables 3 and 4.These basic models are intended to provide logical but strongly averaging age-dependent predictions for different sites and locations of stands.3) In practice, the expectations and error variances for unknown parameters and unknown stand characteristics are calibrated with the set of known stand characteristics in order to get the final, improved prediction for the stand structure.This method is called linear prediction and it is presented in section 2.5.The idea is that the predictors (calibrating variables) are not fixed beforehand; instead, arbitrary set of included stand characteristics can be used for calibration.Linear prediction application utilises a cross-model error variance-covariance matrix for calibration (i.e.estimating the best linear unbiased predictor).The estimated matrix (along with the correlation coefficients) is given in two parts in the results, in Tables 5 and 6. 4) Given examples and the final model evaluation (Tables 7-10) looks for the efficient stand characteristics for calibrating the unknown parameters of the included dbh distributions and height curve.

Modelling of Tree Dimensions
The diameter distribution was characterised by Weibull and Johnson's SB distribution functions.The two-parameter Weibull function for tree diameter d has the following probability density function (PDF) (e.g.Bailey and Dell 1973): where b > 0 is the scale and c > 0 the shape parameter.According to size-biased theory, Gove and Patil (1998) showed that weighting the initial two-parameter Weibull frequency distribution with tree basal area leads to a standard gamma distribution, where the parameter k = (1 + 2 / c) of the gamma function Γ(k).Similarly, if the fitted basal area-dbh Weibull distribution is returned into the dbh-frequency distribution, it leads to a gamma distribution with parameter k = (1 -2 / c).These features were convenient when the samples were taken either from frequency distribution representing young stands or from basal area distributions representing advanced stands in the model evaluation.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.One of them, SB distribution, has been used for size distribution in forestry (e.g.Hafley andSchreuder 1977, Mønness 1982).It has the following PDF (Eq.2) which applies the bounded transformation function shown in Eq. 3.
where γ and δ are shape parameters, ξ and λ are location and scale parameters, and d is the diameter observed in a stand plot.The height-diameter relationship was characterised through the use of the height curve (Eq.4) by Näslund (1936) fitted in the linearised form (Eq. 5).
where β 0 and β 1 are the estimated parameters.
The residual e z from Eq. 5 was assumed to have homogenous and normally distributed standard error (see Näslund 1936, p. 53).Distribution functions were fitted to each stand separately through the maximum likelihood (ML) method (see Bailey andDell 1973, Schreuder andHafley 1977).The Weibull dbh-frequency distribution was fitted with the NLIN procedure in the SAS software application (see the appendix of Cao 2004) and both of the basal area-dbh distributions through the use of the author's Fortran programs as in the work of Siipilehto (1999).A three-parameter form of the SB distribution was applied by fixing the minimum ξ = 0.In order to avoid huge maximum endpoints (outliers in the modelling point of view), the parameter λ was restricted to be two times the maximum observed dbh at most.Näslund's height curve was fitted to each stand separately by tree species in SAS by means of the NLIN procedure, convenient for saving the estimated stand-specific parameters and the standard errors to the output data.The error term was used to generate random variation around expected tree heights for those tallied trees missing a height in the test data set.For more details on the error structure see Siipilehto (2000).The estimated parameters (see Table 1) were regarded as true values for stands, when the models for these parameters were formulated against stand age and site factors.

Construction of the Basic Models
The structure of the basic models was assumed to be multiplicative as with previous models by Nuutinen (1986), Eid (2001) and Siipilehto (2006) for stand characteristics and models by e.g.Kilkki et al. (1989), Mykkänen (1986) and Siipilehto et al. (2007) for distribution parameters.The models were fitted by means of multiple linear regression after logarithmic transformations in order to linearise the models and homogenise the residuals.The following 17 dependent variables were fitted simultaneously: 1) total basal area (G, m2 ha -1 ), 2) total stem number (N, ha -1 ), 3) arithmetic mean diameter (D, cm), 4) basal area median diameter (d gM , cm), 5) dominant diameter (D dom , cm), 6) mean height (H, m), 7) basal area median height (h gM , m), 8) dominant height (H dom , m), 9) parameter β 0 and 10) parameter β 1 of Näslund's height curve, 11) parameters b and 12) c of the Weibull dbh-frequency distribution, 13) parameters b G and 14) c G of the Weibull basal area-dbh distribution, 15) parameters λ 16) and δ of the SB basal area-dbh distribution, and 17) parameter γ from the same.The basic models represented the average development of each dependent variable over the stand total (biological) age.The common structure of the candidate standwise model was as follows: ln(Y + m) = a 0 + a 1 (T+l) k + a 2 ln(DDY) + a 3 Origin + a i S i + a j Thinn j + ε (6 where m = 1 in the model for D, d gM , and D dom , and for parameter b of the Weibull function and λ for the SB function, and otherwise m = 0; T = total age (in years); the candidate power k is -1,-0.9,…,-0.5;constant l is either 0 or 10; DDY = degree days (i.e. the average annual sum of the temperature on days with a mean temperatures above 5 °C); Origin is a dummy variable (with a value of either 0 or 1) for artificial regeneration methods; S consists of dummy variables associated with a certain site (i) defined as forest types by Cajander (1925), and supplementary site characteristics such as stoniness and paludification; Thinn consists of dummy variables associated with a thinning treatment (j) being either intermediate thinning or heavy thinning preparing for natural regeneration; a 0 to a j , are estimated parameters of the model; and ε is the random error.Dominant trees represented the hundred thickest trees per hectare.The models for the sum characteristics were formulated for the whole stand as a potential of the site (total stem number and total basal area).The proportion P (0 < P ≤ 1) of Scots pine was given adding ln(P pine ) to the model with a fixed (restricted) parameter of 1.0.Thus, if the stand was a pure pine stand, the effect was zero.The inverse of the parameter c was found to be the best transformation for the Weibull frequency distribution.The inverse transformation can be motivated by the known percentile estimator (Dubey 1967) and previous studies by Mabvirura et al. (2002) and Robinson (2004).However, a model for ln(c G ) was preferred to the inverse of c G for the basal area-dbh distribution because of its slightly superior behaviour in advanced stands.The models were multivariate; i.e. several models were fitted simultaneously (systems of equations) to the same data set.As a result of the correlation between the errors of different models, fitting the individual models separately would not be an efficient approach to estimation.Instead, the error structure of the ordinary least square (OLS) fit was utilised in the re-estimation of parameters of correlated response variables using the seemingly unrelated regression (SUR) method (Zellner 1962).Estimation was performed with the SYSLIN procedure in SAS (see SAS/ETS User's Guide 1993), giving the estimated parameters of the first step OLS models, the re-estimated parameters of SUR models and the cross-model error variance-covariance matrix (and its transpose and inverse matrix) as an output.

Linear Prediction Application and Model Validation
The error terms of statistical models are random variables.In this study, the error term consisted of the residual term only.The cross-model error variance-covariance matrix is valuable when calibrating the expected value (µ 1 ) with known variables x 2 using linear prediction theory (e.g.Lappi 1991Lappi , 1993)).The best linear unbiased predictor (BLUP) for variable x 1 is where x 1 is a scalar (dependent unknown variable) and x 2 is a vector (known stand variables), σ 12 is a row vector (covariances between unknown dependent and known variables), and Σ 22 is the variance-covariance matrix of x 2 (between known variables).The variance of the prediction error after calibrating the dependent variable x 1 is: where σ 11 is a scalar (initial variance of the residual error of the dependent variable), and σ 12 ´ is a transpose of the row vector σ 12.When assuming the residual errors of the logarithmic models to be multinormally distributed, half of the error variance (s ε rithmic model back to the original scale.With the inverse transformation (1 / c), the bias correction term is s x ε 2 1 2 / ˆ (see Lappi 1993, p. 91-93).Thus, in the applied method, the variance had to be recalculated whenever the vector of known stand variables (x 2 ) for prediction was changing.An example of the BLUP estimation for Näslund's parameters is given in the appendix.
In development and validation of the structure of the candidate models, the RMSE for the fitted model was examined, as was the residual variation of each model with respect to its expected value, temperature sum and stand age.In addition, models were checked visually to ensure logical behaviour of the models when one is using a wide range of stand ages.The effect of known stand characteristics on BLUP estimation was examined with respect to RMSE -i.e. the square root of the variance shown in Eq. 8 for the final models.Theoretically, any combination of the stand variables presented could be used for linear prediction.However, the linear prediction application was briefly checked, predicting basal area-weighted characteristics (d gM , h gM , and G) with arithmetic stand characteristics (D, H, and N), and vice versa as in Siipilehto (2006).As for the predictability of the size distributions and height curve, additional stand characteristics were tested as calibrating variables in order to study their ability to improve reliability in the unknown parameters.Some of the combinations are more realistic than others, because the dominant height is sometimes included in the characteristics assessed for FMP, and stem number is required additionally for basal area and basal area-weighted variables (Kuvioit-tainen… 1998).
The final validation for the family of models generating individual trees (i.e. the dbh distribution and height curve) was performed by means of bias and RMSE for sum and dominant tree characteristics as well as for total, saw timber, pulpwood and waste wood volume.Stem total and assortment volumes were calculated by means of volume and taper curve functions by Laasasenaho (1982).Waste wood is potential energy wood and therefore an interesting characteristic at present (e.g.Tahvanainen et al. 2007, Maltamo et al. 2007).In addition, accuracy in terms of waste wood indicates the accuracy of the left tail of the predicted distribution.On the other hand, accuracy in dominant tree characteristics was examined in order to yield an idea of the reliability of the height curve and the right tail of the distribution.These results were calculated from systematic samples of 40 trees taken from the predicted distribution.
For effective utilisation of the data, young sapling stands were also included in the data.This means that in the case of a mean height of, let us say, less than 4 m, a considerable proportion of trees can be less than 1.3 m high.Thus, an auxiliary logistic model ( 9) was estimated in order to describe the probability of a tree being above breast height: where a 0 = -3.530and a 1 = 2.3241.The model was fitted in SAS with the NLMIXED procedure, which fits non-linear mixed models by maximising an approximation to the likelihood integrated over random (stand) effects.The estimated proportion was used to correct the total number of stems to apply to trees above breast height when, for example, basal area and volume (volume of trees h < 1.3 m was ignored) were calculated.Without this correction, they could be severely overestimated in sapling stands.Indeed, according to the estimated model ( 9), 23%, 76%, 97% and 99.7% of Scots pine stems were higher than 1.3 m when the mean height H was 1, 2, 3 and 4 m respectively.In advanced stands, distributions were scaled according to the basal area instead of stem number.The effect of each stand variable on the vector of unknown dependent variables was focused by including them one by one as a known variable in the family of BLUP models.Some examples are given to illustrate the model application.

Estimated Models
The SUR method proved to be advantageous, by improving the significance of some of the independent variables from the original OLS fit, Using mainly the same variables and transforma-  -0.6 ) and height characteristics (T -0.5 ) helped to achieve a reasonably small dbh value for the short trees (i.e.logical behaviour just above breast height).Degrees of determination were quite high for dbh and height characteristics (76-84%) but also considerable for basal area and stem number, at 66 and 53% (Table 3).
The most typical site for Scots pine was Vaccinium-type (VT) sub-xeric heath (Cajander 1925), representing the reference level of the models.The effects of site fertility on stand variables were characterised using dummy variables.For example, on Oxalis-Myrtillus-type grove-like sites (OMT sites), the heights (H dom , h gM , and H) were about 28-40% greater than with VT, while on Calluna-type xeric heath (CT) these heights were about 15% and for Cladonia-type (ClT) about 40% lower (Table 3).The effect of site fertility on the stem number and basal area was obvious and could be explained by the differences in growth potential with respect to site classes.On the other hand, stoniness and paludification reduced stand mean variables by around 3% and 14%, respectively.Artificial regeneration showed a significant effect on stand characteristics and the effect of planting was greater than that of sowing.Intermediate thinning (Thinn) was significant for sum and height characteristics, while heavier thinning in preparation for natural regeneration (Thinn R ) was significant for each stand characteristic (see Table 3).
Site factors, stand origin and past thinnings had an effect on the height curve as well as on dbh distributions.Note that scale parameter b of the Weibull functions and the maximum λ of the SB function had the same model structure as dbh characteristics because of the close correlation between them (see Table 4).The logical dependencies of D < b, d gM < b G , and D dom < λ were fulfilled.Degrees of determination were considerable for the parameters modelled (32-78%), except for the shape parameters of the SB distribution (9-13%).Note that the adjacent dummy variables have been combined where they have equal estimated parameters in Table 4 (e.g.OMT and MT in the model for ln(β 0 )).
Covariances and correlations of the crossmodel residual errors of the stand characteristics are given in Table 5.It was not surprising that correlation between many stand characteristics was very strong.A correlation coefficient above 0.7 was found in 14 out of 28 cases (see Table 5).The highest correlation, 0.96, was between h gM and H dom .Covariances and correlations between the residuals of the modelled parameters of the size distributions and height curve and stand characteristics are shown in Table 6 and Fig. 1.Analysis revealed relatively high absolute values of correlation coefficients between the height curve parameters and some stand characteristics (e.g.r = -0.4 between β 0 and G and H, and

The Effect of BLUP Estimation on the RMSE
BLUP is a flexible approach -any of the known stand variables can be used for predicting (cali-brating the expected value of the basic model using Eq.7) the unknown dependent variable (see Lappi 1993Lappi , 2006)).Similar to Siipilehto (2006), stand age was the only stand variable required in the basic models.The prediction efficiency of the individual stand variables in terms of RMSE (i.e. the square root of Eq. 8) has been given for stand characteristics for Norway spruce (see Table 4 in Siipilehto 2006).Models for Scots pine behaved in astonishingly similar ways.For example, calibrating d gM with N, D and H resulted in RMSE of 10.8 and 10.5% for spruce and pine, respectively.However, calibration with a sum and a mean dbh characteristic decreased the original RMSEs of 43.2 and 44.5% for G and N to only 13.6% and 20.3% respectively.Thus, the calibrated sum characteristics for pine proved more accurate than those for spruce (18.0 and 24.3%) even though the original RMSEs for pine were higher.Nevertheless, due to overall similarity, the other detailed results between stand characteristics are not worth giving here, but they can be easily calculated using the Eq. 8 and the covariance matrix in Table 5.
The effect of a stand characteristic in improving the precision of the size distribution and height curve parameters varied clearly between the characteristics utilised.It was obvious that knowing both one dbh and one height characteristic (D and H or d gM and h gM ) was relatively efficient for calibrating the parameters of the height curve in terms of reduction in RMSE (Fig. 2: A, C).In addition, known D dom and H dom increased the accuracy in β 0 and β 1 but the additional sum characteristics, either N or G, had no effect.Calibration of the parameters of the Weibull and SB distributions was relatively efficient with D or d gM .Adding mean height and sum characteristics (H and N or h gM and G) affected the RMSE for size distribution parameters only marginally.However, adding dominant tree characteristics clearly reduced the RMSE for each parameter.As for the sum characteristics, additional knowledge of the basal area did not improve the accuracy (Fig. 2: A, B), but the additional stem number was relatively efficient in reducing the RMSE for each size distribution parameter (Fig. 2: C, D).The most efficient calibration was achieved with the additional mean diameter.Thus, including both D and d gM in the BLUP models reduced the RMSE for each parameter as shown in the bar furthest to the right in Fig. 2.
Fig. 2 proves that the dbh-frequency distribution could be more efficiently calibrated with stand variables that are directly associated with frequency distribution (i.e.D, H, and N).Stand variables that are typically assessed by means of relascope sample plots (i.e.d gM , h gM , and G) were surprisingly inefficient for calibrating the shape parameter of the Weibull distribution.Knowledge of the additional stem number seemed quite essential for the shape parameter of W N and W G but simultaneously also enhanced b of W N (see Fig. 2: C, D).Mean and dominant diameter and stem number were significant additional variables for calibration of the models for the SB G parameters, especially when we were dealing with the stand characteristics on the arithmetic scale (Fig. 2: B).On the other hand, when starting with the basal area-based characteristics (Fig. 2: D),

Model Behaviour with Respect to Stand Density
The following example illustrates the effect of stand density on the stand characteristics and dbh distribution.Let three given hypothetical stands represent a 25-year-old Scots pine stand sown on a sub-xeric (VT) site in central Finland (annual temperature sum of 1200 °C), and assume that the density of these Scots pine stands is 2000 ha -1 , 4000 ha -1 , and 8000 ha -1 , respectively.Simultaneously with increased density, the shapes of the distributions shifted from more or less symmetrical bell shapes to being more and more skewed to the right as shown in Fig. 3.In addition, dbh distributions moved to the left, resulting in mean diameters of 8.0, 6.5, and 5.3 cm, respectively.The above mean diameters from the predicted Weibull distributions followed the corresponding density-calibrated BLUP estimates for D -namely, 7.8, 6.3, and 5.1 cm.Also, the calculated basal areas (11.7, 15.8, and 21.3 m 2 ha -1 ) followed quite a similar pattern to that of the density-calibrated BLUP estimates for G (11.2, 15.3, and 20.9 m 2 ha -1 ).Similar comparisons were also made for 25-year-old naturally regenerated stands and 20-year-old planted stands.According to these comparisons, the BLUP estimation for D and G differed only slightly from the respective characteristic calculated from the calibrated Weibull distribution.

Calibration of the Dbh Distribution to a Real Stand
The following example shows the calibration effect of the known stand characteristics on the predicted dbh distribution against one existing stands.Size distributions of the basic model represent the average for a particular site, stand age and origin.Thus, because of natural variation and varying management practices, the observed structure and the expected structure can be considerably different.The example, plot no.1017, represents a 21-year-old planted stand on an MT site in southern Finland (DDY: 1245 °C).Prediction for this stand was performed with the W N model.The observed N, D, and H were 950 ha -1 , 11.0 cm, and 7.1 m, respectively.The respective expected values of the basic models were as follows: 1224 ha -1 , 9.1 cm, and 7.2 m.The additional stand characteristics checked for calibration were G, d gM , D dom and H dom , which had values of 9.6 m 2 ha -1 , 12.4 cm, 15.7 cm and 7.7 m.The respective expected values were as follows: 8.6 m 2 ha -1 , 11.8 cm, 15.9 cm and 9.6 m.The expected W N distribution for the basic models (W(0)) was to the left of the observed distribution, because of the lower expected D (see Fig. 4).In the second step, prediction was conducted with N, D, and H as common known variables for young stands.The resulting dbh distribution was situated comfortably over the observed distribution, but the tail toward the thickest trees was too long, resulting in about a 5 m 3 ha -1 overestimation in total volume.The additional dominant tree characteristics increased the peakedness of the distribution, which clearly improved the fit.In this case, the change was resulted in mainly by H dom because the observed D dom was close to its expectation.Calibrating with additional G slightly decreased but with d gM significantly increased the fit of the Weibull distribution as compared with the three common calibration variables (see Fig. 4).The two bestfitting distributions were as follows: 1. calibrated with dominant tree characteristics, which resulted in an error in total volume of 1.0 m 3 ha -1 and 2. calibrated with d gM , which resulted in an error in total volume of -1.8 m 3 ha -1 .If the same stand's prediction was done with G, d gM , and h gM (common for advanced stands) the resulted distribution differed only slightly from the expected one.Thus, calibrating the Weibull frequency distribution with the above variables did not improve the fit as efficiently as the arithmetic characteristics did and the dominant tree characteristics were also less efficient.Instead, calibration with the additional N or D resulted in a good fit.

Final Validation of the Family of Models
The linear prediction application provided a consistent basis for comparison of alternative size distribution models.The family of models was evaluated through use of the combinations of stand characteristics for the linear prediction proceeding from the common FMP variables.In young stands, the bias in total volume showed 1.3-4.6%underestimation with the three independent variables (see Table 7).In general, when additional predictors were included in BLUP models, the W N model yielded the least biased volume estimates.Inclusion of G-reduced biases most efficiently with basal area-dbh dis-tribution models -i.e.SB G and W G -whereas H dom was efficient with the dbh-frequency W N model (see Table 7).Note that whenever D dom was included, it resulted in clear underestimation in volume characteristics with the SB G and W G models for young stands.It was noted that the greatest underestimates were found when the difference between D and D dom was great (D dom > 2 × D).
When the W N model was used, the effect of each additional variable on the RMSE for volume and assortments was clear, most obviously when d gM was included, and smallest when G was included (see Table 8).However, the calibration of the basal area-based distribution models reduced the RMSE only slightly, most when H dom or d gM was included.All in all, d gM reduced the RMSE value for total and pulpwood volume by 3-12% and 6-16 percentage points, respectively.The combination for the best two additional variables was evidently d gM and h gM for W N and SB G models.Apart from those, H dom and D dom obviously enhanced the accuracy in the validated characteristics with the W N model while G and H dom did the same with the W G and SB G models.
In terms of bias and RMSE in the validated characteristics, W N was by far the best model for young stands.It was superior in terms of bias and RMSE in volume characteristics and basal area.However, dominant height was least biased with the SB G model.On average, the W G model generated the most accurate dominant diameter, providing the least bias for five cases and the smallest RMSE in six out of 11 cases (see Tables 7 and 8).
In advanced stands, total volumes were practically unbiased and the differences in RMSEs were rather small between the basal area-dbh distribution models whereas W N resulted in about 1% underestimation and systematically higher RMSEs than the basal area-dbh distribution models (see Tables 9 and 10).The bias and RMSE in commercial assortments were generally low, with standard three stand variables with the W G and SB G models (see Tables 9 and 10), but they were still reduced with the addition of D or N, whereas H dom reduced RMSE only (Table 10).Also, when the W N model was used, a considerably high bias of about 10% in timber assortments was reduced to 6-7% with additional D or N (Table 9).RMSE  in waste wood volume was about 30%, indicating inaccuracy in the left-hand tail of the distribution.
Including D improved the accuracy in the waste wood fraction greatly, especially with the W N and SB G distribution models.The W G model provided practically unbiased dominant tree characteristics, while the W N and SB G models produced slight biases (< 3%) in D dom .Biases in stem numbers were typically less than 2%.Surprisingly, D was superior to N in reducing RMSE in the stem number (see Table 10).
When the models were compared with each other in terms of bias, the W G model seemed slightly superior in the validated stand characteristics.However, SB G also performed well in volume estimates and W N yielded the least biased waste wood volume.The smallest RMSEs for stem number and waste wood was provided by the W N model.The two basal area-dbh distribution models produced nearly equal RMSEs in dominant tree characteristics and were superior to the dbh-frequency distribution model.

Discussion and Conclusions
One of the main advantages of the applied approach is the flexible prediction with varying sets of known stand variables.The essential predictors were the age of the stand, site factors, and the average temperature sum corresponding to the location of the stand.The temperature sum for a stand can be calculated according to the location in terms of latitude, longitude and altitude (Ojansuu and Henttonen 1983), but tabulated munici-pal average temperature sums are also available in MOTTI.The formulation of the prediction equations for eight stand characteristics was only slightly modified from the corresponding models for Norway spruce that were presented by Siipilehto (2006).It is worthy of mention that the models for Norway spruce were not sensitive to errors in stand age, if correct stand variables were available for linear prediction.Indeed, an assumed 25% bias in stand age resulted in the bias of only 2-5% in the five stand characteristics  Siipilehto 2006).
The expectations of the basic models provided reasonable development for stand characteristics by site and stand origin.In fact, dominant heights at a total age of 100 years for different site types (Cajander 1925) and stand origins seemed comparable with the existing Finnish H 100 site index curves, for naturally regenerated stands by Gustavsen (1980) and for artificially regenerated stands by Vuokila and Väliaho (1980).However, the above-mentioned site index curves showed slower early (T < 40 years) development in comparison with the basic model for H dom in this study.Similarly, the initial dominant height development for young artificially regenerated Scots pine stands as described by Varmola (1993) had a more pronounced sigmoid form and produced smaller H dom values for the first ten years, but the later development was comparable with our results.
In a MOTTI simulator, we have user-defined initial stand as an option.It is convenient for a non-experienced user to include the possibility of a simulator to generate expected age, species and site-dependent values for standard stand characteristics rather than force the user to guess reasonable values for them.However, a typical MOTTI application is to simulate young stands with various user-defined initial stand densities.The brief application for sapling stands showed that the effect of stand density on dbh distribution was reasonable and consistent with the densitycalibrated BLUP estimates for diameter-based stand characteristics.In addition, the effect of stand density on the stand characteristics was generally in line with the models of Varmola (1993) and Huuskonen and Hynynen (2006).In contrast, stand density had a much smaller effect on mean and dominant diameter or mean height in the models by Huuskonen and Miina (2007).
The two-parameter Weibull function was applied for both dbh-frequency distribution (W N ) and the basal area-dbh distribution (W G ) models motivated by the simplicity in weighting (Gove and Patil 1998) but also because Maltamo et al. (1995) reported that the two-parameter Weibull provided systematically better volume estimates than the three-parameter Weibull did.In line with this study's findings, Gobakken and Naesset (2004) and Palahi et al. (2007) noticed that determining two-parameter Weibull distributions for advanced stands as a basal area-dbh distribution instead of a dbh-frequency distribution results in more accurate stand characteristics.On the other hand, Maltamo et al. (2007) did not notice much difference between these two approaches.Furthermore, Gobakken and Naesset (2005) were surprised that the two-parameter Weibull could produce results comparable with the flexible percentile-based distributions.It is well known that the SB function is more flexible than the Weibull function (e.g.Hafley and Schreuder 1977).Nevertheless, the superior flexibility had a minor impact because the shape parameters of the SB distribution were not calibrated efficiently with the stand characteristics that were included.
Typical Finnish combinations of the variables assessed proved to be relatively efficient, but, concurrently, at least three variables were required for reliable stand volume estimates.For example, RMSE in the total volume for advanced stands was 39% with the basic model for W N distributions and it was reduced to 10% with two calibrating variables and further to 2% in calibration with the common three variables.In line with the findings of Siipilehto (1999), Kangas andMaltamo (2000a, 2000b), and Mabvurira et al. (2002), additional stem number knowledge in combination with G, d gM , and h gM produced variability in the shapes of predicted size distributions, which generally improved the fit.Knowledge of the stem number reduced the RMSE for N to 12-14% from the predicted distributions when the basal area was simultaneously correct.In comparison, distribution models by Kangas and Maltamo (2000b) and Siipilehto et al. (2007) including N as an additional predictor yielded a 5% RMSE in N with the independent INKA test data set.Maltamo et al. (2000) have shown the ability of medians of frequency (d M ) and basal area-dbh distributions (d gM ) to predict the shape of the Weibull frequency distribution and also multimodality in the case of predicted percentile-based distribution for Scots pine-dominated heterogeneous stands.When D and d gM were combined with other common FMP variables in this study, the accuracy in volume characteristics was considerably enhanced, especially with the W N model.In advanced stands, RMSEs of about 30% in waste  10).Furthermore, D had an even greater ability to improve the accuracy in the stem number generated than the inclusion of N itself did.Dominant tree characteristics showed a rather confusing calibration effect.On one hand, they showed an ability to reduce RMSEs for the modelled parameters of the distribution functions and height curve (see Fig. 2).Applying them for calibration resulted in reduced RMSEs for each validated stand characteristic in the test data when W N was used for prediction involving young stands (see Table 9).Furthermore, derived from the model by Siipilehto (2009), 1 / c of the dbhfrequency Weibull distribution could be assumed to be closely correlated with ln(D dom / D), which, in turn, can be written as ln(D dom ) -ln(D).However, D dom as a calibrating variable resulted in clearly biased volume estimates when basal areadbh distribution models were applied for young stands (see Table 8).In advanced stands, despite their ability to improve accuracy for dominant tree characteristics of the predicted distributions, their effect on the accuracy of stand volume was not as evident.The reason was probably related to the fact that dominant tree characteristics were too closely correlated with basal area medians so could not provide as useful additional information.The ratio between dominant and mean characteristics has been used previously for predicting the Weibull distribution by Siipilehto (2009), while Sarkkola et al. (2005) used a similar formulation but utilised the maximum diameter instead of the dominant diameter.
In a study by Siipilehto (1999), parameter β 1 of Näslund's function was constructed to be linearly dependent on d gM and h gM , and β 0 was solved for such that the curve crossed the point (d gM ,h gM ).This height model has been used in numerous studies since (e.g.Kangas and Maltamo 2003) and recently, Korpela and Tokola (2006) found a slight bias when using it.Note that the more representative data of the present study clearly showed that the dependence between these variables was non-linear in the original scale, and linearisation between the parameters of Näslund's curve and stand characteristics required logarithmic trans-formation on both sides (see fig. 2 in Siipilehto 2011).The BLUP estimated height curves in this study passed close to the point of means (D,H) or medians (d gM ,h gM ), while the knowledge of D dom and H dom was useful in characterising the bending of the curve (see Eerikäinen 2003).In contrast, stand density did not show any effect on the accuracy of the height curve parameters (see Fig. 2C).
Naturally, the models presented for stand characteristics can be used independently of the size distribution models given, in combination with any existing size distribution model.In that case, we may need the prediction for d gM , h gM , and G for young stands, or N for advanced stands, in order to be able to use the prediction model presented, for example, by Siipilehto et al. (2007) or Kangas and Maltamo (2000a).The height distribution model describes the complete distribution, whereas that of the dbh distribution describes the right tail of the truncated distribution until all the trees reach breast height.When a mean height of 4 m was reached, practically all of the crop trees had reached breast height.If dbh distribution is used for a stand in its early state, instead of the total number of stems, the number of stems above breast height should be used when one is calculating stand basal area.However, it is recommended to apply the height distribution model by Siipilehto (2009) up to a mean height of about 4 m because of its continuous feature.
The approach in this study shows one option for flexible description of stand structure.It is flexible in the sense that size distributions and height-dbh relationships can be predicted in 255 different ways, if the eight stand characteristics that were modelled in this study can be combined freely as explaining (calibrating) variables.The flexibility in combining the stand characteristics may become of great interest whenever new techniques (e.g.laser scanning) find new combinations of variables to be more efficient or accurate to assess as compared with the stand characteristics traditionally collected in the field.
According to this study, it would appear essential to choose the size distribution model such that the available data and the model applied are based on the same scale.Thus, if the basis lies in the common FMP stand characteristics, frequency distribution is superior to basal area distribution for young stands and, vice versa, basal area distribution is superior to frequency distribution in advanced stands.In general, considerable improvements in the accuracy of the predicted distribution is possible when arithmetic mean and basal area-median diameters are simultaneously included in the prediction model.It is worth noting, however, that some combinations can also result in severely biased estimates of volume, as was shown when the W G and SB G models were applied for young stands with D dom used in addition to the common variables N, D, and H.This was probably due to violation of the assumed linear dependence of the crossmodel errors.Such violation was not noticed with respect to the W N model for young stands or any model for advanced stands.Thus, this was indirect evidence that the assumed linearity between cross-model errors held in general.
The presented BLUP models are intended to replace the previous generation of models (e.g. Siipilehto 2006) in the MOTTI simulator.One methodological drawback was the fact that we did not manage to include the effect of repeated measurements in the model, most probably due to the large number of simultaneous models.However, according to Liu et al. (2004) the influence can be considered small.In addition, Siipilehto et al. (2007) showed the comparable results with SUR and mixed-effect SUR models when predicting SB distribution, whereas they were both superior to OLS and separately fitted mixed models.
The accuracy of the present model in terms of total stand volume and timber assortments should be carefully compared with the other models available in order to find the most reliable combinations of models.The linear prediction application presented provides a theoretically consistent basis for comparing alternative size distribution models with alternative combinations of known stand characteristics.Also, the option of utilising methods for compatibility in the known moment or percentile in order to improve the fit of the presented dbh distribution is worth exploring (see Merganič andSterba 2006, Gobakken andNaesset 2004).Obviously, the optimal distribution model and the best set of variables changes with the stand's stage of development and with respect to the specific stand characteristic of interest (Kangas and Maltamo 2003).Calibration yielded an error variance of 0.0478.The corresponding error variance for β 1 was 0.0062 for the calibrated model.
The corresponding OLS model for the parameters was estimated by including ln(G), ln(d gM + 1), and ln(h gM ) in the basic model for β 0 and β 1 .The degrees of determination of the OLS models were 52% and 87%, and the error variances were just above those of the BLUP model, at 0.0483 and 0.0064 for β 0 and β 1 , respectively.The OLS estimates were quite similar to the BLUP estimates, as shown below.
a) β 0 , β 1 , parameters of the Näslund's height curve; b, c, parameters of the Weibull function for frequency distribution; b G , c G , those for Weibull basal area-dbh distribution; λ, δ, γ, parameters of Johnson's SB function for basal area-dbh distribution.
).Because parameter b represents the 63rd percentile of the Weibull distribution, b of W N was highly correlated with mean D (r = 0.97) and b G of W G with basal area median dbh, d gM (r = 0.99).There was no such close correlation between the residuals of the shape parameter c and individual stand characteristics.The highest correlation was found between c and D (r = -0.53).According to correlation scatter plots in Fig.1, the assumption of the linear correlation of BLUP held quite nicely.

Fig. 1 .
Fig. 1.Correlations of the cross-model residual errors between the modelled stand characteristics and the parameters.Stand characteristics from left to right were G, N, D, d gM , H, h gM , D dom and H dom , and the parameters from top to bottom were β 0 and β 1 of the Näslund's curve, parameters b and c of the W N , b G and c G of W G and λ and δ of the SB G distribution.

Fig. 2 .
Fig. 2. The effect of BLUP estimation on the RMSE for the modelled parameters.Calibration was based on stand characteristics which were either mainly arithmetic (A, B) or basal area-weighted (C, D).The parameters are β 0 and β 1 of the Näslund's curve, b and c of the Weibull frequency distribution or that of basal area-dbh distribution (W G ) and λ and δ of the SB distribution.Variables used in calibration are shown in the legend.

Fig. 4 .
Fig. 4. Example of the BLUP estimation for the Weibull frequency distribution for a 21-year-old planted Scots pine stand on MT site in southern Finland.W(0) represents expected Weibull distribution of the basic model.W(N,D,H) was calibrated with N, D and H and W(+D dom ,H dom ) with additional dominant tree characteristics, etc.

Table 1 .
Stand variables and estimated parameters a) for the height curve and dbh distribution for Scots pine stands.P denotes the proportion of Scots pine out of the stand total basal area.

Table 2 .
The mean and range of the validated stand characteristics for young (H < 9 m) and advanced (H ≥ 9 m) stands from validation data.
D < d gM < D dom and H < h gM < H dom ), even though some of them crossed in checking the first step OLS models.Using the transformation ln(Y + 1) for diameters and the various powers in the transformation of age for dbh characteristics (T

Table 5 .
Cross-model error covariances above the diagonal, variances on the diagonal (in bold) and correlation coefficients below the diagonal (italics) for the modelled stand characteristics.
r = -0.6 between β 1 and h gM and H dom

Table 6 .
Cross-model error covariances and correlation coefficients (italics) between modelled stand characteristics and parameters.Correlations | r | > 0.5 are highlighted in bold.
the prediction for λ benefited each of the stand characteristics included but the RMSE for δ was reduced only marginally until the additional stem number or additional mean dbh was included.It is worth noting that the effect of stand charac-teristics on the estimate for γ of the SB distribution was quite marginal (figure not shown).The RMSE of 80% for γ was reduced somewhat only when additional D dom (78-64%), N (61%), D or d gM (60%) was included in the BLUP model.

Table 7 .
The

Table 8 .
The RMSE (%) for stand characteristics using the alternative distribution models for young stands.The first row represents SOLMU variables (N, D, H) and below that are the additional stand variables used for predictions.The best validating criterions are highlighted in bold.

Table 9 .
The bias (%) in stand characteristics using the alternative distribution models for advanced stands (H ≥ 9 m).The first row represents SOLMU variables (G, d gM , h gM ) and below that are the additional stand variables used for predictions.The best validating criterions are highlighted in bold.

Table 10 .
The RMSE (%) for stand characteristics using the alternative distribution models for advanced stands.The first row represented SOLMU variables (G, d gM , h gM ) and the variables below represented the additional information used for predictions.The best validating criterions are highlighted in bold.
SiipilehtoLocal Prediction of Stand Structure Using Linear Prediction Theory in Scots Pine-Dominated Stands in Finland wood volume indicated inaccurate left tails of the predicted distributions.Additional D improved the accuracy of waste wood volume most efficiently (Table