METSÄNTUTKIMUSLAITOKSEN TIEDONANTOJA 929, 2004 FINNISH FOREST RESEARCH INSTITUTE, RESEARCH PAPERS 929, 2004 Predicting Stand Characteristics Using Limited Measurements Lauri Mehtätalo JOENSUUN TUTKIMUSKESKUS JOENSUU RESEARCH CENTRE Predicting Stand Characteristics Using Limited Measurements Lauri Mehtätalo ACADEMIC DISSERTATION To be presented, with the permission of the Faculty of Forestry of the University of Joensuu, for public criticism in auditorium K1 of the University, Yliopistokatu 2, Joensuu, on 24 th September. 2004 at 12 o'clock noon. Joensuu 2004 Mehtätalo, Lauri. 2004. Predicting stand characteristics using limited measurements. Metsäntutki muslaitoksen tiedonantoja 929. Finnish Forest Research Institute, Research Papers 929, 39+72 p. ISBN 951-40-1934-2. ISSN 0358-4283. Publisher The Finnish Forest Research Institute Joensuu Research Centre Accepted by Kari Mielikäinen, Research Director, in September 2004 Supervisors Prof. Annika Kangas University of Helsinki Dr Juha Lappi Finnish Forest Research Institute Prof. Matti Maltamo University of Joensuu Reviewers Dr Jeffrey H. Gove USDA Forest Service Prof. Timothy G. Gregoire Yale University Examiner Dr Risto Ojansuu Finnish Forest Research Institute Copyrights of original publications Papers I and II © Society of American Foresters Papers 111 and V © NRC Research Press Paper IV © Finnish Society of Forest Science and Finnish Forest Research Institute The papers are reproduced by kind permission of the copyright owners. Front cover Distributions of sample order statistics (coloured lines) when the sample was drawn from a percentile-based diameter distribution (black line). Photo: Erkki Oksanen/FFRI Author's contact information Address: Finnish Forest Research Institute, Joensuu Research Centre, P.O. Box 68, Fin-80101 Joensuu, Finland Phone: +358-10-211 3051 email: Lauri. Mehtatalo@metla.fi Fax: +358-10-21 1 31 13 Yliopistopaino Joensuu 2004 Abstract This thesis reports a new system for the production of static stand description in an inventory by compartments. The stand description includes stock density, diameter distribution and height diameter (H-D) curve. The diameter distribution of the stand is expressed with percentiles. Firstly, expected percentiles are predicted with regression models using measurements of stand variables. Secondly, the pre dicted percentiles are localized for the stand using order statistics of horizontal point sample plots (HPS-plots) (i.e. quantile trees), which are interpreted as measured percentiles of the stand. Thirdly, the obtained localized percentiles are adjusted in order to ensure compatibility with the measured stem number. The expected H-D curve of the stand is predicted using the measured stand variables. Furthermore, it is localized for the stand using height sample trees. The longitudinal character of the model makes it possible to use measurements from several points in time. Both the localizations of diameter distribution and of the H-D curve are based on the prediction of random effects of the models with the best linear unbiased predictor using sample measurements of the response. The individual components are combined as a new system for the prediction of stand description. A key feature of the system is its ability to utilize different amounts of input information. Further more, measurement errors of stand variables are utilized to some degree. The minimum input of the system, which can be obtained from one HPS-plot, consists of measurements of the basal area, basal area median diameter (DGM), stand age and site fertility class. Additional UPS-plots, stem number measurements) from fixed plot(s), old or new height sample trees and quantile trees can also be utilized. The system makes it possible to take more measurements from stands with a high accuracy requirement than from stands with a low accuracy requirement. The system was utilized in estimating a model of expected errors of predicted volume and saw timber volume using different measurement strategies in different stands. The prediction error depended on the basal area and DGM of the stand and on the number of UPS-plots, height sample trees and quantile trees. Furthermore, a height measurement from a previous inventory decreased the prediction errors slightly. The measurements of stem number did not significantly improve the accuracy of volume and saw timber volume predictions. The estimated models were used as objec tive functions in a constrained optimization problem, where the object was to find an optimal measurement strategy for a single stand in an inventory where measurement time is limited. Key words: stand structure, diameter distribution, height-diameter, order statistics, linear predic tion. longitudinal analysis, optimization, forest planning, calibration, adjustment, localization. Preface This work was strongly based on the previous work of my supervisors. Matti Maltamo, Annika Kangas and Juha Lappi. Their ideas and earlier findings have had a strong influence on all stages of this work. The work of Annika Kangas and Matti Maltamo on diameter distributions and stand wise inventory, and the work of Juha Lappi on applications of mixed models and linear prediction in forestry, provided a good toybox for a young researcher to play with. My work has actually only consisted of taking different kinds of building blocks from this toybox, building new wobbly tow ers and trying to get them to stand upright. The supervisors have always been available for discus sions and comments about the work being done. I am grateful for this unfailing support, which I have received thoughout the various stages of this work. I also wish to thank my pre-examiners Jeffrey H. Gove and Timothy G. Gregoire for their valuable comments on the work done. The funding for this work was provided by the Academy of Finland (decision numbers 73392 and 200775). Annika Kangas, Jyrki Kangas and Kari T. Korhonen applied for the funding and managed the research projects. The work was carried out at the Finnish Forest Research Institute, Joensuu Research Centre, where Tuula Nuutinen, Jari Parviainen and Jyrki Kangas have provided me with good facilities for carrying out this research. My position as an "observing associate member" of the MELA-group, led by Tuula Nuutinen, has given me an insight into the development of a forest planning system in practice and provided me with the contacts to help me with practical problems. The datasets used in this study have been collected by Timo Pukkala and FFRI. Pekka Leskinen, Juha Alho, and a group of researchers on forest assessment at the University of Joensuu, have given me constructive comments about origi nal manuscripts. The course on Statistical Inference, held by Jukka Nyblom in 2003, was a source of inspiration for me. Without this course, Paper I would hardly have been written. Lisa Lena Opas- Hänninen revised the language of all parts of this thesis. Our telephone calls during the process were valuable private English lessons for me. Coffee breaks with Eero Muinonen, Ulla Mattila and Reetta Lempinen let me forget my research for a moment. The office employees of the research centre helped me always when I forgot my key at home or in my room. Leena Karvinen and Seija Partanen helped me with the makeup of this thesis. I am grateful to all these people for their contri bution to this work. However, my life is mostly outside the office. I give my warmest thanks to my wife. Hilkka, for her understanding and support during this work, and to my children, Juho and lida, who use real building blocks for their towers, for reminding me about where the real life is. Their bright eyes and joyful voices remind me every day who and what I am living for. Joensuu. August 2004 Lauri Mehtätalo List of articles This thesis comprises the present summary and the following articles, which are referred to in the text by Roman numerals I-V. I Mehtätalo, L. Forthcoming. Localizing predicted diameter distribution with sample information. Revised manuscript (Forest Science). II Mehtätalo, L. 2004. An algorithm for ensuring compatibility between estimated percentiles of diameter distribution and measured stand variables. Forest Science 50(1): 20-32. HI Mehtätalo, L. 2004. A longitudinal height-diameter model for Norway spruce in Finland. Canadian Journal of Forest Research 34(1): 131-140. IV Mehtätalo, L. Forthcoming. Height-diameter models for Scots pine and birch spe cies in Finland. Submitted manuscript (Silva Fennica). V Mehtätalo, L. and Kangas, A. Forthcoming. An approach to optimizing field data collection in an inventory by compartments. To appear in Canadian Journal of For est Research. Paper V was planned jointly with Kangas. I built the system, calculated the results and was mainly responsible for writing the article. Mehtätalo 6 Contents 1 Introduction 7 2 Components of stand description - a literature review 9 2.1 Diameter distribution 9 2.1.1 Approaches 9 2.1.2 Distribution families 10 2.1.3 Compatibility of diameter distribution 11 2.2 Height-diameter models 11 2.2.1 Approaches 11 2.2.2 Model forms 12 2.3 Other components 13 2.4 Other approaches 13 3 Percentile-based diameter distribution 15 3.1 Distribution function and density 15 3.2 Moments of the distribution 16 3.3 Relationships between percentiles and stand characteristics 17 3.4 Order statistics 17 3.5 Considerations on the PPM with the percentile-based approach 17 4 Datasets of this study 19 4.1 Mixed forests data (I, II) 19 4.2 INKA-data (m, IV, V) 19 5 Methods 20 5.1 Mixed models (111, IV) 20 5.2 Prediction of random variables (I. 111. IV. V) 20 5.3 Constrained optimization (H, V) 21 6 Summary of results 23 6.1 Localizing predicted diameter distribution with sample order statistics (I. V) 23 6.2 Adjusting the predicted percentiles to obtain a compatible stand description (11, V) 24 6.3 H-D models from longitudinal data (111, IV. V) 25 6.4 A system for producing stand description in an inventory by compartments (V) 25 6.5 Optimizing data collection in an inventory by compartments (V) 26 7 Discussion 27 7.1 The system for producing a stand description 27 7.2 The use of sample information 28 7.3 Optimal allocation of stand measurements 29 8 Reducing the costs of the inventory for forest planning 31 References 32 Introduction 7 1 Introduction The aim of forest planning is to find a man agement strategy for a forest area that maxi mizes the utility for the forest owner (Pukkala 1994). The traditional primary unit of forest planning in Finland is a forest stand. The forest plan includes management suggestions for every stand of the forest area under considera tion. In order to collect the data for the plan, an inventory by compartments is carried out (Poso 1983). A Finnish inventory by compartments is based on few horizontal point sample plots (HPS-plots. i.e. relascope sample plots, angle count sample plots, Bitterlich plots). The plots are established subjectively by the person carrying out the inventory at locations that seem representative for the stand. Using measurements and visual assessments of the plots, the most important characteristics, including basal area, basal area median diameter (DGM), height of a DGM-tree, stand age and site fertility class are assessed from each stratum of the growing stock in the stand (Paananen et ai. 2000). The simulator of a forest planning system utilizes the stand wise characteristics to produce a stand description, including stock density, diameter distribution and height-diameter (H-D) curve of the stand. In the simulator of the MELA-system (Redsven et al. 2004), which is the core of most forest planning systems in Finland, they are predicted using the models of Mykkänen (1986), Veltheim (1987), Kilkki "et al. (1989), Siipilehto (1999) and Kangas and Maltamo (2000b). The obtained stand description is used to generate a set of representative trees for each stand for the prediction of growth (Hynynen et al. 2002) and cutting removal in alternative management schedules. The accuracy of the input of the Finnish stand simulation system based on partially visually assessed stand characteristics is re ported to be rather low (Poso 1983, Mähönen 1984. Laasasenaho and Päivinen 1986, Pussi nen 1992, Pigg 1994, Kangas et al. 2002. Haara and Korhonen 2004, Kangas et al. 2004), and prediction errors of growth models decrease it further (Kangas 1997, 1998 a, 1998b, 1999). The data users are, however, satisfied with the accuracy, but would not like to see it decline (Uuttera et ai. 2002). The aim of Finland's National Forest Programme 2010 (1999) is to increase the cover of forest plans from 50% to 75% by 2010. This requires decreasing the costs of planning per unit area (Heikinheimo 1999, Paananen 2002, Saramäki et ai. 2003). which means that the current level of accuracy should be retained at a lower cost than previ ously. One possibility to respond to these needs is to develop methods and models that provide as accurate assessments as the current inventory system but at a lower cost. Many studies have investigated possibilities to carry out the inven tory from the air using aerial photographs (e.g. Pitkänen 2001, Anttila 2002 a, Anttila and Le hikoinen 2002, Korpela 2004), satellite im agery (e.g. Hyvönen 2002, Saksa et al. 2003) or laser scanning approaches (e.g. Holmgren et al. 2003. Maltamo et al. 2004 a). Even if these approaches are promising, they are not yet real alternatives to the inventory by compartments (Uuttera et al. 2002) and their development will take time. In other studies, the information of the previous inventory has been updated by utilizing growth models and information about treatments from forestry databases (Hyvönen and Korhonen 2003) or aerial photographs (Anttila 2002b) in order to lengthen the interval of two inventories. These approaches can be regarded as practical variations of the approach of Stälil (1994), where it was suggested that if the expected utility of an inventory was greater than its costs then a new inventory of a forest stand should be earned out. Another possibility is to search for possible new variables to be measured in the field that would provide more information at a lower cost than the currently measured stand vari ables. hi many studies carried out in Finland, using accurate stem number in the calculations has been found to be useful (Siipilehto 1999, 8 Mehtätalo Kangas and Maltamo 2000b). However, sev eral studies (e.g. Kangas et ai. 2004) have re ported large errors in measurements of stem number. The calibration algorithm of Kangas and Maltamo (2000 a) made it possible to util ize different stand variables in the prediction of diameter distribution and they found that of several potential new variables the unweighted median, maximum and minimum diameters were the most promising alternatives for new measurements. In addition, new measurement equipment is being developed in order to make it possible for a single person to measure a large sample of diameters and heights rapidly in the field (Laasasenaho et al. 2002, Koi vuniemi 2003). Eid (2000) reported that losses caused by mistiming of harvests are most serious in young stands and in stands that are close to their economically optimal rotation age, while in middle-aged and over-mature stands the losses are smaller. Thus, it was suggested that the inventory data should be most accurate in stands where expected losses are greatest. In addition. Kangas and Maltamo (2002) pro posed that different variables should be meas ured in different kinds of forests. Furthermore, Holmström et al. (2003) reported that in large stands, intensive field sampling should be car ried out while in smaller stands a less intensive inventory might be satisfactory. All these stud ies indicate that money could be saved by vary ing the forest inventory strategy from stand to stand. This study aims at reducing the costs and improving the accuracy of the traditional in ventory by compartments by looking for new measurements, making the use of the collected data more efficient and allocating the meas urement resources to measurements and stands where the utility is greatest. This requires a calculation system that is able to utilize various input information. The current system for pro duction of stand description requires a fixed set of variables from each stand. Furthermore, it assumes that the input information is measured without errors, even if it is well known that the errors of stand measurements are large. Thus, in this study a new system for predicting stand description was developed. In predicting forest characteristics, the use of regression models has become very popular. Some easily measurable stand variables, the most common of which are basal area, mean diameter, stand age and site type, are used to predict other stand characteristics that cannot be measured, or at least are difficult to meas ure, in the field. Examples of these are diame ter distribution (e.g. Rennolls et al. 1985), tree height (e.g. Fang and Bailey 1998), stand growth (e.g. Woollons 1997), log volume re duction (e.g. Mehtätalo 2002), damages (e.g. Jalkanen and Mattila 2000) and the production of berries (e.g. Ihalainen et al. 2003). However, the most natural means to estimate any stand characteristic is to measure it. Thus, if meas urement of the target variable is possible, it should be preferred over the measurements of covariates of a regression model (see e.g. Lappi 1997). This study applies this principle to di ameter distributions and H-D curves. The ran dom parameter approach and linear prediction theory provide effective tools for carrying it out (Lappi and Bailey 1988. Lappi 1991). The aims of this study were • To develop tools for producing a stand de scription. These tools include methods for predicting diameter distribution and models for the H-D relationship. The tools should be able to utilize different kinds of information measured at different levels of accuracy. Special attention is paid to the effective use of sample tree information. (I - IV.) • To utilize the tools developed in order to construct a new system for producing stand description in forest planning in Finland (V). • To utilize the system in the optimization of data collection in an inventory by compart ments (V). Components of stand description - a literature review 9 2 Components of stand description - a literature review Stand description includes the information of the stand measurements in such a form that it can be used in stand simulation. An example of a very coarse description includes the name of the main tree species while a very detailed description may include a complete tree map of the stand with measured tree taper. However, stand description is always a simplification of the reality and the extent of simplification depends on the purpose of the stand description and the data available. The stand description of this study includes the stock density, diameter distribution and H-D curve. The stock density is described either by the number of stems or by the basal area. In the Finnish system, the basal area is used because it can be measured with higher accuracy than the stem number and it is more strongly corre lated with the total volume than the number of stems is. The measurement of the basal area of a sample plot is straightforward with an angle gauge and the basal area of the stand is calcu lated as the mean of the plot wise basal areas. The measurements of diameter distribution and the H-D relationship, on the other hand, are rather laborious to carry out in the field. Thus, these components are not measured in the field; instead, approaches for predicting them from stand and sample tree measurements are used. The next two subsections consider the ap proaches reported in the relevant literature. 2.1 Diameter distribution 2.1.1 Approaches The diameter distribution is the basis of the stand description. Many approaches have been used to construct the diameter distribution of a stand. In the following, the approaches are divided into three main approaches: (i) those based on a sample of diameters, (ii) those based on prediction or recovery of the parame ters of an assumed theoretical distribution model and (iii) those based on known diameter distributions of similar stands (imputation methods). Approach (i). The most natural way to obtain the diameter distribution is to measure a diame ter sample from the stand. If the sample is large enough, it can be used as such in the simulation (e.g. Pienaar and Harrison 1988, Nepal and Somers 1992, Tang et al. 1997). If the sample is small, it can be smoothed (e.g. Droessler and Burk 1989, Uuttera and Maltamo 1995) or a theoretical distribution function can be fitted to it (e.g. Bailey and Dell 1973, Zarnoch and Dell 1985, Van Deusen 1986, Shiver 1988, Lindsay et al. 1996, Zhou and McTague 1996, Scolforo et al 2003, Zhang et al. 2003). The theoretical distribution function can be fitted using the maximum likelihood method, the method of moments, methods based on linear regression or by utilizing properties of certain percentiles or stand variables. Approach (ii). The measurement of a diame ter sample is too time-consuming in many inventories. In these cases, the diameter distri bution can be predicted with some easily measurable stand variables. Traditionally, these methods are divided into parameter prediction methods (PPM) and parameter recovery meth ods (PRM) (Hyink and Moser 1983). In the parameter prediction method, the parameters of the assumed distribution function are predicted with some measured stand variables using estimated regression models (e.g. Schreuder et al. 1979. Little 1982, Rennolls et al. 1985. Kilkki and Päivinen 1986. Kilkki et al. 1989. Maltamo 1997, Siipilehto 1999. Temesgen 2003, Robinson 2004). In the PRM the parame ters are recovered from some stand variables using known relations between the stand vari ables and distribution parameters (e.g. Ek et al. 1975. Burk and Newberry 1984. Magnussen 1986. McTague and Bailey 1987. Kuru et al. 1992). The stand variables used in PRM may be. for example, percentiles or moments of the diameter distribution. The recovery is. how ever, possible only for as many parameters as there are measured stand variables that are 10 Mehtätalo linked with the diameter distribution. If the number of parameters is greater, partial recov ery can be used, i.e. as many variables as pos sible are recovered and the other parameters are predicted (e.g. Kilkki and Päivinen 1986, Kan gas and Maltamo 2000b, I, 11. V). Approach (iii) The third approach is to use known diameter distributions of similar stands as the predicted diameter distribution of the stand (e.g. Haara et ai. 1997, Maltamo and Kangas 1998). These methods have also been called imputation methods (Ek et al. 1997, Temesgen 2003, Temesgen et al. 2003). The similar stands are selected from a neighbor hood that is defined with a distance function. The imputation methods used in predicting the diameter distribution are the method (Altman 1992) and the most similar neighbor method (Moeur and Stage 1995). The approach of Nanos and Montero (2002), where an interpolated surface is used to carry infor mation from geographical neighbors to the target stands, also belongs to this class of ap proaches. The above approaches can also be combined. For example, the neighbors of the third ap proach may be smoothed diameter distributions instead of true distributions, thus combining approaches (i) and (iii) (Maltamo and Kangas 1998). Furthermore, sample information can be used to improve a predicted distribution; ex amples of this are the Bayesian approach of Green and Clutter (2000) where the prior in formation of neighboring stands (iii) is com bined with sample information (i) and the ap proach of Paper I in the present thesis where sample information (i) is used to improve pre dicted diameter distribution (ii). The approach of Maltamo et al. (2003 a) combining empirical distributions of large trees identified from a digital video imagery with predicted distribu tion of small trees is a combination of (i) and (ii). 2.1.2 Distribution families No theoretical results have been presented regarding which distribution family should be used as a diameter distribution. Hence, many distribution families have been used. The most important factors affecting the goodness of fit of a distribution family is the number of free parameters and the flexibility of the distribu tion to represent all possible distributional forms. The most commonly used distributional fam ily is the Weibull distribution (Bailey and Dell 1973), either in a three parametric form with parameters for location, scale and shape, or in a two-parametric form, where the location pa rameter has been given a value of zero, i.e. the minimum diameter of the stand is assumed to be zero. An unrealistic property of the Weibull distribution is that it has no upper limit, i.e. the maximum diameter of a stand is infinity. This problem can be solved by truncating the distri bution at a point that is regarded as the maxi mum possible diameter of the tree species. Another possibility is to use a reversed trun cated Weibull distribution, so that the location parameter receives an interpretation of maxi mum diameter and the minimum diameter is zero (e.g. Kuru et al. 1992, Robinson 2004). Other distribution families used are. for exam ple, Gram-Charlier (Cajanus 1914), normal (e.g. Nanang 1998), lognormal (Bliss and Re inker 1964), gamma (Nelson 1964), beta (e.g. Hafley and Schreuder 1977), the Johnson's distribution families (e.g. Zhou and McTague 1996), the Chaudhry-Ahmad family (Chaudhry and Ahmad 1993, Nanos and Montero 2001) and the exponential distribution (e.g. Cancino and Gadow 2002). Many of these families also need to be truncated in applications because of unlimited maximum and minimum diameters. All the distribution families presented below are. however, too rigid in some stands. This happens, for example, in stands where the shape of the distribution is irregular or multi modal. Thus, approaches that combine several distributions have been presented. These ap proaches are the segmented distribution ap proach (Cao and Burkhart 1984), finite mixture approach (e.g. Zhang et al. 2001) and the per centile-based approach (e.g. Borders et al. 1987). The segmented distribution is constructed from pieces of the selected distribution [Weibull in Cao and Burkhart (1984)]. and the number of pieces and locations of cutting Components of stand description - a literature review 11 points are fixed before the fitting procedure. In the finite mixture approach, the density func tion is a weighted sum of several density func tions of the selected distribution family. In the reported studies, the number of distributions has been two and the theoretical distribution used has been either the Weibull distribution (Zhang et al. 2001, Liu et al. 2002) or the bivariate normal distribution (Zucchini et al. 2001). The percentile-based approach ex presses the distribution using a fixed number of percentiles, corresponding to predefined values of the distribution function. The continuous function is obtained by interpolation, either by linear (e.g. Borders et al. 1987, I. II) or by spline interpolation (Maltamo et al. 2000, Kan gas and Maltamo 2000 c). When linear interpo lation is used, the obtained distribution consists of pieces of a uniform distribution (I). Thus, it can be regarded either as a segmented distribu tion approach or as a finite mixture approach using uniform distribution (See Section 3). 2.1.3 Compatibility of diameter distribution Compatibility of stand description means that all stand variables calculated from the diameter distribution are equal to their meas ured variables. It is a desired property of the stand description and incompatibility may indicate that the information of the stand meas urements is not utilized effectively. In PRM, the compatibility can be guaranteed with re spect to as many stand variables as there are parameters in the distribution family used. However, it usually leads to a very complicated set of equations, whose solution does not nec essarily exist in closed form. Compatibility of basal area weighted and unweighted diameter distributions can be guar anteed through size-biased distribution theory (Gove and Patil 1998. Gove 2000, Gove 2003), where the relationships between the ordinary and basal area weighted forms of the diameter distribution are derived analytically. It makes it possible to derive the parameters of the ordi nary distribution from the parameters of the basal area weighted distribution and vice versa. Thus, it provides tools for the recovery of the parameters of weighted distributions. In PPM, the compatibility is very hard to en sure. Thus, many studies have presented algo rithms that aim at a compatible stand descrip tion by adjusting the diameter distribution (Nepal and Somers 1992, Cao and Baldwin 1999, Kangas and Maltamo 2000 a). In these algorithms, the adjustment is applied to the frequencies of the stand table and no new di ameter classes are established in the adjust ment. Furthermore, the stand variables calcu lated using the adjusted frequencies are re quired to equal the measurements exactly. These requirements may be too strict in prac tice, where measurements include errors. An alternative for these approaches is presented in Paper 11. 2.2 Height-diameter models 2.2.1 Approaches The H-D model is used to predict heights for trees with given diameters. The H-D relation ship varies considerably between stands (Lappi 1997, Hökkä 1997. Jayaram and Lappi 2001, Eerikäinen 2003. Calama and Montero 2004, 111, IV) and accuracy of the predicted H-D curve has a considerable effect on the accuracy of the stand volume estimate. As with diameter distributions, the approach used in height pre diction depends on the data available and the possible approaches could be classified in a manner similar to the classification of the pre vious section, i.e. into (i) approaches based on a sample of heights, (ii) approaches predicting the parameters of the H-D curve without height measurements and (iii) approaches based on imputation. The studies using approaches be longing to the second (ii) category can further be divided into two classes: (ii-a) approaches which assume that a large forest area can be divided into stands, each of which has its own H-D curve, and (ii-b) approaches where a common H-D curve is assumed for larger ar eas. for example for states or regions. Approach (i). The most natural approach in predicting the H-D curve of a stand is to fit an assumed curve to observed H-D data from the target stand. There are numerous studies con cerning fitting different kinds of curves on 12 Mehtätalo observed H-D data (Curtis 1967, Omule and McDonald 1991. Arabazis and Burkhart 1992, Flewelling and de Jong 1992, Fang and Bailey 1998). The aim of these studies has been to find an appropriate functional form for the H-D curve and to show how the estimation should be done. A recent study of Zhang et al. (2004), which models the spatial variation within stand using geographically weighted regression, provides a new view to this approach. Approach (ii-a). The variation of H-D curves between stands is usually taken into account by modeling the parameters of the H-D curve using stand specific variables as predictors (Veltheim 1987, Borders and Patterson 1990, Parresol 1992, Lynch and Murphy 1995, Fang and Bailey 1998, Knowe et al. 1998, Wang and Hann 1998. Hanus et al. 1999, Zeide and Vanderschaaf 2002). Furthermore, these mod els have been localized using measured heights, for example by re-scaling the model so that the measured mean height is obtained when the diameter equals the measured mean diameter, thus combining approach (ii-a) with approach (i). However, even though this ap proach of localization works rather well when the mean diameters and heights are accurate, it does not take into account the within-stand and between-stand variances of tree heights. In recent studies, in addition to the use of stand specific predictors, the hierarchy of the data has been taken into account through a mixed model approach. It provides an effective and theoretically justified means for localizing the curves for a given stand using measured height(s) (Lappi 1997. Hökkä 1997, Jayaram and Lappi 2001, Calama and Montero 2004). The localization is possible even using just one measured height and the ratio of within-stand and between-stand variances determines how close the expected curve is to the measured height. The mixed model approach also pro vides natural approaches for taking into ac count the temporal development of the H-D curves (Lappi 1997, Eerikäinen 2003, 111. IV). Approach (ii-b). In models belonging to this approach, single values of the parameters of an assumed H-D model are estimated from a large dataset consisting of measurements from sev eral locations of the target area (Huang et al. 1992, Zhang 1997, Peng 1999, Huang et ai. 2000. Zhang et ai. 2002, Colbert et ai. 2002). When the model is used for prediction, the estimated parameter values are applied to the whole target area. Approach (iii). With the exception of Mal tamo et al. (2003b), imputation methods have not been used in the prediction of H-D curves, even though this could be done concurrently with the imputation of diameter distributions. The reason for this seems to be data related rather than methodological: the datasets used in the imputation of diameter distributions have not included tree heights. 2.2.2 Model forms Many studies carried out on H-D curves have compared different functional forms for the H- D relationship. The number of functional forms tested exceeds 30 and no single form has been found to be superior. However, some forms have been among the best ones in many com parisons either in a linearized or nonlinear form. The three most frequently used functions are the allometric function, also called the power function (Greenhill 1881, Curtis 1968, Zeide and Vanderschaaf 2002, Eerikäinen 2003. Zhang et al 2004) Meyer's equation (Meyer 1940. Stout and Shumway 1982. Farr et al. 1989. Huang et al. 1992. Fulton 1999) and the Korf curve, also called the Lundqvist or exponential function where H is tree height, D is diameter and a, b and c are parameters. Based on the biological growth pattern of a tree, Yuancai and Parresol (2001) recom mended the use of functions with an inflection H=aD" , (1) H=a( l-e 6D ) (2) H=ae~ bD , (3) Components of stand description - a literature review 13 point, which means that the curve is S-shaped. and has an upper asymptote. The allometric equation lacks both these properties. However, it has a strong mechanical basis: if the exponent b is given a value of 2/3, the stem is equally resistant to bending at dif ferent heights (Greenhill 1881. p. 66-73; see Zeide and Vanderschaaf 2002). The Meyer equation has an upper asymptote of a, but it lacks an inflection point. It can be expanded to the Weibull type function by add ing a positive power parameter to D, and to the Chapman-Richards type function by adding a power parameter to the whole expression in the parentheses. These expansions naturally also have an upper asymptote and the Chapman- Richards type of function also has an inflection point. Both of these expansions have been recommended and used as H-D curves, the former by Huang et al. (1992). Zhang (1997) and Ishii et al. (2000) and the latter by Huang et al. (1992), Zhang (1997) and Zhang et al. (2002). The Korf curve has both an asymptote and an inflection point. The Korf curve has been used both in the form where c= 1 (Curtis 1967, Zakrzewski and Bella 1988, Arabazis and Burkhart 1992, Calama and Montero 2004) and with other positive values of c (Huang et al. 1992, Lynch and Murphy 1995, Lappi 1997, Hökkä 1997, Jayaram and Lappi 2001. Colbert et al. 2002. 111. IV). hi addition to these func tions. the Schnute function (Schnute 1981). which has both an inflection point and an upper asymptote, has been recommended in many studies (Huang et al. 1992. Zhang 1997. Yancai and Parresol 2001). 2.3 Other components The diameter distribution and H-D curve of a stand are regarded as the most important com ponents of stand description in forest manage ment planning. Other important components are. for example, taper curves, spatial pattern and age distribution. It is well known that tree taper varies from stand to stand (Lappi 1986. Ojansuu 1993) and spatial patterns may be very different in different stands (Lin 2003). In the future, information about tree taper may be obtained from harvester measurements by using imputation methods and the spatial pat tern may be determined by using high resolu tion remote sensing imagery (e.g. Uuttera et ai. 1998). However, these approaches are not currently in use, and field measurements of these stand characteristics are too time consum ing in inventories for forest management plan ning. Thus, the taper curves of Laasasenaho (1982) are assumed to apply to all stands and the spatial pattern within stands is assumed to be random. The diameter distribution and H-D curve produce a static description of the stock struc ture, which is the basis of growth and yield prediction. The next important component of the stand description is a set of models that predicts the development of the stand (e.g. Hynynen et ai. 2002). However, because the estimation of stand growth is out of the scope of this study, growth models will not be dis cussed any further. 2.4 Other approaches An alternative to the separate prediction of diameter distribution and H-D curve is the use of a bivariate distribution as the joint distribution of heights and diameters. Ever since Schreuder and Hafley (1977) proposed it, Johnson's SRB distribution has been widely utilized (Hafley and Buford 1985, Knoebel and Burkhart 1991, Siipilehto 1996, Tewari and Gadow 1999, Siipilehto 2000). Zucchini et al. (2001), however, found that a mixture of two bivariate normal distributions, which has more free parameters than the bivariate Johnson's S BB distribution, fitted their Central European beech stand better. Furthermore, the trivariate Johnson's S BBb distribution has been used in the estimation of the joint distribution of height, diameter and age (Schreuder et al. 1982). However, although these approaches are theoretically appealing, their utility is not self evident when compared to the approach based on a univariate diameter distribution and the H- D curve (Knoebel and Burkhart 1991, Siipilehto 1996, 2000). Furthermore, a bivariate joint distribution of heights and diameters is obtained also by using the estimated diameter distribution and error distribution of the H-D curve. 14 Mehtätalo Because diameter measurements are much easier to obtain than height measurements, the current practice is to predict tree height from its diameter. Laser scanning (e.g. Maltamo et al 2004b) and digital photogrammetry of trees from aerial photographs (Korpela 2004). which are promising alternatives to the inventory by compartments in the future, provide quite accu rate height measurements, while diameter measurements are hard to obtain from the air. Thus, the development of these methods to realistic approaches in forest inventories may reverse the roles of height and diameter in the future (Maltamo et al. 2004b). Percentile-based diameter distribution 15 3 Percentile-based diameter distribution This study uses the percentile-based ap proach in the prediction of diameter distribu tions. The percentile-based approach of diame ter distributions was first presented by Borders et al. (1987) and has later been used in Borders and Patterson (1989), Maltamo et al. (2000), Kangas and Maltamo (2000b) and Eerikäinen and Maltamo (2003). The percentile-based diameter distribution was introduced because of its ability to reproduce multi-modal stand tables and the simplicity of the mathematics needed (Borders et al. 1987). More generally, it is much more flexible than the traditional pa rametric distributional families, e.g. Weibull and Johnson's SB , because the implied assump tions about the form of the diameter distribu tion are weak (Maltamo et al. 2000, Kangas and Maltamo 2000b). This study regards the percentile-based dis tribution as a piecewise defined uniform distri bution. The interpretation of the residuals of the percentile models as the horizontal errors of the percentiles is emphasized, and methods for effective use of the error variances in calcula tions are developed (I, 11, V). Furthermore, it is shown that measured sample order statistics can be regarded as measured percentiles, which can be plotted onto the figure of the c.d.f. and used in localizing the predicted percentiles for a given stand (I, V). Finally, exact analytical formulas for relationships between different stand variables are derived without transform ing the distribution to a stand table. The formu las are utilized in adjusting the predicted per centiles in order to ensure compatibility of the stand description (11. V). The next five subsec tions summarize and complete the properties of the percentile-based approach derived in Pa pers I, II and V. 3.1 Distribution function and density Let the diameter distribution of a stand be described with a strictly increasing vector of diameters, d=(dl ,d 2, ■■■,dk )\ corresponding to predefined values, p=(p,,p:,...,pk)'. of the cu mulative distribution function (c.d.f.) where p,=o and pk = 1. Assuming that the c.d.f. be tween consecutive percentiles is linear, the c.d.f. of tree diameter Y is (I): where (The capital letter Y is used for random variable and the lower case letter y for its realization; a notation that is commonly used in statistical literature.) The notation F v (y\A m ) emphasizes that vector p is a predefined constant vector that defines the distribution family used and vector d,„ includes the parameters of the distri bution in stand m. The density is obtained by differentiating (4): 0 V<£?, F„ (>'|d,„)=- a,+b,y d d t b = Sm.—El and a = p.-bd . ' d M -d, 16 Mehtätalo Table 1. Equations for the calculation of some stand characteristics from the percentile-based distribution of form (4). Note. The equations for the basal area weighted diameter distribution are on the left and for the unweighted diameter distribution on the right. The vertical lines in the equations of median mean the value of d that satisfied the equation on the right hand side of the vertical line. The notation /' means the density of the weighted diame ter distribution and/ 1 the density of the unweighted distribution. 3.2 Moments of the distribution The expectation or mean of tree diameter in a stand is calculated as the expectation of the percentile-based diameter distribution (Eqs 4 and 5). It follows straightforwardly from the definition of the expected value that (Casella and Berger 2002. p. 55) For calculating variance, the second moment of the distribution, defined as E(Y 2 ). is needed (Casella and Berger 2002. p. 59). It follows again from the definition of the expected value that and the variance is obtained using the well knovvn formula (Casella and Berger 2002, p. 60) dk £(>')= J.i/po'|d,„>/v • (6) =Z ] >'b.dy = z;Y. b,{dJ - d?) ,=l d £ i=l E{Y 1 )= }v 2/ P (v|d„,Xv (7) = Z jy^dy^^idj-d,3 ) ,=\ di I=l weighted unweighted k ] -ln<) }2>,K, 2 -4 2 ) • 1 1 i- =1 unweighted mean k-\ I '\ A £ z i=i weighted mean 1 2 &K.2 i=i i -<-) 3 2», 4 &(4.,' - d;) 1=1 unweighted median d l dsovo Z", i=i ' l U" 1 li ' 2 2 4 J weighted median ds 0% d~ )rf(,)dt=\ X b\d,j-d?y 2 i=i I 1 1 Zt>,(dj-J,i) quadratic mean J k-1 | 2>' i=l 1 r 1 1 1 U" d,+\ y 1 3 stem number/basal area K4 b, h n < 1 ,7" 1 ] 1 / 12 Percentile-based diameter distribution 17 3.3 Relationships between percentiles and stand characteristics Many computations with the theoretical dis tributions used in the description of forest structure are very complicated and lead to results that do not exist in closed form. A usual solution to this problem is to use a stand table approximation of the distribution (e.g. Nepal and Somers 1992). With the percentile-based diameter distribution (Eqs 4 and 5), many stand characteristics can be rather simply derived analytically (See the Appendix of Paper II). The density of Y r:„ is The most important stand characteristics are given in Table 1 for both basal-area weighted and unweighted diameter distributions. 3.4 Order statistics The theory of order statistics is of great im portance with the percentile-based diameter distribution because a measured sample order statistic is an unbiased estimator of a certain percentile of the underlying distribution (I). Denote the r lh order statistic in a sample of size nby Yr:„. If the sample is drawn from a distribution with a c.d.f. of the form (4). the exact distributions of order statistics can be derived rather easily. and the joint density of two order statistics Yr]:„ and Yr2:„ is where (Reiss 1989.1). The expectation and variance of a single or der statistic follow from the same definitions and general rule that were used in Equations 6, 7 and 8 (Casella and Berger 2002. p. 55-60. 1). The covariance of two order statistics can be calculated using the corresponding definitions and rule (Casella and Berger 2002. p. 144. 170) for bivariate distribution, as shown in Paper I. Furthermore, the measured sample order statis tic is an unbiased estimator of the 100/? lh per centile of the diameter distribution of the stand, where and Fp is the c.d.f. of the stand. 3.5 Considerations on the PPM with the percentile-based approach Assume that the percentiles of the diameter distribution in stand m follow the model var(T) = £(}' 2 )-[£(y)] 2 . (8) f , lA = r ~'(l-a,-blyy''ifd i ys~\aj +b>y2-a - b=¥]. var(x2)=V2, and cov(xi,x2') = Vi - Using the notation of McCulloch and Searle (2000, p. 247), this can be written as Assume that we have observed random vec tor x 2 and we want to predict vector X! The best predictor (BP) of X! is the conditional expecta tion (McCulloch and Searle 2001. p. 248). The best predictor can usually not be calculated because it requires the distribution of x,|x 2. An estima tor that requires only first and second moments is obtained by limiting the consideration to linear unbiased predictors. The Best Linear Unbiased Predictor (BLUP) of Xj is with the prediction variance of (McCulloch and Searle 2001. p. 250). Thus, if the expectations and variance-covariance ma trices of two random vectors are known and either of them is observed, the other can be predicted. The variance of the prediction error can be calculated using Equation (20). If x follows the multinormal distribution. BLUP is also BP. The standard theory of mixed models utilizes BLUP to predict the realizations of the random effects in the modeling data (McCulloch and Searle 2001, p. 254-258). It is a special case of prediction, where the correlations between the observed random vector. x 2=y-Xp. and the vector of random effects, x,=b. are generated by the structure of the data (Eq. 15), i.e. V 2=ZDZ'+R and Vl2=DZ\ If measurements of the response variable from a new stand are available in applications, the realizations of the random parameters in that stand can be pre dicted (Lappi and Bailey 1998. Lappi 1991. III). Similarly, if we have several models that have correlated residual errors and the correla tions are known, an observed response of any single model can be utilized in predicting the responses of the other models in that stand (Lappi 1991, I). Thus, the realized random effects and/or random errors of the models can be predicted in order to localize the models for that stand. In practice, the variance-covariance matrices used are replaced with their estimates. A pre dictor obtained using these estimates is some times called EBLUP, Estimated Best Linear Unbiased Predictor (e.g. McCulloch and Searle 2001. p. 257). 5.3 Constrained optimization (11, V) This study utilizes constrained optimization in adjusting the percentiles and stand variables to obtain a compatible stand description (11. V) and in searching for an optimal measurement strategy for a single stand (V). A general con strained optimization problem is defined as follows (Bazaraa and Shetty 1979. p. 2): minimize subject to where z=(: l - 2, :„)' is a vector of decision variables and /, g;,...,g„, and h, /?, are func tions of z. The function / is called an objective function, functions g/ g,„ inequality con straints and functions hi,...,hi equality con straints. The solution to the above problem is a value of z that minimizes the function / and x = M, (16) x 2 M rrAjjv, v u Ti x : V u * V,J BP(\ t )= £(x,|x,) (18) BLUP(x,) =i,= n, + V; V l2 "'(x : - ) (19) var(i,-x t)=y-V; : V,-Iy2' 1 y 2 ' (20) /(*) (21a) g,(z)•)]"' [l - Fr (y)]"" , (2) frn (}')=Po(r - n )b, [a, +b,y\ 1 [l-tf, -&, v]" r , (3) dt E( Yr„) =Ä, (''•") \yfr ,\y)dy ■ £( }'„) = A('"-")Z \yfr,:(yHv ■ (4) '= ' d, P = Ft [E(Yrn )]. (5) E{Yj)=/)o (r,n)£ jy 2 fr ,:(y)dy '= 1 d, var(7 r „) = £(7j)-[£(f r„)] 2 . (6) fr, „r,r(V, , V, ) = Mri* r2 -n)fr ( V,)/, (V 2 ) X (7) [wrw^rNwr n ' PM- r i- n ) = -, w rr: • { n ~ri)-\ r 2 ~r \ -l)!(ri "I)! ,!' J {>\-y2)=Mrr'2-n)b,bx (g) (a{ +bty[ )' \a j +b J y 2 -at -b iy^ 2 '' (l-a J -b j v 2 )' Localizing predicted diameter distribution 5 gives Using this result and the expectation of single order statistic, (4), the covariance of two order statistics is calculated as Since the densities (3) and (8) are n- Ith1 th and n-2 th order polynomials, respectively, exact calculation of the integrals may lead to degen eration of computing accuracy with large val ues of n and need to be calculated numerically in applications. Predicting localized percentiles Assume that models predicting k diameter percentiles are available and the 100*/?, 01 per centile in stand m follows the model d m = fijm +em for i=1,2,...,k. Writing all per centiles as a vector, the model for the whole set of percentiles in stand m is where d, „={dlm, d2m d k,„)\ Hkm y and e m=ielm e2m ...,ekm y with E(e,„)=o. Assuming that the model is correct and the parameters are known, vector fim includes the conditional expectations of the percentiles given the values of the stand variables, £(d|x,„). These are later referred to as expected percen tiles of the stand. Vector em includes the stand effects of stand m, i.e. the deviations of the stand-level percentiles from their conditional expectations. Estimation of the model (Equa tion 10) using a simultaneous regression tech nique (SUR) produces an estimate of the vari ance-covariance matrix of stand-level devia tions, var(e,„), denoted by V)k/k . Writing the predefined p-values into vector p=(pi,p2,...,Pk)'■ the diameter distribution based on expected percentiles of stand m can be writ ten as Perc p (H,„) and the diameter distribution of stand m, correspondingly, as Percv(\i,„ +e,„). Since the following calculations apply to one stand, the stand index m is dropped hereafter. The measured order statistics are used to predict the vector of stand effects, e* x i., using the standard linear prediction theory (see e.g. Lappi 1986. 1991. 1997). Assume that we have q measured order statistics from a stand in vector d*=(Yr: .„2,..., Yrq:mi )'. The meas urements follow the model where the measured diameters are in vector d* r/ xi, their conditional expectations in vector the stand effects in vector e* r// ] and the sampling error in vector For the random part of the model E(e*)=E(£)=o and cov(e*, e')=o holds. At this stage we assume that the actual distribution of the stand, Perc v {\i +e), is known (I return to this later). The /rvalues corresponding to the measure ments are obtained using formulas (4) and (5) and written into vector p* 9 xi. The conditional expectations of measurements are then calcu lated as n*=F'(p*), where F 1 is the inverse of the distribution based on expected percentiles, Perc p(\i). For predicting the realized stand effects of model 10. the sampling variance-covariance matrix of order statistics, var(e), denoted by Rqxq, the variance-covariance matrix of the stand effects of model 11, var(e*), denoted by D * qxci and covariance matrix of stand effects of models 10 and 11, cov(e.e*'), denoted by Ck , q are needed. The variances of order statistics are calculated using formula (6) and written on the diagonal of R. If one has measured several order statistics from the same sample plot, their covariances are calculated using formula (9) and written in the corresponding cells of matrix R. The covariances of order statistics from different plots are naturally zero. Matrices D* and C are obtained from matrix D by interpo lating it for the values of p*. In localization, the unobserved random vec tor e is predicted using the observed random vector e*+£=d*-fi*. For the prediction, we [t/„ dl=i)*[dJ,dj+l ) separately for /.y-1 1 P\ (l -")ZZ I j . V'l-V;/ij ~.r: „' J (V,, v2Kv2^, 7=l '=' d, d. cov(}; „)=£(>; „)-£(>; „)E(J; „). (9) d = u +e . (10) /»/ * »/ m * v / d* = n*+e*+s, (11) 6 Mehtätalo Figure 1. The distributions based on the expected percentiles, µ (dashed line) and on the localized percentiles, µ+e (solid line) obtained using the ob servation Y3:13=8.5 cm (￿). need to derive the variances and covariances of these random vectors using the known matrices D, D*. C and R. The variance of stand effects of model 10, var(e), is straightforwardly D. Since cov(e*, e)=o and cov(e, e)=o, var(e*+e)=D*+R and cov[e.(e*+£)']=cov(e.e*')=C. Using the nota tion of McCulloch and Searle (2001, p. 247), these results can be expressed as The Best Linear Unbiased Predictor (BLUP) of e is calculated as with the prediction variance of (McCulloch and Searle 2001, p. 250). The stand level diameter percentiles are obtained by adding the predicted stand effects to the ex pected percentiles (see model 10). The stand-level diameter distribution would already be needed in the calculation of p* and R. Since it is the result of the prediction and is not known when p* and R are calculated, the solution is searched iteratively. At the first iteration step, the expected diameter distribu tion of the stand, Perc v(\x), is used as the stand level distribution to approximate p* and R. Subsequently, these approximations are used to predict the stand-level diameter distribution, Perc p (fi+e), with BLUP. Furthermore, this prediction is used to calculate new approxima tions of p* and R and the prediction of stand effects is carried out again. Repeating this until the predicted stand-level percentiles converge gives the final predicted stand-level diameter distribution. In some cases, we may obtain two observa tions of the same percentile. This happens when there are two measurements of the same percentile from different plots, i.e. the number of tally trees is the same on two sample plots from the same stand and diameters with the same rank are measured from both plots. This means that the same row (and column) is in cluded in matrix D* twice. It is therefore singu lar and the calibration cannot be carried out. The non-singularity of matrix D* can be guar anteed by treating the several measurements of the same order statistic as one measurement in the calculations. In this case, the mean of measured diameters is used as the observed percentile and the elements of matrix R are calculated using general rules for the variances and covariances of sums. Numerical example This section presents a numerical example of the calibration algorithm. The model of Kangas and Maltamo (2000 a) was used to predict the conditional expectations of A=ll percentiles of basal-area diameter distribution. The predictors of the model are basal area median diameter (DGM), age, basal area and soil type. The models predict the oth0 th 10 th , 20 th , 30 th , 40 th , 60"\ 70 th . 80lh , 90 lh , 95 th and 100th percentiles, i.e. In addition, the known DGM is used as the 50 th percentile. The models have been estimated using seemingly unrelated regression and the estimated variance-covariance matrix D was "el r To") Fd C ?| D*+R_JJ e = C(D * +R) ' (d * -(i *) (12) var(e-e) = D-C(D*+R) ' C' (13) p=(o 0.1 0.2 0.3 0.4 0.6 0.7 0.8 0.9 0.95 1.0)'. B Localizing predicted diameter distribution 7 obtained from the original SUR-fit (Table 1). The model of percentiles is estimated on a logarithmic scale, i.e. the model is of the form (cf. Equation 10) In order to give a numerical example of the proposed algorithm, the diameter distribution of a Norway spruce stand was predicted, where DGM is 20 cm, basal area is 22 m 2 /ha , stand age is 64 years and site fertility class is mesic. The predicted logarithmic percentiles using these values (i.e. conditional expectations given the values of stand variables) were In addition to the known stand variables, the third smallest tree of a horizontal point sample of 13 trees has been observed to be 8.5 cm in diameter, i.e. Y 3 ]3 has been observed to be 1n(8.5)=2.140. The expectation of the sample order statistic Y3:13 of the distribution Perc v (\i) is (Equation 4) E(Yi u )= 2.595. Equation 5 gives the value p*= 0.182, which means that the measured order statistic is a measurement of the 18.2 th percentile of the diameter distribu tion. The observed percentile has been plotted onto Figure 1 at the location (2.140,0.182). Note that we have £=ll percentiles of prede fined percentage values and <7=l measured percentiles. Thus, p*. e*, e, R and D* are sca lars and C is a vector with a length of 11. The next step is to predict the stand effects of the percentiles, e (Equation 14). In order to be able to do it, we need the variances of sampling errors (R), the variances of the stand effects (D*) and the covariances between the stand effects of model 14 and the stand effect of the 18.2 th percentile (C). The variance of sampling error is (Equation 6) R=var(e)= 0.0577 and the variance of stand effect is obtained by interpo lation of matrix D (Table 1) for the value of 0.182. Linear interpolation gives D*=var(e*)= 0.0746+(0.0293-0.0746)/0.10*(0.182-0.1)=0.0376. The covariances (C) are also obtained by linear interpolation of matrix D (Table 1) as The stand effects of model 14 are predicted as (Equation 12) The localized percentiles are obtained by adding the stand effects to their conditional expectations (see Figure 1). The next step in the calculations would be to calculate the values of p* and R again by as suming that the sample has been drawn from the localized distribution and to iterate this until convergence. However, in this example the iteration is not carried out. ln(d) = n + e. (14) f 1.63 2.35 2.65 2.79 2.91 H= 3.08 . 3.15 3.24 3.31 3.38 13.53J " 0.050 1+(0.0223-0.0501)/0. 10*(0.182-0.1) "l f 0.0274 N c= ; : -0.00903+(-0.00774+0.00903 )/0.10*(0.182-0.1) j [-0.00798, ' 0.0274 f -0.13 T e = ; (0.0376 + 0.0577) ' (0.0274 -0.00798)( 2.140 - 2.595) = : -0.00798 J [0.0381 y 8 Mehtätalo Table 1. The within-stand variance covariance matrix (D) of the percentile models of Kangas and Maltamo (2000 a) for Norway Spruce (Kangas. A., personal communication). Testing with real data Arrangement of the test The calibration was tested in a small dataset consisting of 43 fixed rectangular sample plots from mixed Scots pine-Norway spruce stands. For the test, Norway spruce trees were se lected, because the number of Norway spruce sample trees per stand is much greater (42-222) in the data than the number of Scots pine trees (13-168). Furthermore, the Norway spruce data is more challenging than the Scots pine data, because it includes various forms of the diame ter distribution, including symmetric, skewed, extremely wide, mound-shaped, bimodal and multi-modal forms. The dataset has been originally collected by Pukkala et al. (1994) for productivity studies and it has been further used by Kangas and Maltamo (2000b) and Mehtätalo (2004) in testing the performance of diameter distribu tion prediction algorithms. The size of the sample plots in the data varied from 600 to 3000 nr. the number of tally trees from 42 to 222, the basal area of Norway spruces from 1.54 to 24.07 m 2 /ha and DGM fVom 5.5 to 33.9 cm. Because of the fairly small number of trees in each plot, it was assumed that a slightly smoothed distribution describes the distribution of the stand better than the actual measured one (Droessler and Burk 1989, Maltamo and Kan gas 1998). Hence, the actual distribution was smoothed with a Gaussian kernel (Härdle 1990, p. 15-20) using bandwidth determined by the function where wm is the width of the actual distribution and Nm is the number of sample trees in stand m (see Mehtätalo 2004). The smoothed tree stock of the original sample plot is subse quently referred to as a stand. To illustrate the effect of calibration on the accuracy of the predicted diameter distribution of the stand, a varying number of horizontal point samples were simulated in each stand and two sample trees were randomly selected from each plot as sample measurements of sample order statistics. The number of sample plots in a stand was varied systematically from 1 to 6. thus resulting in the number of sample trees from any one stand being 2, 4, 6. 8. 10 or 12. The trees belonging to a sample plot were se lected by generating a uniformly (0.1) distrib uted random number for each tree of the stand and selecting those trees for which the random number was less than the sampling probability of the tree. The localized predictions of diame ter distributions were calculated by predicting the stand effects of model 14 using these meas urements. The models of Kangas and Maltamo (2000 a) were used in the test. In this test the exact val i u'„i 6.5 d, d : ds d< d5 d6 d? ds d, d/o du d, 0.162 d: 0.0501 0.0746 d 3 0.0223 0.0349 0.0293 (symm) d4 0.0107 0.0156 0.015 0.0142 ds 0.00689 0.00877 0.00935 0.00933 0.0098 d 6 0.00021 -0.00269 -0.00265 -0.00153 -0.00093 0.00319 d7 -0.00274 -0.0051 -0.00395 -0.00245 -0.00133 0.00301 0.00592 ds -0.00548 -0.00729 -0.00592 -0.00357 -0.00249 0.00296 0.00584 0.00868 d, -0.00655 -0.00814 -0.00656 -0.00447 -0.00328 0.00311 0.00579 0.00835 0.011 dm -0.00672 -0.00818 -0.00699 -0.00479 -0.00351 0.00305 0.00597 0.00837 0.011 0.0135 du -0.00982 -0.00903 -0.00774 -0.00541 -0.00366 0.00281 0.00625 0.00863 0.0109 0.0138 0,025 Localizing predicted diameter distribution 9 ues of the stand characteristics were used in predicting the expected percentiles to guarantee an equally accurate predictions regardless of the number of sample plots. To ensure safe interpolation of matrix D, linear interpolation was used; and to guarantee logical behavior of the interpolated variances and covariances, the interpolation of variances was carried out be fore interpolating the covariances (see the nu merical example presented above). Since the models were originally estimated on the loga rithmic scale, the calibration was carried out on the same scale. Thus, the logarithmic percen tiles and diameter measurements were used instead of the arithmetic ones in all calcula tions. When transforming the logarithmic per centiles onto the arithmetic scale, the bias caused by the logarithm transformation was corrected by adding half of the prediction vari ance to the predicted logarithmic percentiles before applying the exponential function (see Lappi 1991). The correction terms were ob tained for the expected percentiles from the diagonal of D and for the localized percentiles from the diagonal of var(e-e) (Equation 13). In some stands, the calibration algorithm failed for two reasons. Firstly, the iteration did not converge before the maximum number of iterations was reached and secondly, there was some iteration step after which the distribution was not monotone and thus the iteration could not be continued. These situations were cir cumvented in the test by repeating the sam pling until the calibration did not fail. How ever. to get a reliable figure about the perform ance of the method in practice, the proportion of those cases where the first attempt failed for one of the first two reasons was calculated. To compare the predicted distributions with the actual one. bias and relative root mean square error of volume (in per cents) were calculated. The volume was calculated in 1 cm classes, using the volume functions of Laasasenaho (1982) with tree DBH as the pre dictor. In addition, an error index proposed by Reynolds et al. (1988) was used in the com parisons to measure the goodness of fit of the distributions. The error index was calculated in 2 cm classes using the basal area as weight. Thus, the error index was the sum of the abso lute differences between the actual and pre dicted basal areas of the diameter classes. The calculation was repeated 100 times to decrease the effect of sampling error in the results. The calculations were carried out in R, an environment for statistical computing (www.r project.org, Venables and Ripley 2002). In addition, the numerical integrations needed in the calculation of expectations, variances and covariances of order statistics were carried out with IMSL-subroutines(lMSL 1997), which were linked into R. Test results The effect of calibration on the accuracy of the predicted diameter distribution was consid erable. Even two sample trees per stand de creased the RMSE of volume from 2.71t0 2.52 and the absolute bias of volume from 0.84 to 0.66. The mean of the error index also de creased clearly, from 4.67 to 4.37. Further more. increasing the number of sample trees improved the accuracy of the localized distri butions steadily (Figure 2) with the result that 12 sample trees reduced the RMSE of volume by 29 per cent and the absolute bias by 62 per cent when compared to the distributions based on expected percentiles. The prediction results were relatively accurate because the exact stand variables were used in the prediction of expected percentiles and all calculations of volume were based only on tree diameter, not on tree height. To demonstrate the effect of iteration. Figure 2 presents the results for localized distributions after the first iteration step and converged itera tion. The iteration had a considerable effect on the accuracy of volume prediction. With two sample trees, the accuracy of volume predic tion after the first iteration step was even lower than that of the distributions based on expected percentiles, but with a larger number of sample trees, the first iteration step also improved the accuracy. However, regardless of the number of sample trees, the advantage of iteration for the accuracy of volume predictions was re markable. The effect of iteration on the error index, on the other hand, was somewhat con fusing. With two sample trees, the values of the error indices were almost the same after the 10 Mehtätalo first iteration step and converged iteration. With a large number of sample trees, on the other hand, the iteration increased the value of the error index when compared to the value after the first iteration step. I will return to this. Visual examination of the localized distribu tions showed that they are realistic (Figure 3). In most cases the localized distributions were more consistent with the form of the actual distribution than the distributions based on expected percentiles were. Furthermore, the calibration algorithm could produce skewed and multimodal distributions in cases where the distribution based on conditional expectations was rather symmetric and unimodal (Figure 3. stands a and b). In addition, calibration usually moved the predicted minimum and maximum diameters in the right direction. The prediction errors of the localized percentiles were smaller than the errors of the expected percentiles (Figure 4). In some cases the short distance between consecutive percentiles caused quite high peaks in some diameter classes (e.g. stand a in Figure 3). However, the frequencies of the diameter classes beside the peaks were usually low, which compensates for the error caused by peaks. Hence, peaks do not increase the RMSE and absolute bias of stand volume remarkably but they may have an effect on the error index. Visual examination of the localized distribu tions showed that the peaks were in many cases higher after converged iteration than after the first iteration step. In addition, increasing the number of sample trees increased the occur rence of high peaks. This may be an explana tion for the result that if there were a large number of sample trees the error indices after the first iteration step were lower than those after converged iterations. The peaks are in many cases in percentile intervals just below or above the 50 th percentile (=DGM) and may be a consequence of using the measured DGM as the 50 lh percentile. A peak arises when the measured DGM is far from the midpoint be tween the localized 40 th and 60 th percentile. For convergence it was required that the sum of the absolute differences of the percentiles at subsequent iterations was less than 10" 6 . The Figure 2. The effect of increasing the number of sample trees on the relative RMSE (a) and bias (b) of volume and the basal area weighted error index (c). The results are presented both using the predic tions after the first iteration step (￿) and after con verged iterations (•). maximum number of iterations was 50. A sam ple with which the calibration was successful could easily be found for all stands. In the following, the proportions of failed calibrations are presented for the first sample attempted. The number of sample trees had a clear effect Localizing predicted diameter distribution 11 Figure 3. Examples of the predicted and localized distributions with 2. 6 and 10 measured sample trees from 1. 3 and 5 sample plots, respectively. The upper graphs are the cumulative distributions and the lower ones the densities. The dotted lines are the distributions based on expected percentiles, dashed lines the localized distri butions after the first iteration step and solid lines the localized distributions after converged iteration. The actual distribution of the stand is the stepped line in each of the upper graphs and the histogram in the lower ones, hi the upper graphs, the marks are the measurements (the same symbols are used for measurements from the same sample plot). on the convergence: with two sample trees per stand the iteration converged in 98% of the stands and the average number of iterations in these stands was 8.28, while with 12 sample trees the proportion of converged iterations was only 88% with an average number of iterations of 12.0. Of the unsuccessfully localized distri butions. the iteration did not converge in 7%, and the percentiles were not monotone in some iteration step in the other 93%. The above calculations were carried out with the Norway spruces of the test data. Because Norway spruce is a shade tolerant tree species, its diameter distribution may have various forms. The same calculations were carried out also with the Scots pines of the test stands. Scots pine is a shade intolerant tree species, which usually has a unimodal diameter distri bution. Thus, the diameter distribution of Scots pine is much easier to predict than that of Nor way spruce. With Scots pines, RMSE and ab solute bias of volume were lower than in the Norway spruce data when using the expected percentiles in the prediction. In addition, the effect of calibration on the RMSE and absolute bias were clearly stronger. The error index was, however, greater in the Scots pine data than in the Norway spruce data, but the effect of cali bration was. again, stronger. For example, in the Scots pine data, two sample trees decreased the RMSE of volume from 1.31 to 1.10, abso lute bias from 1.03 to 0.63 and mean of the error index from 5.76 to 5.41. Thus, the cali bration seems to be clearly useful in both Scots pine and Norway spruce data. 12 Mehtätalo Figure 4. The effect of calibration on the accuracy of the predicted percentiles in stand b of Figure 3. The dotted line is the distribution based on expected percentiles of the stand, the solid line is the localized distribution, and the marks are the sample measure ments. The lengths of the horizontal lines at the measurements and predicted percentiles show the standard deviation of the sampling and prediction error, respectively. The standard deviations are obtained from the diagonal elements of the matrices D, var(ê-e) (eqn 11) and R (eqn 6) for the ex pected. localized and measured percentiles, respec tively. Discussion This study presented an algorithm for com bining the information from sample order sta tistics and the predictions of expected diameter percentiles. The test results showed that using this method the accuracy of predicted diameter distribution could be improved considerably with a small number of sample tree measure ments. Furthermore, by increasing the number of sample trees, the predicted distribution seemed to approach the actual distribution of the stand. Because of its flexibility, the percentile based diameter distribution has been found to be a good alternative in stands with complex diameter distributions (Maltamo et al. 2000). This study showed that using measured order statistics further utilizes the flexibility of the percentile-based method. Thus, the proposed method seems to be a promising alternative for the prediction of diameter distributions in com plex stands. A clear advantage of the method is its flexi bility in producing predictions of different accuracy. In other words, using this method, a realistic prediction of the diameter distribution can be achieved with quite a small number of measurements (i.e. even with no sample order statistics measured) and the accuracy of the predicted distribution can be improved steadily by increasing the time used in measuring sam ple trees. Hence, any realistic predefined accu racy requirement can be obtained with a suffi cient number of diameter measurements and. on the other hand, any allotted measurement time can be used effectively for improving the accuracy of prediction. Therefore the method can be useful in various inventories with dif ferent accuracy requirements: for example, for tactical forest planning one or two sample trees may produce a sufficiently accurate prediction, while for the pre-harvest measurement of a stand some more measurements may be needed. The aim of this study was to explain the methodology, demonstrate its use and show how the proposed method works with accurate measurements of order statistics from i.i.d. samples. No comparisons with other methods were carried out, since no methods that utilize the same input information are available. How ever, the models used in the prediction of ex pected percentiles have been thoroughly tested (Kangas and Maltamo 2000b). They were found superior when compared to the parame ter prediction methods based on the Weibull distribution. Furthermore, the percentile mod els of Scots pine were found to produce ap proximately as accurate predictions as the k nearest neighbor method of Maltamo and Kan gas (1998). This study showed that predictions based on expected percentiles can further be remarkably improved by the use of measured order statistics. When deriving the expectations, variances and covariances of order statistics, the sample was assumed to be independent and identically distributed (i.i.d.), i.e. the stand was assumed to be spatially homogeneous. The method was tested with simulated i.i.d. samples, where selecting sample trees was not based on the location of the tree in the stand. This approach Localizing predicted diameter distribution 13 was selected to test the method in a situation where the assumptions behind the method are valid. In practice, however, stand characteris tics are spatially correlated within a stand. Thus, the test results of this study are an exam ple of localization in a somewhat ideal situa tion. In a practical situation, where samples are not i.i.d, the advantage achieved by calibration may be smaller than in the test situation of this study. In addition, it may be profitable to dis tribute the measurements to several plots in a non-homogenous stand to get a spatially more representative sample. A more reliable view of the advantage of the method would be achieved by testing the method in a dataset with real measured sample plots. Figure 3 showed that the localized distribu tions may be of various forms, including bi modal and skewed forms. However, for exam ple in stand c of Figure 3, all combinations of sample order statistics did not produce a bi modal distribution. Thus, attention should be paid to the selection of measured order statis tics. For example, in a stand where distribution seems to be bimodal, the measurements of maximum diameter of the lower peak and minimum diameter of the upper peak would probably provide more information about the form of the distribution than two randomly selected sample order statistics. Furthermore, measurement of the smallest possible sawtim ber tree as a sample order statistic would probably lead to more accurate estimation of sawtimber volume than measurement of a ran domly selected tree. Thus, in the future, the usefulness of different strategies in the selec tion of sample order statistics should be stud ied. Since the percentile models used predict the diameter distribution weighted by basal area, the sample plots of this study were horizontal point samples. However, the approach pre sented can be used as such with an unweighted diameter distribution. In this case, fixed-radius sample plots should be used instead of horizon tal point sample plots Since the measured order statistics are re garded as measured percentiles, the approach of this study requires the percentile-based di ameter distribution. However, the approach could be generalized also to other distribution families, when the parameter prediction method is used. This requires derivation of the expectations, variances and covariances of order statistics using the distribution family used. Furthermore, the covariances between the measured percentiles and the parameters of the distribution family should be known. Estimates of them would be obtained by fitting simulta neously models for the parameters of the dis tribution family and percentiles of the distribu tion. The required covariances of the measured percentiles and distribution parameters would be obtained by interpolation of the estimated within-stand variance covariance matrix for the />values of the measurements. No measurement error was assumed for the sample plot measurements. In practice, how ever, measurements always include error. When measuring sample order statistics, the error may be either a measurement error of diameters or an error in the determination of the rank of the tree on the plot. The measure ment error of diameter is quite easy to take into account by adding the variance of measurement error on the diagonal of matrix R (see Lappi 1991). The error in the determination of rank, on the other hand, is somewhat more compli cated. It can be regarded as a measurement error of diameter, the magnitude of which is related to the difference between the diameters of subsequent trees in the sample plot. Hence, the measurement error of rank should increase the diagonal elements of R proportionally to the absolute difference between subsequent diameters on the sample plot. However, most errors in determining the rank presumably happen with trees that have diameters close to each other. In this case the effect of measure ment error is quite small and should not cause serious problems. In practice, the predictors of percentiles in clude measurement or sampling error, which may be large. Therefore, in addition to the error caused by the between stand variation, the predictions of expected percentiles also include error variation and bias caused by the meas urement error of the predictors. Taking this error into account would require derivation of prediction error variances and covariances of 14 Mehtätalo the expected percentiles as well as some kind of bias correction on the expected percentiles. Since the aim of this study was merely to in troduce the concept of combining the two types of studies, the work required for taking such errors into account was not tackled. A problem with the proposed method was that the localized distribution sometimes had peaks around DGM. A simple solution to the peak problem would be to use linear interpola tion between the localized 40 th and 60 th percen tiles instead of using DGM as the localized 50 th percentile. Another, more elegant solution would be to use the measurement error vari ance of the DGM as the within-stand variance of the 50 th percentile and derive its effect on matrix D. Thus, taking the effect of measure ment errors of stand variables into account would probably solve the peak problem. If a practical situation had been simulated, the expected percentiles should have been predicted using stand variables calculated from the sample plots. However, it would have af fected the accuracy of the expected percentiles since the number of sample plots varied in different simulations. To make the test show only the effect of calibration on the accuracy of the predicted distribution the exact stand char acteristics were used instead of the means of the sample plots. The procedure used with the unsuccessfully localized distributions was to pick a new sam ple from the stand and try again, in the test situation of this study. In applications, a proce dure for the situations where the calibration fails is needed. A simple approach would be to drop one or more sample trees from among the measured sample trees and to try to find a solu tion by using the rest of the sample trees. How ever. this would waste expensive field informa tion. Because the main reason for failure was the non-monotony of the percentiles, the pro portion of failures could be decreased remarka bly by ensuring the monotony through formula tion of the percentile models. The calculation of exact expectations, vari ances and covariances of order statistics with numerical integration was rather slow, espe cially when the number of measured sample order statistics was high. For example, the calculation of the 4300 sample plots with two sample trees took 1.5 hours with a modern personal computer and increasing the number of sample trees to 12 increased the time re quirement to almost 4 hours. In these calcula tions. IMSL subroutines DCONG and DQDAG (IMSL 1997) were used. With a view to speed ing up the computations, the calibration was also carried out using asymptotical expecta tions and variances of sample order statistics (Reiss 1989, p. 109-110), integrals of which can be calculated analytically. With this method the accuracy of the localized distribu tions was clearly better than that of the distri butions based on expected percentiles, but nevertheless clearly worse than using exact results. This could be expected, since the asymptotical results are quite far from the exact ones with such a sample size that is obtained using relascope factor 1. The time requirement could probably be decreased by adjusting the parameters of the numerical integration algo rithm. Furthermore, more effective algorithms might be available and the development of faster computers will eliminate this problem in the future. However, even though the results of this study represent a somewhat ideal situation and there are many things that could be taken into account in order to improve the method, it is a promising method that brings together the good properties of the parameter prediction method and of the use of sample information. Acknowledgements This work was part of the project "Statistical modeling for forest management planning", which was carried out at the Finnish Forest Research Institute and funded by the Academy of Finland (decision number 73392). I would like to thank Dr Juha Lappi for the clarifying discussions we had during this work. In addi tion. I am grateful to Professor Annika Kangas, Professor Matti Maltamo, Dr Pekka Leskinen, the two anonymous referees and the associate editor for their valuable comments on the manuscript. Finally. I wish to thank Dr Lisa Lena Opas-Hänninen for revising my English. Localizing predicted diameter distribution 15 Literature cited Borders, 8.E., Souter, R.A., Bailey, R.L. and Ware. K.D. 1987. Percentile-Based Distribu tions Characterize Forest Stand Tables. For. Sci. 33(2): 570-576. Borders, B.E. and Patterson. W.D. 1990. Pro jecting Stand Tables: A Comparison of the Weibull Diameter Distribution Method, a Percentile-Based Projection Method, and a Basal Area Growth Projection Method. For. Sci. 36(2): 413-424. Burk, T.E. and Newberry, J.D. 1984. A Simple Algorithm for Moment-Based Recovery of Weibull Distribution Parameters. For. Sci. 30(2): 329-332. Cao, Q.V. and Burkhart, H.E. 1984. A Seg mented Distribution Approach for Modeling Diameter Frequency Data. For. Sci. 30(1): 129-137. Casella, G. and Berger, R.L. 2002. Statistical Inference, Second Edition. Duxbury Ad vanced Series, Pacific Grove. 660 p. Droessler, T.D. and Burk, T.E. 1989. A Test of Nonparametric Smoothing of Diameter Dis tributions. Scand. J. For. Res. 4: 407-415. Hyink, D.M. and Moser, J.W. Jr. 1983. A Gen eralized Framework for Projecting Forest Yield and Stand Structure Using Diameter Distributions. For. Sci. 29(1): 85-95. Härdle, W. 1990. Smoothing Techniques with Implementation in S. Springer Series in Sta tistics. Springer-Verlag, New York. 256 p. IMSL 1997. Fortran Subroutines for Mathe matical Applications. Math/Library. Vol umes 1&2. Visual Numerics. 1218 p. + Ap pendices. Kangas, A. and Maltamo, M. 2000 a. Percen tile-Based Basal Area Diameter Distribution Models for Scots Pine, Norway Spruce and Birch Species. Silva Fenn. 34(4): 371-380. Kangas, A. and Maltamo, M. 2000b. Perform ance of Percentile Based Diameter Distribu tion Prediction and Weibull Method in Inde pendent Data Sets. Silva Fenn. 34(4): 381- 398. Kilkki, P. and Päivinen, R. 1986. Weibull function in the estimation of the basal area dbh-distribution. Silva Fenn. 20(2): 149-156. Laasasenaho, J. 1982. Taper curve and volume functions for pine, spruce and birch. Comm. Inst. For. Fenn. 108. Finnish Forest Research Institute, Helsinki, Finland. 74 p. Lappi, J. 1986. Mixed linear models for analyz ing and predicting stem form variation of Scots pine. Comm. Inst. For. Fenn. 134: 1- 69. Lappi, J. 1991. Calibration of Height and Vol ume Equations with Random Parameters. For. Sci. 37(3): 781-801. Lappi, J. 1997. A Longitudinal Analysis of Height/Diameter Curves. For. Sci. 43(4): 555-570. Lindsay, S.R. Wood, G.R. and Woollons, R.C. 1996. Stand table modelling through the Weibull distribution and usage of skewness information. For. Ecol. and Manag. 81: 19- 23. Liu, C.. Zhang, L., Davis. C.J., Solomon, D.S. and Gove, J.H. 2002. A Finite Mixture Model for Characterizing the Diameter Dis tributions of Mixed-Species Forest Stands. For. Sci. 48(4): 653-661. Maltamo, M. and Kangas, A. 1998. Methods based on k-nearest neighbor regression in the prediction of basal area diameter distribution. Can. J. For. Res. 28: 1107-1115. Maltamo, M., Kangas, A., Uuttera, J., Torniai nen, T. and Saramäki, J. 2000. Comparison of percentile based prediction methods and the Weibull distribution in describing the di ameter distribution of heterogenous Scots pine stands. For. Ecol. and Manag. 133: 263- 274. McCulloch, C. E. and Searle, S. R. 2001. Gen eralized. Linear and Mixed Models. John Wiley & Sons. 325 p. Mehtätalo. L. 2004. An algorithm for ensuring compatibility between estimated percentiles of diameter distribution and measured stand variables. For. Sci. 50(1): 20-32. 16 Mehtätalo Nepal, S.K. and Somers, G.L. 1992. A Gener alized Approach to Stand Table Projection. For. Sei. 38(1): 120-133. Pienaar, L.V. and Harrison, W.M. 1988. A Stand Table Projection Approach to Yield Prediction in Unthinned Even-Aged Stands. For. Sci. 34(3): 804-808. Pukkala. T. Vetteenranta, J, Kolström, T. and Miina, J. 1994. Productivity of mixed stands of Pinus sylvestris and Picea abies. Scandi navian Journal of Forest Research 9: 143- 153. Reiss, R.-D. 1989. Approximate Distributions of Order Statistics With Applications to Nonparametric Statistics. Springer Series in Statistics. Springer-Verlag, New York. 355 P- Rennols, K., Geary, N. and Rollison, T.J.D. 1985. Characterizing diameter distribution by the use of Weibull distribution. Forestry 58(1): 57-66. Reynolds, M.R. Jr., Burk, T.E., Huang, W-C. 1988. Goodness-of-fit Tests and Model Se lection Procedures for Diameter Distribution Models. For. Sci. 34(2): 373-399. Tang, S., Wang. Y.. Zhang, L. and Meng. C-H. 1997. A Distribution-Independent Approach to Predicting Stand Diameter Distribution. For. Sci. 43(4): 491-500. Van Deusen, P.C. 1986. Fitting Assumed Dis tributions to Horizontal Point Sample Di ameters. For. Sci. 32(1): 146-148. Venables. W.N. and Ripley. B D. 2002. Mod ern Applied Statistics with S. Fourth edition. Spring er-Verlag, New York. 495 p. II 20 Forest Science 50i 11 2004 An Algorithm for Ensuring Compatibility Between Estimated Percentiles of Diameter Distribution and Measured Stand Variables Lauri Mehtätalo ABSTRACT. It is difficult to formulate a diameter distribution model that is compatible with many stand variables. In previous studies, compatibility of diameter distribution has been ensured with the aid of calibration (adjustment) based on making small changes to the predicted frequencies of diameter classes. In these methods, the minimum and maximum diameters cannot be changed, and the measurement error of the stand variables is not taken into account. In this study, two calibration methods based on minimizing deviations from predicted percentiles were developed. Because minimum and maximum diameters were among the predicted percentiles, there were no problems in changing them in the calibration. In the first method, the measurement error of the stand variables was not taken into account. In the second method, in addition to deviations from predicted percentiles, deviations from the measured stand variables were allowed, and the measurement error was taken into account by weighting each term of the objective function inversely by its error variance. The methods were tested in a dataset of Finnish mixed coniferous forests. Both methods were found to be better than the reference method used because the minimum and maximum diameters could be changed. Even if the measurement error was large, the second method was still advanta geous, while the other methods were of no use. For. Sci. 50(1):20-32. Key Words: Diameter distribution, compatibility, calibration, nonlinear optimization, mea surement error, percentile. Diamktkr distribution is an essential stand character istic when assessing, for example, volume, basal area, or stem number. It enables the calculation of the volumes of trees between certain diameter limits, thus providing a tool for assessing, for example, the volume and monetary value of logs obtained in a harvest. The diameter distribution can also be projected into the future, enabling the prediction of future incomes and the effect of different harvest schedules on these. The results can then lie compared with each other in a forest management planning process. The accurate description of the stand structure is very important when simulating stand development for forest management planning, and in long-term simulations the small trees are also of interest. Different theoretical distributions have been used for describing the diameter distribution of a stand. These in clude, for example, the lognormal (Bliss and Reinker 19641, Weibull (Bailey and Dell 1973), beta (Loetsch et al. 1973. p. 48-61) and Jonson's SB (Hafley and Schreuder 1977) distri butions. In recent decades, the Weibull distribution has been the most commonly used because of its flexibility, the fairly straightforward interpretation of its parameters, and the closed Lauri Mehtätalo, Researcher, Finnish Forest Research Institute, Joensuu Research Centre, P.O. Box 68, FIIM-80101 Joensuu, Finland —Phone: +358-10-2113051; Fax: +358-10-2113113; E-mail: lauri.mehtatalo@metla.fi. Acknowledgments: This work was part of the project "Statistical modeling for forest management planning," which was carried out at the Finnish Forest Research Institute and funded by the Academy of Finland (decision number 73392). I would like to thank Juha Lappi, Annika Kangas, Matti Maltamo, and the three anonymous referees for their valuable comments on the earlier drafts of the manuscript. Professor Kangas also provided the computer program for the reference method. Finally, I wish to thank Lisa Lena Opas-Hänninen for revising my English. Manuscript received May 5, 2002, accepted May 16, 2003. Copyright © 2004 by the Society of American Foresters 21 Forest Science 50( I) 2004 form of its cumulative distribution function (Bailey and Dell 1973). Many methods for estimating the parameters of the Weibull distribution have been presented. The maximum likelihood method or the method of moments can be used if a sample dbh-distribution from the stand is available (Van Deusen 1986. Lindsay et at. 1996). If a sample is not avail able. the parameters can be predicted with the parameter prediction method (PPM) or estimated using the parameter recovery method (PRM). In PPM. the parameters are esti mated by developing regression models with stand character istics as independent variables, and parameters of a target stand are predicted using measured stand variables. In PRM, the parameters of the distribution are recovered from some known or predicted stand characteristics. Recovery can be based on the moments of the distribution (Burk and Newberry 1984). known characteristics of certain percentiles (Bailey and Dell 1973, Shiver 1988) or some measured stand vari ables (Eketal. 1975, Hyinkand Moser 1983). When predict ing the future distribution of a stand, the development of the parameters, stand variables, moments or percentiles is mod eled as a function of time. It has been demonstrated that theoretical distributions work well in even-aged unthinned forests: but in more hetero geneous forests, distributions which are purely parametric have been found to be too rigid (e.g., Cao and Burkhart 1984, Maltamo and Kangas 1998) and more adaptable methods are needed. In a percentile-based method (Borders et al. 1987), the diameter distribution of a stand is described with percen tiles of the distribution. The percentiles can be predicted using known stand characteristics, and the development of the percentiles with respect to stand age can also be modeled (Borders and Patterson 1990). The continuous distribution function between the percentiles can be interpolated either with linear interpolation (Borders et al. 1987) or with an adequate spline function (Maltamo et al. 2000). If a measured sample of tree diameters is available, the sample distribution can be used as such in calculating present volume. If the sample is small, the sample distribution can be smoothed, for example with a kernel method, to obtain a distribution that better describes the underlying population distribution (Droessler and Burk 1989. Maltamo and Uuttera 1998). Fitting a theoretical distribution to the sample distri bution can have a detrimental effect: if the population distri bution does not follow the theoretical distribution used, some of the information on the form of the distribution is lost (Borders and Patterson 1990. Nepal and Somers 1992). A good estimate of the future distribution may be obtained by projecting the present sample diameters and stand variables into the future with individual tree growth models (Pienaar and Harrison 1988. Nepal and Somers 1992.Tang etal. 1997) and adjusting the projected diameter distribution to obtain a stand description compatible with the projected stand vari ables. The incompatibility of the distribution may be a prob lem in those cases where the distribution is projected into the future with models (Nepal and Somers 1992, Cao and Baldwin 1999). or if all the measured stand variables have not been used in the prediction of the diameter distribu tion, or if the diameter distribution model used does not guarantee compatibility with all the stand variables used as predictors (Kangas and Maltamo 2000 a). In such cases, calibration (adjustment) of the diameter distribution can be useful. In the calibration, the predicted frequencies of the diameter classes are adjusted slightly to obtain a distribution compatible with all the measured or projected stand variables. This approach has two drawbacks. The first is that new diameter classes are not constructed, and problems arise if the minimum and maximum diameters of the uncalibrated distribution are far from the true diam eters. The other drawback is that the stand characteristics used are assumed to be accurate, i.e.. no measurement or prediction error is assumed. In Finland, forest management planning is based on stand-wise inventories carried out approximately every tenth year, and stand development is predicted with indi vidual tree growth models using representative trees picked from the predicted diameter distribution of the stand ( Kilkki et ai. 1989). The costs of the inventories are kept low to make forest planning inexpensive. Therefore, rather than measuring a sample diameter distribution to predict the stand diameter distribution, the assessed stand character istics are used. The basal-area weighted form of the diam eter distribution (basal area diameter distribution) (see Gove and Patil 1998) is used to ensure the accurate de scription of the large, valuable trees. In the inventories, an angle-gauge is used to determine the basal area and basal area median diameter of a few. subjectively chosen sample plots, and the stand variables are calculated using the sample plot measurements. Due to within-stand variance and measurement error, the total error in the measurement may be as high as 20% of the mean values of the variables (Poso 1983. Laasasenaho and Päivinen 1986. p. 9). Depending on site features, damage, the forest owner's activity, etc.. Finnish forest stands may differ markedly from each other. One stand may have dense natural undergrowth, while in another there may be no undergrowth at all: a stand may be a mixed forest or a one-species forest: it may comprise many age cohorts or be even-aged: it may be naturally regenerated, planted, or sowed, etc. In this situation, it might be advantageous to change the set of assessed variables according to the type of forest stand. However, it would require a large dataset and much work to model diameter distributions with all measurement combinations forall types of stands in order to determine which variables work in any given type of stand. It would also be very hard to formulate models that produce compatible distributions with all mea surement combinations. A simpler approach is to predict the distribution with commonly measured stand variables and to lake other variables into account through calibration. The aim of this study was to develop a calibration method that could be used in Finnish forest management planning. This meant that (1) the method could not prohibit the formu lation of new diameter classes and ( 2) the measurement error of the stand characteristics had to be taken into account. The calibration was based on the percentiles of the diameter distribution instead of the frequencies of the diameter classes. 22 Forest Science 50( 11 2004 At the end of the study, the effect of calibration in the predicted diameter distributions is demonstrated with a case study. The Method A diameter distribution compatible with the measured basal area median diameter and basal area is obtained by scaling the basal area diameter distribution with the mea sured basal area of the calculation unit. i.e.. by multiplying the two. If additional stand characteristics have been mea sured. e.g., stem number or mean diameter, the resulting distribution may not be consistent with all the measurements. This study proposes a method for calibrating the predicted percentiles. The derivations are presented for the basal area diameter distribution, but they can be easily reformulated for a stem number diameter distribution. Calibration of the Percentiles of a Basal-Area Diameter Distribution Let us assume that the diameter distribution is described using percentiles of the basal area diameter distribution and the cumulative distribution is linear between one percentile and the next (Figure 1). Denote the ordered set of M percen tiles by [dj, d-, dm ) and the corresponding percentage values by {p,, p1 pm }. where p { = 0 and pm = 100. (The usual notation, e.g.. is used to show exactly which percentile is under consideration.) Suppose that we have predictions of the percentiles and denote them by {dv d-,,--- d M ]. In addition, we may have measurements of some stand variables. When presenting the calibration methods, it is assumed that the basal area G. basal area median diameter d Gn, and stem number N have been measured. Denote these measurements by G, dGm and N , because they are from a sample and are therefore estimates of the true values. Let us assume that the measurements are unbiased, i.e.. E(G) = G. E{dc,„)=dc„, and £(jv) =W. We want to modify the predicted percentiles to make the predicted distribution consistent with the measured stand variables. Here it is also assumed that the stand variables are known exactly, i.e.. they are assumed to be measured without Figure 1. An example of a cumulative predicted percentile-based distribution and the predicted and calibrated fourth percentile. error. The problem is. in fact, to find a calibrated distribution that is as close to the predicted distribution as possible and concurrently satisfies the calibrating equations. This leads us to formulate the calibration as an optimization problem, where the distance between the predicted and the calibrated distribution is minimized, and the measured stand character istics are included in the optimization problem as constraints. First, for each predicted diameter, we define a deviational variable sr which implies the deviation of the calibrated ith percentile from the predicted percentile. The calibrated zth percentile is d: =d, +sr The optimization problem is minimize subject to 4 G where is the deviational variable corresponding tour„ urn».. "o*. 4.06 0.61 0.48 0.35 0.28 0.00 -0.14 -0.18 -0.20 -0.20 -0.25 -0.21 "lO9 » 2.95 0.83 0.69 0.59 0.00 -0.18 -0.27 -0.28 -0.27 -0.29 -0.21 "20° o 2.12 0.84 0.70 0.00 -0.15 -0.29 -0.30 -0.29 -0.28 -0.19 "30». 1.49 0.81 0.00 -0.13 -0.28 -0.30 -0.31 -0.29 -0.22 "40° o 1.15 0.00 -0.12 -0.21 -0.23 -0.23 -0.21 -0.14 "50° o 0.00 0.00 0.00 0.00 0.00 0.00 0.00 "60°. 0.83 0.66 0.53 0.41 0.42 0.29 "70°. 1.21 0.75 0.68 0.64 0.48 "8(1° o 1.79 0.83 0.77 0.60 cL,()o o 2.60 0.90 0.72 "95® o 2.97 0.77 " 1 ()()».> 4.39 (b) ui» dur..-d- 1° 'A,,-i... ämr.-ä-m. 10°o™ I) 0 » 3.26 -0.07 -0.04 -0.07 -0.18 0.01 -0.04 0.02 0.02 0.10 0.06 d-.ff.-d 1.6') 0.15 0.06 0.15 0.13 0.03 0.05 0.05 0.13 0.03 ?()•„"" "*0% 1.18 0.12 0.23 0.11 0.12 0.08 0.03 0.07 -0.04 "40° »"".Wo 0.87 0.07 0.06 0.20 0.11 0.14 0.05 0.07 " 5()° 40" „ 1.15 0.12 0.17 0.13 0.13 0.02 -0.01 "60 # o~"5<>°. 0.83 -0.04 0.13 0.08 0.14 0.01 "70°."" 60° . 0.90 0.04 0.32 0.00 0.08 "X(l o o~"70°. 1.19 0.09 0.03 0.11 "•>00 o~"X0°o 1.49 -0.03 0.17 1.32 -0.02 2.82 24 Forest Science 50( I) 2004 distribution may have abnormally high peaks in some diam eter classes. Therefore it is better to calibrate the distribution only approximately to the measured values. In this section the optimization problem is modified to take the measurement error into account. In addition to the deviational variables of the percentiles, deviational variables for the measured stand characteristics are introduced. The calibrated basal area and stem number are defined as G=G + sG and N= N + sN , where the v-vari ables are the corresponding deviational variables. The optimization problem is now minimize subject to where and are the measurement error variances of the stand charac teristics. The changes to the previous formulation, i.e.. Equations (1)—(5). are adding three terms to the objective function (8), substituting the calibrated basal area and stem number for the measured ones in the constraint (11). dropping the basal area median diameter constraint (5) from the problem formulation (no longer needed since we no longer assume that the calibrated median is equal to the measured median) and adding the nonnegativity constraints of the calibrated stand variables (12)—(13). Additional Stand Variables in the Calibration In addition to basal area, stem number, and basal area median diameter, we may also have other measured stand variables. This section presents some examples of how additional stand variables are included in the optimization problem. In most cases, the new stand variables can be included in the optimization problem as additional calibration con straints. However, a nonnegativity constraint for the new calibrated stand variable should be added [cf. Equations (12) and (13)]. The calibration constraints are formulated by calculating analytically the corresponding variable from the diameter distribution and equating it with the calibrated stand variable, i.e.. the measured stand variable plus the corresponding deviational variable. Some ana lytically calculated stand variables are presented in the Appendix. The measurement error of the new variables is taken into account by adding the corresponding term into the objective function, as was done in objective function (8) [cf. objective function (7)]. In the following, there are three examples of new constraints. The constraint for the mean diameter of the stem number diameter distribution is formulated by equating the expected value of the distribution (A 10) with the corresponding mea sured mean diameter as in and the constraint for the mean diameter of the basal-area diameter distribution is formulated by equating the expected value of the basal area diameter distribution (A 9) with the corresponding measured mean diameter as in The constraint for the quadratic mean diameter is obtained similarly by equating the calculated quadratic mean diameter (Al2) with its calibrated value as in In Equations < 14)—(1 6>. äi is the slope of the calibrated diameter distribution between the calibrated percentiles d and djM [see Equation (A5)l: _ Y (*'■+' ?i ) . ?50'j" . SC~ s\~ -~L, — + —r +~+ —r (8) '=> °/.i+l °5(W> CG C,V <5, +.?; +£ < di+l +si+r i = 1.2,..., M -1, (9) (/, +.?,+£>(). (10) 4(G+?c) 100 Jt M-l I 1 il xy Pi+l ~ P-- 1 1 ,=i (rf,+i +i,+i) -(4 +■'/■) (4 +■*;) (4+i + vi+i) = N +sN (11) G+ sg >O. (12) iV + .sN >o. (13) °5O - var(dc.m ~ lU',m ) • o G ~ =var(C-C) c v " =var(,\'-A'j >.[, ,1 (14. i r - •> - ■>! =d G (15) v-.r i , y LU-J- . (i 6, )wL U ''-ii a> 100 d i+l -di ' /-U M'• /=1 li (K-K) / KV = jv • ,22 > Standard deviation of the error term (% of the true value) Scenario Gf. Or. i , V . y a . ±_J_ (Ai 2) 1 tr'U III 131 Can. J. For. Res. 34: 131-140 (2004) doi: 10.1139/XO3-207 © 2004 NRC Canada A longitudinal height-diameter model for Norway spruce in Finland Lauri Mehtätalo Abstract: A height-diameter (H-D) model for Norway spruce (Picea abies (L.) Karst.) was estimated from longitudi nal data. The Korf growth curve was used as the H-D curve. Firstly, H-D curves for each stand at each measurement time were fitted, and the trends in the parameters of the H-D curve were modeled. Secondly, the trends were included in the H-D model to estimate the whole model at once. To take the hierarchy of the data into account, a mixed-model approach was used. This makes it possible to calibrate the model for a new stand at a given point in time using sample tree height(s). The heights may be from different points in time and need not be from the point in time being pre dicted. The trends in the parameters of the H-D curve were not estimated as a function of stand age but as a function of the median diameter of basal area weighted diameter distribution (d° m). This approach was chosen because the stand ages may differ substantially among stands with similar current growth patterns. This is true especially with shade tolerant tree species, which can regenerate and survive for several years beneath the dominant canopy layer and start rapid growth later. The growth patterns in stands with a given cf' m , on the other hand, seem not to vary much. This finding indicates that the growth pattern of a stand does not depend on stand age but on mean tree size in the stand. Resume : Une equation de prediction de la hauteur par le diametre (H-D) pour repinette de Norv£ge (Picea abies (L.) Karst.) a et£ calibr£e ä partir de donn£es longitudinales. La courbe H-D est representee par la courbe de croissance de Korf. Des courbes H-D ont d'abord £te ajustees pour chaque peuplement et chaque prise de mesures. et les variations dans les parametres de la courbe H-D ont ete modelisees. L'equation H-D a ensuite ete reajustee apr£s y avoir incor pore les modeles de variation de ses parametres. La procedure d'estimation par modele mixte a ete utilis£e pour tenir compte de la hierarchie des donn£es. Cette procedure permet de calibrer le module pour un nouveau peuplement h un moment donne dans le temps en utilisant un £chantillon des hauteurs d'arbre. Les hauteurs peuvent provenir de mesu res prises ä differents moments dans le temps et n'ont pas besoin d'etre prises au meme moment que la prediction. La variation des parametres de liquation H-D n'a pas £t£ modelis£e en fonction de Tage du peuplement. mais plutöt en fonction du diametre median de la distribution des diametres ponderee par la surface terriere (d?™). En effet. des peu plements d'age tr£s different peuvent avoir des patrons de croissance similaires. C'est le cas particuli&rement chez les essences tolerantes ä F ombre qui peuvent se regenerer et survivre pendant plusieurs annees sous couvert et entamer une croissance rapide par la suite. Par contre. les patrons de croissance dans les peuplements ayant une valeur de ,„ = bltd + b Ur \tll + b 2ir x 21 , + .... and p u, = + b Ur \ lh + b 2 , r x 2h + where .v,„. A,/,. ,v 2 „. .v2;,.... are the additional predictors, and blhr b„h. b UI. />,, are fixed parameters to be estimated. These equations were written into eq. 7, and the whole model was estimated using restricted maximum likelihood. Only significant pre dictors at the 1% level of significance were included in the final models; the predictors were found by backward elimination. The additional predictors tested were stand co ordinates. altitude above sea level, cumulative temperature sum. soil type, dummy variable for thinned stands, basal area. age. basal area of spruces, and cf' m of spruces. In Table 2 the estimated models are presented from the crudest model (I) to the broadest model (V). Model I is the basic model with only tf' m as a predictor. In model 11. site fertility class and geographical variables are included. These variables take into account the large geographic variation in the data: the maximum distance between stands in the mod eling data was 845 km. the altitude above sea level varied from 10 to 360 m, and the cumulative temperature sum var ied from 725 to 1342. The geographical variables are always known in practice (they can be derived from the stand loca tion. for example, using generally available maps), and thus they are also used in all other models. A commonly assessed stand characteristic is also the site fertility class, which is in the model as a dummy variable for mesic and poorer sites. These variables do not cause any harm if the H-D curve of a stand is predicted for several points in time, since their val ues do not change over time and thus need not be measured at each point in time. In addition to the predictors of model 11, models lII—IV include different combinations of stand variables. In models 111 and IV. stand characteristics are assumed to be measured only as mean values of the dominant story, not by tree spe cies. The difference between models 111 and IV is that stand age is used in model IV but not in model 111. In model V the stand characteristics (G and cf" n ) are also assumed to be measured by tree species, but as seen in the models, these variables provide only a small amount of additional informa tion: basal area of spruces was not significant at all. and ,„ and i.e., the trend functions of A and B are assumed to be constant, and all variation of the variables A and B is described by the random part of the model. In the other models, this variation is shared between fixed predictors and random parameters. The proportion ex plained by fixed parameters of a given model can be calculated as R : = 1 - (variance of random parameter in the given model/ variance of random parameter in model 0). For example, for stand-level random parameter ak in model I, it is calculated as R2 , =I - var(a ( i)/var(a ( () ). This statistic has been calcu lated for both stand-level and measurement-occasion-level random parameters in Table 2. One should note, however, that the value of R 2 does not behave consistently in all cases, because the correlation between random parameters is high, and the estimated share of the variation between stand and measurement occasion varies slightly among models. In model I. the estimated trend functions explain 85% of the stand-level variation and 95% of the measurement occasion-level variation of the maximum height in the stand (parameter A). With parameter B the proportion is clearly [B] var(e t .„) = cr{ [max(Dt,,,7.s)r 6 ) 2 136 Can. J. For. Res. Vol. 34, 2004 © 2004 NRC Canada Table 2. Estimated parameters of eq. 7 with different predictor combinations. Note: Only significant parameters at the 1% level are included in the models. The predictors are as follows: yk, the north coordinate in the Finnish G, total basal area of the stand (100 nr/ha): (f" m. median of basal area weighted diameter distribution of the stand (100 cm): (f" m of spuces (100 cm): which takes a value of 1 with stands on mesic or poorer sites, a and 5 are coefficients of the variance function (eq. 8). lower, which is due to the large within-stand variation of pa rameter B when compared with that of parameter A (Fig. 2). One can see two clear improvements in the calculated R~ values among the models. The first one occurs between models I and 11, i.e., when geographical variables and site fertility class are included in the model. In this case, only the stand-level R 1 increases markedly because the additional fixed predictors in model II are constant in a certain stand, and thus they do not explain the variation at the measurement-occasion level. The second improvement be comes evident when the most commonly measured stand characteristics are included, i.e., between models II and 111. In this case, the improvement occurs also at measurement occasion-level R 2 because values of the additional predictors vary among measurement occasions. Stand age also im proves the fit slightly, and it should be used when it has been measured. On the other hand, measurement of stand characteristics by tree species is clearly not useful. Application of the model When the model is applied in predicting the H-D curve of a stand, both the fixed and random parts can be utilized. The prediction of the fixed part is obtained simply by placing the measured values of the predictors in the appropriate model of models I-V. If no tree heights have been measured, the expected value 0 is used for all random parameters, and the prediction of the fixed part is the predicted H-D curve. If tree heights have been measured, they can be used to predict the random parameters of the stand and measurement occasion levels. The random parameters of a linear mixed model are predicted using the best linear unbiased predictor (see Searle 1971; McCulloch and Searle 2001). which has been applied, for example, by Lappi (1991. 1997) in the context of linear mixed-effects models in forestry. Assume that we have n measured trees, possibly from different points in time. The measured log-transformed heights are in vector ynx | and their expectation in n„ x] . The measured heights are described by the model where Z„x,„ is the model matrix of the random part of the model. bmxl is the vector of the random parameters, and e„xl is the vector of prediction errors. In addition, we define the variance-covariance matrices of the random parameters and residuals as var(ft) = D mx „, and var(e) = R„ x„, respectively. The number of random parameters to be predicted (m) de pends on how many measurement occasions we have in the stand. For example, if we are predicting the random parame ters of eq. 7in stand kat measurement times t t and t 2. in = 6 and h' = (a k, P( .. a(li , ). A numerical example of predicting the random parameters is presented in Lappi (1991). When random parameters are predicted, the estimates of matrices D and R are needed. The estimate of D (denoted by D) can be built up from the estimated standard deviations and correlations of the random parameters (Table 2). The es timated error variances of the measured trees can be calcu lated using the variance model in eq. 8, with the estimated parameters 5 and CT from Table 2. These variances are placed on the diagonal of matrix R. Because prediction errors of different trees are assumed to be uncorrected, the non diagonal elements of R are zeros. The best linear unbiased predictor of the random parameters is calculated as When the random parameters are calculated, the estimated variance-covariance matrices R and D are substituted for matrices R and I), and the prediction of the fixed part is sub stituted for (he vector (a. For further details see. for example. Lappi (1997) and Searle (1971. pp. 458-462). To obtain calibrated H-D curves, the predictions of the random parameters from the vector b are placed in eq. 7. Note that the predicted stand-level random parameters can also be used to predict the H-D curve at a development stage for which we do not have measurements. In this case, we just place the predicted stand-level random parameters in eq. 7 and use the expectation 0 for the measurement occasion-level parameters. Calibrating the H-D curve with measured heights [9] y = |jl + Zb + e [lo] Population level (fixed part) Stand level Model Pin Phi Pit Plb SD (at) SD (pt ) 0 2.91 0.692 0.278 0.136 1 1.40 2.15 0.384 0.0157 0.108 0.0958 II 9.45 - 0.94 lyk - 0.1 19alt - 0.121 dd + 0.0324soil 2.12 0.656 - 0.028 ldd 0.0171 0.0891 0.0867 III 8.39 - 0.8 lyk - O.lOlalt - 0.0999dd + 0.46G + 0.0213thin + 0.0367soil 1.82 0.646 - 0.0262dd - 0.305G 0.0196 0.0869 0.0842 IV 11.79 - 1.24yk - 0.127alt - 0.138dd + 0.369C - 0.106cPm + 0.0210thin + 0.2287 1.92 2.33 - 0.225yk - 0.0379dd - 0.338G + 0.00 128f 0.0167 0.0813 0.0800 V 11.81 - 1.24yk - 0.1 28alt - 0.139dd + 0.36C - l.OSrf 0"1 + 0.2247" + 0.332d G"u + 0.0201 thin 1.93 2.33 - 0.226yk - 0.0381dd - 0.339G +0.1 28r 0.0167 0.0812 0.0800 Mehtätalo 137 © 2004 NRC Canada uniform coordinate system, 1000 km: alt, altitude above sea level (100 m); dd, cumulative temperature sum (average of the years 1951-1980) (100 x dd): T, stand age at breast height (100 years); thin, dummy variable indicating if the stand has been thinned within the last 10 years: soil, dummy variable Utilizing the longitudinal character of the model When the longitudinal character of the model is applied in practice, the fixed predictors of the model need to be known at each point in time under consideration, i.e.. both at the point! s ) used in predicting random effects and at the point at which the H-D curve is predicted. The longitudinal charac ter of the model can be utilized in two different ways: (/) us ing height measurements from several points in time to predict the H-D curve at a certain point in time and (/'/') making predictions into the future. In case /. the use of the model is straightforward: just place the measured stand variables into the appropriate model and, if heights have been measured, predict the random parameters using eq. 10. In case ii. the situation is a bit more complicated, because usually one wants to predict the H-D curve for some partic ular point in time, not for the time when (/"" reaches a cer tain value. Therefore, we need to predict c/"" first and use it in the prediction of the H-D curve. Figure 3 shows that the variation in cf-' m in stands at any particular age is considerable. Among-stand variation is, however, greater than within-stand variation. In addition, in any one stand, the change in cf' m is stable with respect to stand age. This is a typical situation, in which a longitudinal analysis can be applied (Diggle et al. 1994). and it was used to model the development of cf' m as a function of stand age. In any one stand the growth rate seems to be almost linear except for a slightly concave trend. Hence, a simple power function seems to fit the data well. The power function was linearized with a logarithmic transformation to obtain the model where u and v are fixed population-level parameters. Uku k is a stand-level random parameter, and Tkl is the breast-height age of stand k at time l. The parameters it k and ekl are as sumed to be independent and normally distributed with ex pectations of 0 and variances a;, and a;,, respectively. The transformation linearized the data rather well, with the ex ception of a few very variable young stands. The likelihood ratio test showed no significant difference between the model in eq. 1 I and a model without intercept n: thus, the Fig. 3. Stand weighted median diameter (dGm) versus stand age in the modeling data. The lines connect measurements of any one stand. intercept was dropped from the model. The parameter esti mates were r = 0.727. a;, = 0.372 2 . and a;, = 0.0751 2 . The estimated standard deviations show that the major propor tion of the random part is due to the among-stand variance component as opposed to the within-stand component. In practice, stand age and cf' m are always known at least for one point in time, and they can be used to predict the lo calized stand-level ™ s <* £ t/l CG c/l O cG O CG "§ f 73 a> o -3 » 3 g S =3 I I CG c/3 «, O s >» .-S pii O •' -o I O CG rXS S C Ö rH j 1 ö ** B 9 X •' ° S ® C 3 rn ■ä* © ti .S 3 c/3 CG » r £ B C/l o tS -öb o^U c S > <2 « S « 'g « £ « Z5 c/i -￿—• c/l CG CG lS 2 | §ii i c - £"1 §2^ s 5 S r- s o •§ >> CG ..' . "o S « öfl H 2 2 o tfl 03 a> o CG 2«2 o S CG ■—￿ i-< S S "2 "S S o 3 ~ .s .2 •<—» c/3 CG .2 o S s J | S S o g (D fN Ö OO «o fN ö fN On fN Ö r-» ON «O fN Ö •O rn fN fN OO t-4 O O t-» o r-- © VO cp r-> cp r» ö i S Tl- on o r-- Tf fN 00 in */~> VO T3 C/j o ö o ö O ö o ö © ö T3 C/j OO r» OO (N O ö VO VO OO fN O ö r-- o fN O O Ö f- VO U-) Ö s ON •/"> VO rn O OO OO OO "3- t-» ~d C/j O Ö o ö o ö o ö p ö "3- VO U~) m U~) On fN r-- o fN VO O fN C/j Ö ö Ö Ö Ö 3 Q, u-> rn ON o o o ö o On O O o ö r-~- rn •/"> o o o o •O «n IT) O o o ö oo m o o o ö £ r-» o r-. fN o o l/~) t-» (N o cp On fN o o ö fN O cp' OO OO o ö P Ib VO *o VO Ö -0.04113 +0.0917*v* Tf- fN «/-> O * rr) (N Tf- cp 0.3147 +0.143*; * fN O VO ö ■3fr VO O o ro O ö + * fN VO O O + r*» On fN m ö *" o ro ö + s-/ * ON Tf VO cp * IT) ro OO fN O ö + * Tt- lo o ö + O * «n VO •T) ö + •^1" — r~- (N -t VO VO rf c fN in >o e*-) S fN •—< O * fS *-f o fN Ö + * VO vO O o *' -3- r» rn o o 1 OO CO vq *" r-» o o rn O O f^l >2 * r-- fN On fN O O 5 5 * r-~ vO •O O * On OO ö + »5 *■ OO VO i/~> rn o ö + 3.025 +0.2431*/ ro «/-> On cp >3 * ro VO •<3- ro O O 5 0 Q * ro o > Height-diameter models for Scots pine and birch 7 Table 3. The estimated models for birch. DGMb is DGM of birch, 100 cm. For other notations, see Table 2. (N O C/3 c/n' C O 08 O S .£ o o p £ CJ O o JS CJ £ <4_ O •< o Ci C/3 O Z '5 -0.4535 S6 -0.4546 ?8 -0.4536 '5 -0.4536 rt ON rt rt ö V) to ON 00 NO ro 00 O Vi ro ON O "d O Ö - 0.1265 +0.01159* < O Q C1 ON rt rt •o rt V) *T + > "" 8 Mehtätalo Figure 2. Stand DGM against stand age in a data including all stands of the original data. Subfigure a shows the relation on the arithmetic scale and sub figure b shows the relation between logarithmic DGM and age. The solid line in subfigure b is the expected development obtained with the fixed part of model (12). where A=4.5 for Scots pine and A-9.5 for birch. The estimated nonlinear part of the trend function of A was for Scots pine and for birch. In step 5. the fixed parameters pla and plb were written as linear combinations of the predictors used in the model, i.e. Pxa =K a +b l a X la +K X 2a +- and P\b ~ b X\b ~^"^2b X2b "'"•••? Where Xj a, Xjfo x2a. xy„... are the additional predictors and b()a . b 0b • bla , b,are their coefficients. The estimates of the fixed parameters and variances of the random parameters are in Tables 2 and 3. For both tree species, 5 models with different predictor combinations were estimated. Model I uses only DGM as predic tor. In model 11. variables describing the geo graphical location are included. Models 111 and IV include, in addition to the predictors of model 11, stand characteristics measured from the whole growing stock and furthermore, model V includes those measured by tree spe cies. The difference between models 111 and IV is that model IV includes stand age. Only pre dictors with statistically significant coefficients were included. The significance level used was 1% in Scots pine models and 5% in the models for birch. The significant predictors were searched stepwise, refitting the model and eliminating the least significant parameter until all remaining coefficients were significant. One can see that the variances of the random parameters decrease when the number of pre dictors increases. This happens because the variation explained by the additional fixed predictors belongs to the random part in the more sparse models. Prediction of DGM in the future In the prediction of the H-D curve of a stand at a given point in time, the DGM of the stand at that point in time needs to be known. With current and past points in time this is not a problem, since in Finnish inventories DGM is a very commonly measured mean stand charac teristic. However, in order to be able to predict var(%) = CT2 (max(Dto ,A)) , (8) _ fc= (l (9) Z h (10) Height-diameter models for Scots pine and birch 9 the H-D curve in the future, a model for the DGM was estimated. Since the DGM used in the models is the DGM of all trees species belonging to the dominant canopy layer, one and the same model can be used with Scots pine, birch and Norway spruce. To construct such a model the modeling data comprised both the full data used in this study and that used in Mehtätalo (2004). As seen in Fig. 2a, the DGM of a stand with a given age varies considerably, especially with older stands. However, the change in DGM is stable with respect to stand age. Thus, if we have an observation of stand DGM at some age, we can predict the DGM at some other age rather well. The nonlinear trend in Fig. 2a was linearized by taking logarithms from DGM and age to obtain the model where Tkt is age of stand kat time t, u and v are fixed parameters, uk and vk are stand level random parameters and ek, is the residual. The random parameters and residual error were assumed to be normally distributed with a constant variance. The parameter estimates obtained were w=0.1113(5.e.=0.05411) and v=0.6999 (5.e.=0.01342), var(«*)= 1.060 2 , var( vk )= 0.2681 2 . cov(u k,vk)= 0.2693 and var(ek,)= 0.05280 2 . The linearized model seems to fit the data well (Fig. 2b). When utilizing the model to predict the DGM of a stand with a given age. the stand effects uk and vk are pre dicted using the best linear unbiased predictor (e.g. Searle 1971: 458-462, see Lappi 1993, 1997 and Mehtätalo 2004). Application example To demonstrate the use of the estimated models, H-D curves were predicted with each of the models for one stand selected from the modeling data. The stand was a 35-year-old mixed-species forest with a DGM of 12 cm and a basal area of 22m 2 /ha, of which 17.7nr/ha was Scots pine and 3.4m"/ha birch. Fig. 3 shows the fixed part predictions ob tained with each of the models in Tables 1 and 2. The predictions obtained with the various models differ slightly from each other. How ever, for both tree species all models give too low heights in this stand. Fig. 4a demonstrates the effect of localiza tion on the predictions. Models II and V for Scots pine were localized for the sample stand by predicting the stand and time effects of the H-D models using one measured height sample tree from the stand. The prediction was carried out with the best linear unbiased predictor; the equations have been presented many times before (Lappi 1993, 1997, Mehtätalo 2004) and are thus not presented here. The localized mod els give much better height predictions than the fixed parts only. Note that in Fig. 4a the local ized models are very close to each other even if the fixed part predictions are not. This is be cause the information of one measured sample overrides the information of the additional fixed predictors of model V. This demonstrates the somewhat self-evident fact that it is more efficient to improve the prediction of the H-D curve of a stand by measuring heights and diameters than by measuring the covariates of model V. Fig. 4b demonstrates the projection of the H- D curve of a stand into the future. The meas urements were made at the age of 35 years and the H-D curve is projected to the age of 45 years. This required knowledge about the DGM at the age of 45 years, which was predicted using model (12). The random effects of model (12) were predicted with BLUP using the known DGM at the age of 35 years to obtain a stand-specific DGM-age -curve. It was used to predict the DGM of the stand at the age of 45 years. The projected H-D curves were calcu lated using the predicted DGM. The projections obtained with the localized model are again clearly better than the predictions of the fixed part only but they seem to be underestimates of the height. In fact, all localized models in Fig. 4 seem to give slightly underestimated heights. This is because the expected H-D curve is so far from the observations that one measurement does not move it far enough. Using more than one sample tree would reduce the bias of the predictions. \n(DGM h )=u+vln(7^) +tik + vt ln(rfe ) +e to , (11) 10 Mehtätalo Figure 3 Predicted H-D curves of Scots pine (a) and birch (b) in a 35-year-old mixed pine-birch stand selected from the modeling data. The predictions are calculated with each of the models I-V using only the fixed part of the model. The marks show the observed heights and diameters. Discussion In this study H-D models for Scots pine and birch were estimated for practical use in Finland. The models can be used to predict the H-D relationship of a stand with known DGM. If other stand characteristics than DGM are measured, they can be used through selecting from models I V the model that best suits the situation. Measured height-diameter pairs Figure 4. Predicted H-D curves of the Scots pines in the stand of Fig. 3. Dashed lines are the predictions of the models obtained using the fixed part only and solid lines are localized using one measured height diameter pair at the age of 35 years (the big solid ball in the figures). Subfigure a shows the predic tions obtained using models II and V at the age of 35 years, hi subfigure b. the H-D curves after 10 years growth are predicted with model 11. Small symbols are the observations (• at the age of 35 years and ￿ at the age of 45 years). The DGM of the stand at the age of 35 years was 9.6 cm and the predicted DGM at the age of 45 years was 11.5 cm. can be used in localizing the model for a target stand and. due to the longitudinal character of the model, information from any point in time, i.e. any stage of development, can be utilized. Height-diameter models for Scots pine and birch 11 Furthermore, the H-D curve of a stand in the future can be predicted, but this requires the prediction of DGM at that point in time. These properties of the model are discussed in Lappi (1997) and Mehtätalo (2004). Mehtätalo (2004) observed that the devel opment of the H-D curve of a shade-tolerant tree species (Norway spruce) depends on mean tree size in the stand rather than on stand age. The present study continued this work and showed that the observation made with shade tolerant tree species is true also with shade intolerant tree species. The reason for this is that the site properties affect the development rate of a forest stand, so that stands on poor sites develop more slowly and for longer than stands on rich sites. Because of this, the models predicting the development of the H-D curve of a stand perform better when mean diameter of the stand is used as the variable describing the maturity of the stand instead of stand age. This effect on the performance is probably also true of models predicting other things than H-D curves. Hence, when modeling the develop ment of any stand characteristics, for example, diameter distribution and stand growth, the use of stand age as the only variable determining the stage of development of the stand should be viewed critically. However, the method for taking into account the effect of DGM on the development of H-D curves which has been presented here is not the only correct one; other approaches may also lead to equally good results. For example, the development of parameters A and B can be linked with stand age and the DGM can be used in the prediction of the random effects as Lappi (1997) did in his models 6 and 7. Thus, if both DGM and age are known, it is not obvi ous that the development of H-D curves should be linked with DGM and linking it with age may work as well, if the DGM is taken into account in the models. However, using age only will not lead to as good models as will the use of DGM only and thus, if one of them must be chosen, it would be preferable to use the DGM. Acknowledgements This work was part of the project "Statistical modeling in forest management planning", which was funded by the Academy of Finland (decision number 73392) and carried out at the Finnish Forest Research Institute. The author would like to thank Dr Juha Lappi for com ments and criticism on the manuscript and Dr Lisa Lena Opas-Hänninen for revising my English. Literature cited Arabatzis, A. A. & Burkhart, H. E. 1992. An evaluation of sampling methods and model forms for estimating height-diameter rela tionships in loblolly pine plantations. For. Sci. 38: 192-198. Chambers, J. M. 1998. Programming with data. Springer, New York. Curtis, R. O. 1967. Height-diameter and height-diameter-age equations for second growth douglas-fir. For. Sci. 13: 365-375. Diggle. P.J., Heagerty, P., Liang, K.-Y. and Zeger, S.L. 2002. Analysis of longitudinal data. Second edition. Oxford Statistical Sci ence Series. Oxford university press, Oxford, U.K. 379 p. Eerikäinen. K. 2003. Predicting the height diameter pattern of planted Pinus kesiya stands in Zambia and Zimbabwe. For. Ecol. Manage. 175: 355-366. Fang, Z. & Bailey. R. L. 1998. Height-diameter models for tropical forests on Hainan island in southern China. For. Ecol. and Manage. 110: 315-327. Flewelling. J. & de Jong, R. 1994. Considera tions in simultaneous curve fitting for re peated height-diameter measurements. Can. J. For. Res. 24: 1408-1414. Huang, S., Titus, S. J. and Wiens. D. P. 1992. Comparison of nonlinear height-diameter functions for major Alberta tree species. Can. J. For. Res. 22: 1297-1304. Gustafsen, H. G„ Roiko-Jokela, P. & Varmola, M. 1988. Kivennäismaiden talousmetsien 12 Mehtätalo pysyvät (INKA ja TINKA) kokeet. Suunni telmat, mittausmenetelmät ja aineistojen ra kenteet. Finnish Forest Research Institute, Research papers 292. FFRI, Helsinki, Fin land. 212 p. Lappi, J. 1991. Calibration of height and vol ume equations with random parameters. For. Sci. 37: 781-801. Lappi, J. 1997. A longitudinal analysis of height/diameter curves. For. Sci. 43: 555- 570. Lynch. T. B. & Murphy. P. A. 1995. A com patible height prediction and projection sys tem for individual trees in natural, even-aged shortleaf pine stands. For. Sci. 41: 194-209. Mehtätalo, L. 2004. A longitudinal height diameter model for Norway spruce in Finland. Can. J. For. Res. 34(1): 131-140. Pinheiro, J. C & Bates, D. M.. 2000. mixed effects models in S and S-PLUS. Springer- Verlag, New York, USA. 528 p. Richards 1959. A flexible growth function for empirical use. Journal of Experimental Bot any 10: 290-300. Searle, S. R. 1971. Linear models. John Wiley & Sons. New York. USA. 532 p. Venables, W.N. & Ripley, B D. 2002. Modern applied statistics with S. Fourth edition. Springer-Verlag, New York. USA. 495 p. Zakrzewski, W. T. & Bella, I. E. 1988. Two new height models for volume estimation of lodgepole pine stands. Can. J. For. Res. 18: 195-201. V AN APPROACH TO OPTIMIZING FIELD DATA COLLECTION IN AN INVENTORY BY COMPARTMENTS Lauri Mehtätalo 1 and Annika Kangas 2 1 (Corresponding author), Researcher, Finnish Forest Research Institute, Joensuu Research Cen tre, P.O. Box 68, FIN-80101 Joensuu, Finland 2 University of Helsinki, Department of Forest Resource Management, P.O. Box 27, FIN-00014 University of Flelsinki, Finland Abstract. This study presents models for the expected error of the total volume and saw timber volume due to sampling errors of stand measurements. The measurements considered are horizon tal point sample (HPS) plots, stem numbers from circular plots, sample tree heights, sample order statistics (i.e. quantile trees) and sample tree heights from the previous inventory. Different meas urement strategies were constructed by systematically varying the numbers of these measurements. A model system developed for this study was utilized in a dataset of 170 stands to predict the total volume and saw timber volume of each stand with each measurement strategy. The errors of these volumes were modeled using stand characteristics and the numbers of measurements as predictors. The most important factors affecting the error in the total volume were the numbers of HPS plots and height sample trees. In addition, the number of quantile trees had a strong effect on the error of saw timber volume. The errors were slightly reduced when an old height measurement was used. There were significant interactions between stand characteristics and measurement strategies. Thus the optimal measurement strategy varies between stands. It was demonstrated how constrained optimization can be used to find the optimal strategy for any one stand. Key words: forest planning, response surface, optimization, quantile, diameter distribution, H-D curve, prediction 2 Mehtätalo and Kangas Introduction Forest management planning in Finland is based on inventory by compartments of the target area. In the inventory, basal area, basal area median diameter (i.e., the median diame ter of the basal area weighted diameter distri bution. DGM), height of the DGM-tree, stand age and site fertility class of each stand are estimated with the aid of a few subjectively located horizontal point sample (HPS) plots. Using these data, estimates of the diameter distribution and height-diameter curve (H-D -curve) of the stand are produced. These esti mates are used to build up a set of representa tive trees (i.e. a sample from a diameter distri bution), which are then used in the simulation of alternative management schedules for the stand. The system described above produces, in most cases, a reasonable description of a stand. There are. however, many things that are not taken into account or that could be done better in the system. One weakness is that, due to the subjective location of sample plots and visual assessment of stand characteristics, the sam pling errors cannot be calculated using the standard formulas from sampling theory. Thus, the accuracy of stand measurements is not known and cannot be taken into account in the predictions. Furthermore, no assessment of the accuracy of the predictions can be given. Be cause the system includes chains of models with nonlinear transformations, random errors of the predictors may cause bias in the predic tions (Gertner 1991. Kangas 1997. Kangas 1999). Another weakness is that the combina tion of measurements is fixed. This combina tion may not be optimal and thus reallocation of the measurement resources (i.e.. the meas urement time) to different measurements could improve the accuracy of the stand description and growth estimates. For example. Kangas and Maltamo (2000 c) suggested that substitut ing height sample trees for HPS plots might improve the accuracy of volume predictions. Kangas and Maltamo (2002) studied the use fulness of different measurements in the esti mation of total volume. They utilized a calibra tion (or adjustment) algorithm (Kangas and Maltamo 2000 a) to predict the stand descrip tion using several measurement strategies and estimated models for the prediction variance of stand volume as a function of basal area, DGM and stem number of the stand. The best meas urement strategy varied between stands and depended on stand characteristics. However, their model did not take into account the sam pling and measurement errors of the stand characteristics. According to previous studies, the standard errors (including sampling and measurements errors) of partly subjective and visual stand assessments may even be more than 30% of the true values of the stand char acteristics (Mähönen 1984. Pigg 1994. Kangas et ai. 2002, Kangas et ai. 2004. Haara and Korhonen 2004). In this study, in order to be able to control the errors of stand measure ments, sample plot measurements are assumed to be accurate and the sample plots are located randomly in the stand. Thus, the stand meas urements include only sampling errors, which are controlled by varying the number of meas urements carried out in a stand. Response surfaces (Khuri and Cornell 1987) can be used to study the effect of one or more experimental factors on the outcome of a com plex system. For example, they are used in studying the effect of errors in predictors on the response of a process-based forest growth model (Gertner et al. 1996). The idea is to vary systematically the values of the experimental factors and calculate the outcome of the system using each combination of them. The response surface is a regression model that has the ex perimental factors as predictors and the out come as the dependent variable. The values of the experimental factors are selected so that the model matrix is orthogonal. The approach of this study resembles the response surface ap proach. However, the analysis is not carried out in a single stand but in a dataset consisting of several stands and the stand variables are also treated as experimental factors. Thus, the model matrix of the surface is not orthogonal. However, the estimated surface can be used in the prediction of responses, in a similar manner as the ordinary regression model. This study presents a model that can be used in the optimization of stand measurements in Optimizing field data collection 3 an inventory which is carried out by compart ments and where sample plots are located objectively. The model predicts the expected error in the predicted total volume and saw timber volume of the stand using the stand characteristics and the accuracy of different stand measurements as predictors. The model can be used in comparing different measure ment strategies and in finding the optimal one for a single stand. The limitations of the tradi tional calculation system made it necessary to build up a new one for this study. The aims of this study were: (a) to build up a new system for producing a stand description for forest planning, i.e. a system that is able to utilize different kinds of information measured with different levels of accuracy, (b) to estimate a model for the prediction variance of total vol ume and saw timber volume using stand char acteristics and the accuracy of the stand meas urements as predictors, and (c) to show exam ples of how this model can be used in antici pating the prediction variance of a stand description, i.e. in allocating the measurement time to different measurements optimally. The system for producing the stand description A model system for the prediction of stand description was developed. The input of the system includes the values of the stand vari ables (basal area. DGM, age. site fertility class, information about thinning during the past ten years, stand coordinates, temperature sum, and altitude), the variance-covariance matrix of their errors (including sampling and measure ment errors), any number of height sample trees from any points in time (i.e., trees with known diameter and height) and any number of sample order statistics (i.e., trees with known diameter and rank at the sample plot). The sample order statistics are later called quantile trees. The output of the system is a stand description that includes the basal area, basal area diameter distribution and H-D curve of the stand; this output can be used in the calculation of volume, saw timber volume, dominant height, etc. The R-implementation of the S language was used in implementing the system on a computer (R Development Core Team 2003). However, some parts of the system (numerical integration, constrained nonlinear optimiza tion) could not be done with standard R functions and IMSL-subroutines were used (IMSL 1997). They were linked to R. The production of the stand description comprises the following four stages. Since all parts of the system have been published else where, the algorithms are not presented here in detail. It is not necessary to understand the algorithms in order to grasp the concept as a whole and therefore skipping the following four subsections will not prevent the reader from understanding the rest of the article. Stage 1: Predicting the expected diameter percentiles The expected diameter percentiles are pre dicted with the models of Kangas and Maltamo (2000b). The models use a percentile-based approach (Borders et al. 1987), predicting the logarithmic oth,0 th , 10 th , ..., 80 th . 90 lh , 95 th and 100 th percentiles of the basal area diameter distribution. The fixed predictors of the model are DGM, age. basal area and dummy variable for site fertility. The model for the diameter percentiles of stand m is of the form where B is the matrix of fixed parameters, xm is the vector of predictors in stand m. and e,„ is the vector of residuals in stand m with an ex pectation of 0 and variance-covariance matrix D. Assuming that the model is correct and the parameters are known, the fixed part of the model gives the conditional expectations of the logarithmic percentiles given x„„ and the re siduals e,„ are the deviations of the true loga rithmic percentiles of stand m from their condi tional expectations, later called stand effects. The conditional expectations of the percentiles of stand m are predicted using model (1). The continuous distribution function is obtained by linear interpolation on these predictions. In(d,„ ) = Bx,„ + e,„ = E (lnd|x,„ ) + e,„ . (1) 4 Mehtätalo and Kangas The measurement errors in the estimation data were so small that their effect on the pa rameter estimates is neglected. However, when utilizing the models in real forest inventories, the measurement and sampling errors are much larger. Because of the nonlinear transforma tions, these errors cause bias in the predicted diameter percentiles and affect the variance covariance matrix of the error term, which is needed later. The effect of errors in predictors on the predicted distribution is shown in the Appendix, as are the derivations of the cor rected predictions of the percentiles and the error variance-covariance matrix. The formulas given in the Appendix (Equations A 9 and A 10) were used for prediction. Stage 2: Localizing the diameter distribution for a stand If one has measured the diameter of a sample tree and, in addition, knows its rank r in a sample of size «, the measured diameter is the r Ol order statistic of the sample, D r:„. It can be shown (Mehtätalo 2004 d) that an order statistic of an i.i.d. diameter sample is an unbiased estimate of a certain, say 100/ A diameter per centile of the stand, where the value of p de pends on the diameter distribution of the stand (we will return to this later). In other words, the diameter of a quantile tree. D,.„, is a meas ured diameter percentile of the stand. The measured percentiles can be utilized in localiz ing the expected percentiles of the stand by predicting the stand effects e,„ of model (1). The localization utilizes the variance covariance matrix of the measurement error (i.e. sampling error) of the measured percen tiles and the variance-covariance matrix of stand effects. The former is derived using the theory of order statistics (Reiss 1989: 21, 30- 31, see Mehtätalo 2004 d) and the latter was obtained in stage 1. The predicted percentiles follow model (1). Furthermore, the measured percentiles follow the model where the logarithmic measurements are in vector ln(d* m ), their conditional expectations in vector £(lnd*|xm ). the stand effects in vector e* m and the measurement errors in vector z m. The measured diameters estimate the 100xp* lh percentiles of the stand. The elements of p* are calculated as p* = F x [E(D r: „)\, where F is the cumulative diameter distribution function based on the conditional expectations of the percentiles and E(D r:„) is the expectation of the measured order statistic, which depends on the diameter distribution of the stand [for derivation of E(Dr:n), see Mehtätalo (2004 d)]. At this step, F is used as the diameter distribution of the stand (we will return to this later). Vector £(lnd*|x m ) is obtained by interpolation of the predicted logarithmic percentiles for the values of p*. Vectors e,„, e*,„ and z m will be predicted under the assumptions E(em)=E(e* m)=E(E in)=o, var(e m)=D. var(e* m)=D*. var(£,„)=R, cov (e,„.e* m)=C and cov( e*m, £,„)=o, where matrices D. D*, C and R are known: matrix D is the corrected variance-covariance-matrix of the stand effects of model (1) (Equation A 10) and matrices D* and C are derived from it with linear interpolation. The matrix R includes variances and covariances of measured order statistics, which can be derived using their probability distributions (Reiss 1989: 21. 30- 31, see Mehtätalo 2004 d). The aim is to predict the unobserved random vector e,„ using the observed random vector ln(d*„,)- £(ln d*|x m ). The prediction is based on the theory of linear prediction (McCulloch and Searle 2001: 169. Searle et al. 1992: 269-275). The best linear unbiased pre dictor of vector e m is: and the variance-covariance matrix of the pre diction error is The localized logarithmic percentiles are ob tained by adding the predicted stand effects to In ( d *„) = E (In d* | x,„ ) + e +e,„ , (2) e„, = C(D*+R)"'[lnd* ll -£(lrid*|x m )] (3) var(e,„-e,„) = D-C(D*+R)- I C. (4) Optimizing field data collection 5 the conditional expectations of the percentiles as in Because we have obtained a improved predic tion of the diameter distribution of the stand, new estimates of E(D r; „) and p* can now be calculated using the localized diameter distri bution (Equation 5). Therefore, the localization is repeated until the iteration has no consider able effect on the predicted distribution. In some cases, the localized percentiles are not obtained because the distribution of some iteration step is not monotone or the iteration does not converge (see Mehtätalo 2004 d). In these cases, the localization is carried out using a reduced set of quantile trees, including the same trees as before but dropping one tree randomly. Stage 3: Ensuring compatibility with meas ured stand variables The predicted diameter distribution is com patible with the measured basal area and DGM. If additional stand characteristics dependent on the diameter distribution are measured, the compatibility is no longer guaranteed. Me htätalo (2004 a) presented an algorithm for ensuring compatibility of predicted percentiles, basal area. DGM and stem number. In this study, the stem number of trees with a diame ter of 5 cm or more was known. Thus, a modi fication of the algorithm of Mehtätalo (2004 a) was needed in order to take into account the fact that the stem number was measured only above a fixed threshold diameter of 5 cm. Let d x ,...,d M be the predicted 100*/), 100*pj h diameter percentiles from stages 1 or 2 and F the cumulative distribution function obtained from them with linear inter polation, where Pi,P2,-- Pm-i,Pm are the fixed values given in stage 1 (i.e., 0. 0.10 0.95,1.0). Furthermore, let k be the index of the smallest percentile above 5 cm. The modified algorithm adjusts only percen tiles dk ,...,dM . The adjusted combination of percentiles and stand variables is the combina tion that is compatible with and as close as possible to the original combination with re gard to a specified distance measure. The dis tance (Equation 6a) is defined as the weighted sum of squared deviations from the predicted percentiles and measured stand variables, where the weights are the inverse errors of the percentiles and stand variables. The compati bility requirement is included in the algorithm as a constraint (Equation 6d) in which the stem number based on the adjusted percentiles is equated to the adjusted stem number. Addi tional constraints (Equations 6b, 6c, 6e and 6f) are needed to ensure monotony of the adjusted percentiles and non-negativity of the percen tiles and stand variables. The adjusted set of percentiles and stand variables is found by solving the following optimization problem: ln( 5 , (6c) 4(Grf > 5 (^ +^)~ 5 lv 5 (^ +Si), - 100* „ r f \\ -^+s »" (6d) + g 1 1 ( 5 +sG >o, (6e) + s.y O, (6f) Optimizing field data collection 7 Figure 1. The map of sample plots and trees in one stand. The diameter of the mark is proportional to that of the tree. Data The primary dataset of the study This study used a permanent dataset col lected by the Finnish Forest Research Institute between the years 1976 and 1992 (Gustafsen et al. 1988). The stands were a sample from those stands of the 6 th and 7 th national forest invento ries that were located on mineral soils and on forest land. In each stand. 3 permanent circular plots were measured 1-3 times with 5-year intervals. However, only the first and third measurement occasions were used in this study. The plot size varied between stands and was determined so that at least 120 sample trees per stand were measured in southern Finland and 100 sample trees per stand in northern Finland. The tree species and the diameter at breast height were recorded from all trees belonging to the sample plot. In addition, the height was recorded from approximately one third of the trees, i.e. those trees which belong to the in nermost circle of the plot (see Figure 1). This study used only the Scots pine trees of the data. The dataset of this study included those stands where no regeneration cuts were made between the first and third measurement occa sions. the total volume of Scots pine was more than 10m 3 /ha, the total basal area of Scots pine was more than 4 m 2 /ha and a HPS plot using a basal area factor of 1 could be established on all plots of the first and third measurement occasions (see below). The total number of stands in the data was 170. The stand descrip tion was predicted for the third measurement occasion using new measurements simulated using the data of the third occasion and old measurements simulated using the data of the first occasion. Data preparation In order to calculate the (assumed) true stand volume and saw timber volume of a stand, the heights of all sample trees were needed. The unknown heights were predicted using the model of Mehtätalo (2004 c) (model V), which was localized for each stand using the meas ured height sample trees of the third measure ment occasion. The saw timber volume was defined as the total volume of those parts of stems that include at least a 4-meter long log with a minimum top diameter of 15 cm. In the calculation of volume, the taper curve func tions based on diameter and height were used (Laasasenaho 1982). The within-stand variances of DGM, basal area and stem number were needed in order to calculate the sampling errors of the stand measurements in the simulation. They were estimated from the measurements generated on the three sample plots of the stand. The gener ated measurements of basal area and DGM were made from HPS plots using a basal area factor of 1 and the measurements of stem number from a circular plot with a radius of 3.99 m. It was assumed that the true values of DGM and basal area could be obtained using all sample trees of the stand, i.e. it was as sumed that the expectations of the measure ments were known. The within-stand sampling variance of random variable X was calculated using the formula i[*,-£(x)] 2 . (7) 8 Mehtätalo and Kangas where E(X) is calculated from all sample trees of the stand and X f is the measurement gener ated on the z" 1 sample plot. The covariance of sampling errors of variables X and Y was cal culated correspondingly as The radius of the HPS plot determined by the thickest tree of the circular plot was often greater than the radius of the circular plot. In this case, all trees that would have belonged to the HPS plot were not present in the data, i.e. the circular plot should have been larger in order to include all trees of the HPS-plot. The area of the missing surface was calculated by subtracting the area of the circular sample plot. A f , from the area determined by the maximum diameter of the plot, D max , where A s is the area of the missing surface. To obtain unbiased measurements on the HPS plots, trees were generated on the surface. The stock density and diameter distribution on the surface were assumed to be similar to those on the fixed plot and the tree locations were as sumed to be random. These criteria were ful filled by duplicating trees of the circular plot on the missing surface and placing them at random locations as follows (see Figure 1). A random number (from L'(0.1)) was generated for each tree on the circular plot. If it was less than the ratio AJAf, the tree was included in the surface. If A s was greater than A,, for one or more plots in a stand, the stand was not in cluded in the data. The distances of the dupli cated trees from the center of the plot were generated using the formula + • where ''Dmax and iy are the radii corresponding to the areas in Equation 9 and C/(0,1) is a random number from the uniform distribution between 0 and 1. The result of the above calculations was a stand specific dataset including the true values of the stand characteristics and the within stand variances and correlations of basal area, DGM and stem number on the first and third measurement occasion. In addition, the data included the (assumed) true total and saw tim ber volumes of the third measurement occa sion. A short summary of the data is given in Table 1 and histograms of estimated within stand standard deviations and correlations in Figure 2. In addition to the stand specific dataset, data sets of potential new and old sample trees were generated. The dataset of potential new sample trees included those trees of the third meas urement occasion that belonged to the HPS plots established at the centers of the circular plots. However, only the true height sample trees were included, i.e. the duplicated trees and trees with unknown height were deleted after determination of the ranks and the total number of trees on the HPS plot. For each potential sample tree, diameter, height, rank on the HPS plot and total number of trees on the plot were saved. The dataset of potential old sample trees was generated in a similar manner from the data of the first measurement occa sion but only heights and diameters of the plot specific basal area median trees were included. Table 1. Summaries of the data. Note. DGM is the basal area weighted median di ameter of the stand. t[.Y-£(T)][K-£(}•)] cov(x,y) = . (8) A=*{D m j2) 2 -A f , (9) min. mean max. stand age, years 22 72 180 basal area, nr/ha 4.5 15.7 28.5 stem number. 1/lia 123 1052 3463 DGM, cm 6.3 16.8 32.0 total volume, nr/ha 19.7 107.1 308.8 saw timber volume, nrAia 0 47.4 254.4 Optimizing field data collection 9 Figure 2. Histograms of the estimated within-stand standard deviations of basal area (a). DGM (b). stem number (c) and the within-stand correlations of basal area and DGM (d) on the third measurement occasion. Construction of the data for the response surface models The next step in the analysis was to define the measurement strategies and simulate reali zations of them for each stand. These realiza tions were then used to study the accuracy of predicted stand description with different measurement strategies. The simulated inventory of a stand com prised various amounts of the following meas urements: 1. Measurement of a basal-area and DGM from an HPS plot using a basal area factor of 1 (G-plot). The basal area is then the number of trees on the plot and the basal area median diameter (DGM) of the plot is determined as the median diameter of the sample. 2. Measurement of the stem number by count ing the number of trees above 5 cm from a plot with a radius of 3.99 m (N-plot). The N plots were assumed to be from different lo cations than the G-plots (i.e. no correlation between the sampling errors of N- and G plots was assumed) in order to derive more information from the N-plots and to avoid problems in the simulation of different num bers of G- and N-plots. [see Husch et al. (1982. p. 258-260) for previous work on combining measurements from HPS-plots with measurements from fixed-area plots.] 3. Measurement of diameter and height of one of the trees on the G-plot (H-tree). 1 0 Mehtätalo and Kangas 4. Measurement of diameter and rank of one of the trees on the G-plot (Q-tree). 5. In addition, it was assumed that one 10-year old height measurement from the DGM-tree might be available from the stand (old H tree). A measurement strategy of the stand defines how many G-plots, N-plots, H-trees, Q-trees and old H-trees are measured from a stand. A total of 384 (4x3x4x4x2) measurement strate gies were constructed by varying systemati cally the number of measurements carried out in a stand between the assumed minimum and maximum values that would be used in prac tice. The number of G-plots was 1, 3, 5 or 7, the number of N-plots 0, 2 or 4. the number of H-trees 0, 2, 4 or 6, the number of Q-trees 0, 2, 4 or 6 and the number of old H-trees 0 or 1. When simulating a realization of a meas urement strategy, the measurements of H- and Q-trees were selected randomly from the set of potential new sample trees and the old H-tree measurement was selected randomly from the set of potential old sample trees of the stand. The H- and Q-trees were the same trees when ever possible, which minimized the number of diameter measurements in the inventory. Since the required amounts of G- and N plot meas urements could be more than three, the true sample plots were not used, a Monte Carlo approach was used instead. The measurements were generated by adding errors from a multi normal distribution to the true values of the stand characteristics. The plot measurements were assumed to be accurate and the generated measurements included only a sampling error component. The error variances were obtained by dividing the within-stand variance obtained with Equation 7 by the number of plots. The correlation between DGM and basal area measurements was the within-stand correlation based on Equation 8 and their correlation with stem number was zero since the measurements were assumed to be from different locations. When an old height measurement was used in the prediction of the H-D curve, old DGM and basal area measurements were needed in the calculation of the fixed part of the height model. The old measurements were generated in a similar manner as the new ones, but as suming that the number of sample plots in the previous inventory was 5, i.e. the variance of the generated measurements was the within stand variance of the first measurement occa sion divided by 5. The information about other stand variables (age, site fertility, information about thinning, stand coordinates, altitude, temperature sum) was assumed to be as accu rate as in the data of this study and no meas urement errors were added to them. Ten realizations of each measurement strat egy were produced for each stand and the diameter distribution and H-D curve were predicted using the system described previ ously. Using these predictions, the total vol ume and saw timber volume were calculated for each realization. The root mean squared errors (RMSE) of these characteristics were calculated from the 10 realizations for each stand and measurement strategy, thus resulting in 65280 (=l7O standsx3B4 strategies) RMSE values of total and saw timber volumes. The RMSE of variable V was calculated using the equation where Vt is the i h observation and V,nie is the true value of variable V. Modeling the errors The RMSE of volume and saw timber vol ume were modeled using the number of meas urements (G-plots. N-plots, H-trees. Q-trees and old H-trees) and stand characteristics (basal area and DGM) as predictors. The RMSE using measurement strategy j in stand k was assumed to follow the model where ykj is the RMSE, xkj includes the fixed predictors and b the fixed parameters. Ukuk is a random effect for stand k, v ; a random effect for strategy j and ekj the residual error for strat \w~ vJ RAISE , (10) = x i7b + "t +v j+v n i) Optimizing field data collection 11 egy j in stand k. The variance of the residual error, var(ej ;), seemed to increase as the pre dicted value increased and it was taken into account by using weights in the estimation of the model. The weights were obtained from the inverse of a variance function that was fitted using an approach similar to that of Lappi (1997). The mixed model approach was used to take into account that the observations of a stand are correlated, as are the observations of a strategy. The model was estimated using the MIXED procedure of SAS (Littell et al. 1996). Assuming that the sample plots and sample trees are a random sample from the population, their effect on the RMSE should be propor tional to l/ V number of measurements . Plot ting the observations against the number of measurements showed that this assumption holds rather well in the data and transforma tions of this form were used. First, a main effects model including only significant predictors (using 1% level of sig nificance) was fitted. Plotting the residuals in different groups of fixed predictors showed that the stand variables had significant interac tions with the amounts of measurements. In addition, the number of old and new height measurements interacted. The cross products were included in the model and the non significant ones were dropped one by one using backward elimination until all terms of the model were significant. The estimated models are in Table 2 and Figures 3 and 4 show predictions of the models using different measurement strategies. In the estimation of the model for RMSE of total volume (RV, otal), all 170 stands were used. Since the volume of saw timber was often zero in stands with a small DGM. the models for RMSE of saw timber volume (R Vsm, lmlber ) used only stands with a DGM of more than 13 cm. To provide an idea of the total variation in the data, a vari ance component model with the intercept as the only fixed predictor was fitted to the data. The estimated variance components of this model are shown in the last three lines of Table 2. However, note that the total variation of single realizations is greater than the variance component because the observations of the data are based on ten realizations. The estimated coefficients (Table 2) show that RV, olni and RV sau. ,imber decrease with in creasing amounts of measurements and de creasing basal area and DGM. The old height measurement decreases the RV totah except for situations where the number of new height measurements is high and the basal area is low, and it decreases the RV sm nmher if the number of new sample trees is three or less. These trends are realistic and in accordance with our prior hypotheses. The consistency of the model coefficients was considered by plotting the model predictions against DGM and basal area using different values of other predictors. The plots showed logical behavior with different values of stand variables and stand measure ments. Examples of these plots are shown in Figures 3 and 4. In addition to basal area and DGM, the most important factors affecting RVlola/ were the number of G-plots and H-trees (Figure 3). The number of old H-trees had a slight effect on the predicted R V,„ lai but the numbers of Q-trees and N-plots were not significant predictors at all. The number of N-plots did not affect RV sair nwber significantly either. However, the number of Q-trees was one of the most impor tant predictors in the model of RV smr , lm i,er. along with the numbers of G-plots and H-trees (Figure 4). Note that when stand DGM is small, the measurement of one Q-tree may decrease the RVsmv nmb„ even more than the measurement of seven G-plots. The effect of an old H-tree was, again, significant but slight. The main result with respect to the use of an old height sample tree is. however, that it re duces RVtola, (except for stands with low basal area and some new H-trees) and RV:au nmher (if the number of new sample trees is three or less). Thus, by using old height sample trees, an improvement on the predictions of total and saw timber volume can be obtained with little effort. 12 Mehtätalo and Kangas Table 2. Models of the root mean squared errors of volume (left) and saw timber volume (right). Note: Tlie estimated standard errors are in parentheses. The total variation is obtained from a model with the intercept as the only fixed predictor. The notations are as follows: /)G=number of G-plots; /;//=number of H trees; nHo=number of old H-trees (0 or 1); >/(>=number of Q-trees; G=basal area. m 2 /lia: Z)GA/=basal area median diameter, cm. RV total Vscnvtimber Fixed part Intercept -0.2195 (0.3805) -0.5369 (0.6725) 1/yfnG -42.333 (1.5905) -6.2779 (0.1468) xj-JnH + X -5.5433 (0.2565) -0.2364 (0.1581) l/yjnQ + l - 3.5707 (0.1394) nHo 0.7473 (0.1305) 0.3144(0.09612) nHoj V nH +1 -0.8203 (0.1647) -0.6574 (0.1495) G 2 - 0.008175 (0.001856) DGM/yfnG 1.9529 (0.0503) . G/V^G 0.29 (0.01683) - 1/(DGM*^//^G) 334.35 (11.8928) _ DGM 2 /4nG - 0.0373 (0.000409) G 2 /V^G - 0.003378 (0.00035) DGM/0. (12c) Optimizing field data collection 15 minimized, the resources are distributed be tween G-plots, H-trees and Q-trees (Table 3). The values of stand characteristics have a strong effect on how the resources should be distributed between different measurements. In particular, a slight growth in DGM causes remarkable changes in the optimal combina tion. In order to show an example where the cost vector and maximum measurement time de pend on the state vector, the optimization was carried out again using a cost vector of and c max of 0.25 XG, where Gis the basal area of the stand. The time unit is now the time required for the measurement of one G-plot in a stand with a basal area of 20 m 2 /ha. The solution differs from that in Table 3, but the trends are quite similar in both solutions. How ever, no conclusions about the optimal combi nation should be drawn on the basis of these results because of the ad hoc definitions of the cost vectors. If both total volume and saw timber volume should be predicted accurately, both variables need to be taken into account in searching for a single solution. In this case, the objective func tion can be defined as a weighted average of the RVtoiai and RVsaw umber (Table 5). Table 3. The solution of the optimization problem (Equations 12a-12c) using costs determined by vector c'=( 1,0.3,0.5) and a cmax of sin the four hypothetical stands. Table 4. The solution of the optimization problem (Equations 12a-12c) using costs determined by Equation (13) and a c maN of 0.25 XG in the four hypothetical stands. Table 5. The solution of the optimization problem (Equations 12a-12c) using the objective function r=0.2x RVtotal+ O.8xRVsa w. ,timber costs determined by Equation (13) and a cmax of 0.25xG in the four hypothetical stands. c'=(o.o5xG, O.I+O.OUDGM, 0.033xG) (13) target variable DGM G si S-, S3 v(s,r) R Vtotal 15 15 3.97 3.42 0 8.88 R'san' timber 15 15 2.48 2.92 3.29 5.57 RI total 25 15 4.02 3.27 0 15.6 R'saw timber 25 15 4.25 2.5 0 13.01 R V,total 15 25 3.68 4.39 0 12.15 R y saw timber 15 25 2.4 3.2 3.27 11.16 R Vtotal 25 25 3.84 3.88 0 18.92 RVsmr timber 25 25 3.88 2.62 0.67 19.12 target variable DGM G Sl s- > S3 vo and X2>o and the unbiased measurements Xt and X 2 are lognormally distributed with expectations X x and X 2 and variances var ( Ä- , j = cr, 2 . var lj(2 ) = cr 2 and co v^Xt ,X =p. We want to know £ (in X} j. var (in Ä,) and cov(lnX,,lnX2 j . It follows from the properties of a lognormal distribution (see. e.g.. Johnson and Kotz 1972: 18-20. Hines and Montgomery 1980: 188-191) that the logarithmic measurements are multi normally distributed with the expectations, variances and covariances of and Now assume that we have k unbiased log normally distributed measurements in random vector i = yX l ,Xl ,...,Xk j' with the variance covariance matrix var(x) =C. Formulas (A1)-(A3) give the matrix results and Prediction error of models with erroneous logarithmic predictors Assume a model of the form where ym is a vector of interest variables from observation m, B is a matrix of parameters. \ m the vector of predictors and e„, the vector of residual errors. For example, vector y m could include percentiles of a diameter distribution of stand m. Assume that the measured predictors are lognormally distributed with E (x m ) =\ m and var(x n , )= C. We want to know how errors in the predictors affect the model predic tions. y , =BIn x ~ v ill 111 The expectation of the prediction is Using Equation (A 4) we get £(lnA',) = 21nA', - jln(er, : + A', 2 ) . (Al) var(lnA',) = ln(af +X, 2 )-21n.Yl (A 2) = ln(/7+X 1X2 )-ln(J[' 1A'2 ) . (A3) £(ln i) = 21n x -jln{[l ®(C + xx')]lM } (A 4) var(lnx) = ln(C + xx')-ln(xx') . (A 5) y„, = Blnx m +e„,, (A 6) E(y m ) = £(Blnx„) = B£(lni,„) . (A 7) 20 Mehtätalo and Kangas Thus, the predictions are biased, the last term of (A 8) being the bias. The true value of the bias cannot be calculated since the expecta tions of the measurements are not known. However, approximate bias can be calculated by substituting the expectations of the meas urements with the measured values. An ap proximately unbiased prediction with model (A 6) using measurements (A 7) is obtained by subtracting the approximate bias from the prediction: The variance-covariance matrix of the error term is Using Equation (A 5) we get The first term of (A 10) is the increase of the residual variance caused by the measurement errors of the predictors. £(y,„)=Bf2lnx„,-^ln{[l®(C +^A;)]l}l f j v (A8) =Blnx„, +Bl lnx„, --ln{[l®(C+xmxm')]lj I -Bln(i) ~h{[i(E)(c+v,;)]l}j. (A 9) var(y,„) = var(Blni,„+e,„) = Bvar(lnx ; „)B' + D var(y,„) = B[ln(C + xx')-ln(xx')]B' + D « B[ln(C + xx')-ln(xx')] B'+D ISBN 951-40-1934-2 ISSN 0358-4283