Composite Estimation of Stand Tables
Dissertation
zur Erlangung des akademischen Grades Doctor of Philosophy (PhD) der Fakultät für Forstwissenschaften und Waldökologie der
Georg-August-Universität Göttingen
vorgelegt von Daniel Bierer
geboren in Kilchberg (CH)
Göttingen, 2008
2. Gutachter: Prof. Dr. Christoph Kleinn (Georg-August-Universität Göttingen) 3. Gutachter: Prof. Dr. George Gertner (Universität von Illinois)
Tag der mündlichen Prüfung: 6. März 2008
Abstract
One of the most important parameters for characterizing a forest stand is its stand table. In this thesis, the problem of estimating stand tables in stands with only a few sample plots (i.e. in small areas) is considered. An alternative estimation strategy to the common design-unbiased approach is proposed and it will be shown that drastic savings in mean squared error can be achieved by relaxing the restriction of unbiasedness.
In contrast to previous works on small area estimation of stand tables, a com- posite estimator is suggested that takes into account the mean squared errors instead of only the variances of its component estimators. As component estima- tors serve the common design-unbiased estimator and a synthetic estimator that links external data by kernel smoothing. The weighting scheme of the proposed composite estimator is based on a weight estimator for estimating Schaible’s approximation to the optimal constant weighting scheme.
The results of the thesis are promising in that the suggested composite esti- mator performs well also in stands for which the synthetic component estimator is highly biased. A further appealing property indicated by the results of the thesis is the general applicability of the proposed estimator, i.e. it performs well in mixed and unmixed stands as well as in even and uneven aged stands.
Furthermore, the diameter distribution of the subject stand has not to look like a presumed parametric distribution such as the Weibull distribution, i.e. it can have also, for instance, a multimodal distribution. Finally, the suggested esti- mator is also attractive in terms of simplicity: only (weighted) averages have to be build and no numerical methods for optimization or integrating have to be applied. All these properties make the estimator appropriate for the practical use in forestry.
Zusammenfassung
Einer der wichtigsten Parameter zur Charakterisierung eines Waldbestandes ist dessen Durchmesserverteilung. In der vorliegenden Dissertation wird das Problem der Schätzung der Durchmesserverteilung aus nur wenigen Stichproben (d.h. in Kleingebieten) betrachtet. Es wird eine alternative Schätzstrategie zum üblichen design-unverzerrten Ansatz vorgeschlagen und es wird gezeigt, dass der mittlere quadratische Fehler drastisch reduziert werden kann, wenn man die Ein- schränkung auf Unverzerrtheit lockert.
Im Gegensatz zu früheren Arbeiten zur Kleingebietsschätzung der Durchmesser- verteilung von Waldbeständen wird ein zusammengesetzter Schätzer vorgeschla- gen, welcher die mittleren quadratischen Fehler und nicht nur die Varianzen der Komponentenschätzer berücksichtigt. Als Komponentenschätzer dienen dabei der übliche design-unverzerrte Schätzer und ein synthetischer Schätzer, welcher externe Daten über Kern-Regression einbindet. Das Gewichtungsschema des vorgeschlagenen zusammengesetzten Schätzers basiert dabei auf einem Gewichts- schätzer zur Schätzung von Schaible’s Approximation des optimalen konstanten Gewichtungsschemas.
Die Resultate der Dissertation sind insofern vielversprechend, als dass der vorge- schlagene zusammengesetzte Schätzer auch in Beständen gut abschneidet, in welchen der synthetische Komponentenschätzer stark verzerrt ist. Eine weit- ere attraktive Eigenschaft des vorgeschlagenen Schätzers ist dessen allgemeine Anwendbarkeit, d.h. der Schätzer schneidet in Rein- und Mischbeständen wie auch in gleich- und ungleichaltrigen Beständen gut ab. Desweiteren wird für den Zielbestand keine bestimmte parametrische Verteilung wie z.B. eine Weibull- verteilung vorausgesetzt, d.h. der Zielbestand kann beispielsweise auch eine multimodale Durchmesserverteilung aufweisen. Der vorgeschlagene Schätzer ist zudem auch aufgrund seiner Einfachheit attraktiv: es müssen nur gewichtete Mittel gebildet werden und es werden keine numerischen Methoden zur Op- timierung oder Integration verwendet. Alle diese Eigenschaften machen den Schätzer für den Einsatz in der Forstpraxis geeignet.
Acknowledgements
I would like to thank to Prof. Dr. Joachim Saborowski for supervising this thesis with great interest and for making possible my stay in Göttingen, I will never forget the friendly atmosphere at the Institute of Forest Biometrics and Informatics at the Georg-August-University.
Further thanks go to my supervisor at the Swiss Federal Research Institute of Forest, Snow and Landscape Research in Birmensdorf, Dr. Adrian Lanz, for being always available when I needed his help or advise and for the many in- teresting and inspiring discussions. Moreover, I want to express my thanks to Prof. Dr. Christoph Kleinn (Georg-August-University Göttingen) for acting as co-referee and examiner, to Prof. Dr. George Gertner (University of Illinois) for acting as co-referee and to Prof. Dr. Jürgen Nagel (Northwest German Forest Research Institute) for acting as examiner.
I would also like to thank to Dr. Peter Brassel and Dr. Oliver Thees (both from the Swiss Federal Research Institute for Forest, Snow and Landscape Re- search) for the financial support and to all my other colleagues from Birmens- dorf, especially to my office colleagues Edgar Kaufmann and Jürgen Böhl for their friendship and to Dr. Norbert Kräuchi and Dr. Matthias Dobbertin for making available the test stand data.
Finally, I would like to thank to my parents, to my brother and especially to my wife Angélica, for the great personal support and encouragement.
”A combination of the two is better than either taken alone.”
(Royall, 1973)
List of Figures 10
List of Tables 12
1 Introduction 14
1.1 Problem Description . . . 14
1.2 Indirect Estimation Approach . . . 14
1.3 Differences to Previous Works . . . 15
1.4 Outline of the Work . . . 16
1.5 Notational Conventions . . . 16
2 Sample Designs 18 2.1 Basic Terminology . . . 18
2.2 Target Population . . . 18
2.3 Sampling Frames . . . 19
2.4 Design Stage . . . 20
2.5 Estimation Stage . . . 21
3 Principle of Composite Estimation 24 3.1 Component Estimators . . . 24
3.2 Combining Estimators . . . 25
3.3 Optimal Composite Estimator . . . 25
3.4 Schaible’s Approximation . . . 27
4 Direct Component Estimator 30 4.1 Unbiased Estimation . . . 30
4.2 Horvitz-Thompson Estimator . . . 31
4.3 Infinite Population Approach . . . 35
4.4 Different Views . . . 36
5 Synthetic Component Estimator 38 5.1 Synthetic Estimation . . . 38
5.1.1 Borrowing Strength . . . 38
5.1.2 Estimation Procedure . . . 38
5.1.3 Neighborhoods . . . 39
5.2 Recursive Partitioning . . . 41
5.2.1 Disjoint Neighborhoods . . . 41
5.2.2 Partitioning Algorithm . . . 43
5.2.3 Right-Sized Partitions . . . 44 9
CONTENTS 10
5.3 Kernel Smoothing . . . 50
5.3.1 Overlapping Neighborhoods . . . 50
5.3.2 Nadaraya-Watson Estimator . . . 52
5.3.3 Method of Fukunaga . . . 52
5.3.4 Nearest-Neighbor Bandwidths . . . 54
5.3.5 Selecting Neighborhood Size . . . 54
5.4 Model Selection . . . 55
5.4.1 Measure of Performance . . . 55
5.4.2 Double-Layered Cross-Validation . . . 56
5.5 Monte-Carlo Experiments . . . 60
5.5.1 Mean Relative Errors . . . 60
5.5.2 Bias-Variance Decomposition . . . 65
6 Composite Estimators 68 6.1 Principle of Selecting Weights . . . 68
6.1.1 Optimal Weights . . . 68
6.1.2 Estimating Mean Squared Errors . . . 68
6.2 Different Approaches . . . 69
6.2.1 Approximation Approach . . . 69
6.2.2 Ordinary Approach . . . 72
6.3 Monte-Carlo Experiments . . . 75
6.3.1 Mean Relative Errors . . . 75
6.3.2 Bias-Variance Decomposition . . . 79
7 Discussion 84 7.1 Direct Component Estimator . . . 84
7.2 Synthetic Component Estimator . . . 86
7.3 Composite Estimators . . . 87
7.4 Comparison with Previous Works . . . 88
7.4.1 Linking Data at the Stand Level . . . 88
7.4.2 Precision-Weighted Composite Estimation . . . 88
7.4.3 Multinomial-Dirichlet Approach . . . 89
7.5 Applicability in Forest Practice . . . 92
A Appendix 93 A.1 Description of the Data Sets . . . 93
A.1.1 Test Stands . . . 93
A.1.2 Learning Data Set . . . 101
A.2 Percentile Intervals . . . 103
A.2.1 Synthetic Estimator . . . 104
A.2.2 Composite Estimators . . . 108
A.3 Mahalanobis Transformation . . . 116
A.3.1 Jordan Decomposition . . . 116
A.3.2 Centering Matrix . . . 116
A.3.3 Standardization and Decorrelation . . . 117
Bibliography 119
Curriculum Vitae 125
3.1 Normalized mean squared errors (with respect to the mean squared error of the synthetic estimator) as a function ofφ whenF = 1 (after Schaible (1978b)). . . 29 3.2 Normalized mean squared errors (with respect to the mean squared
error of the synthetic estimator) as a function ofφ whenF = 2 (after Schaible (1978b)). . . 29 3.3 Normalized mean squared errors (with respect to the mean squared
error of the synthetic estimator) as a function ofφ whenF = 6 (after Schaible (1978b)). . . 29 5.1 Disjoint neighborhoods (after Hastie et al. (2001)) . . . 41 5.2 A regression tree (after Breiman et al. (1984)). . . 44 5.3 Principle of 10-fold cross-validation (after Hastie et al. (2001)) . 47 5.4 Estimated relative generalization errors (with standard errors) as
function of tree size and application of the 1 SE pruning rule. . . 48 5.5 Optimal partitioning of the learning data set by using the 1-
SE rule of Breiman et al. (1984) (NPH=number of stems per hectare, BASFPH=basal area per hectare, EST=stage of devel- opment (see Appendix A.1.2)). . . 49 5.6 Disjoint and overlapping neighborhoods . . . 50 5.7 Side-length of a subcube needed to capture a fraction of the vol-
ume of the data for different dimensions. In ten dimensions 80%
of the range of each coordinate is needed to capture 10% of the data (after Hastie et al. (2001)). . . 51 5.8 Estimated relative generalization error as function of the nearest-
neighbor bandwidth for different methods of transforming inputs previously to smoothing (estimates are connected by dashed lines). 55 5.9 Folds for model- and parameter selection by double-layered cross-
validation (layer L I contains the folds for model selection and layer L II the folds for parameter selection). . . 56 5.10 Estimated generalization error curves for recursive partitioning
based on the learning data setsL(−1),L(−2), ....,L(−10). . . 57 5.11 Estimated generalization errors of the Nadaraya-Watson estima-
tor based on the standard method for the learning data setsL(−1), L(−2), ....,L(−10). . . 58 5.12 Estimated generalization errors of the Nadaraya-Watson estima-
tor based on the Fukunaga method for the learning data sets L(−1),L(−2), ....,L(−10). . . 58
11
LIST OF FIGURES 12
5.13 Comparing performances ofθˆdirandθˆsyn. . . 62
5.14 Stand A . . . 62
5.15 Stand B . . . 62
5.16 Stand C . . . 64
5.17 Stand D . . . 64
5.18 Comparing bias-variance decompositions ofθˆdir(D) andθˆsyn (S). 67 6.1 Comparing performances ofθˆcom1andθˆcom2. . . 78
6.2 Comparing bias-variance decompositions ofθˆcom1 andθˆcom2. . . 83
A.1 Graphical representation of the stand tables . . . 100
A.2 P90% intervals ofθˆ(k)syn(k=1..20) in stand ”Lausanne” . . . 105
A.3 P90% intervals ofθˆ(k)syn(k=1..20) in stand ”Neunkirch” . . . 107
A.4 P90% intervals ofθˆ(k)com1 (k=1..20) in stand ”Lausanne” . . . 109
A.5 P90% intervals ofθˆ(k)com1 (k=1..20) in stand ”Neunkirch” . . . 111
A.6 P90% intervals ofθˆ(k)com2 (k=1..20) in stand ”Lausanne” . . . 113
A.7 P90% intervals ofθˆ(k)com2 (k=1..20) in stand ”Neunkirch” . . . 115
5.1 Reduction of the RSS by the two ancestor nodes relative to the RSS of the parent node LRRR provided by the two best splits. . 49 5.2 Bandwidths that minimizes the estimated generalization error. . 55 5.3 Selected model parameter values . . . 59 5.4 Estimated RGE for the competing regression methods . . . 59
13
Chapter 1
Introduction
1.1 Problem Description
Advanced forest management planning often depends on detailed information at the stand level (v. Gadow and Schmidt, 1998). One of the most important characteristics of a forest stand in forest management planning is its stand table which is a table that contains the number of trees per unit area by diameter class (see Green and Clutter (2002) and van Laar and Akca (2007)). Stand tables are needed for several purposes such as:
• Use of growth models
• Ecological assessment of stands
• Calculation of timber assortments
Usually, stand tables are estimated by direct estimators, i.e. estimators which are based exclusively on data from the stand of interest (Green and Clutter, 2002). The problem of direct estimation is that the corresponding estimators can be highly imprecise when the class width of the stand table is chosen very small or in stands with only a few sample plots. Under such conditions, the stand is regarded as a small area. The term ”small area” is generally used when direct estimates of adequate precision can not be produced, otherwise, i.e. when direct estimates of adequate precision can be produced, the stand is regarded as a large area (see Rao (2003))1. In the scope of this thesis, indirect estimation (see next section) is proposed as an alternative estimation strategy to estimate stand tables in small areas. The work can be seen as a contribution to the implementation of a stand-oriented inventory system as it was suggested by v. Gadow and Schmidt (1998).
1.2 Indirect Estimation Approach
In this thesis, indirect estimation is proposed as estimation strategy in stands with only a few sample plots. In contrast to a direct estimator, an indirect estimator is an estimator that uses (exclusively or not exclusively) values of
1For the sake of clarity, it should be noted that accuracy has to be adequate as well.
14
the variable of interest from stands other than the stand of interest and thus increase the ”effective” sample size. These external values are brought into the estimation process through a linking model that provides a link to sample plots in stands that are similar to the subject stand through the use of auxiliary information related to the stand table, i.e. availability of good auxiliary infor- mation and determination of a suitable linking model (see Chapter 5) is crucial to the formation of an indirect estimator (Rao, 2003). The indirect estimator constructed in this thesis is a combination of a design-unbiased direct estimator (Chapter 4) and a synthetic estimator (Chapter 5) which is an estimator that estimates the stand table of the subject stand by synthesizing the information from many stands that are similar to the subject stand. The two estimators are combined by building a weighted average of the estimates that result from the two estimators for each realization. The resulting estimator, i.e. the estimator whose realizations are the weighted averages, is called a composite estimator (Chapter 3). For each realization, the performance of the two estimators is estimated from the observed data in order to give more weight to the estimate which results from the estimator whose estimated performance is better, i.e. not only the two estimates are random but also the weights of the weighted averages (Chapter 6).
1.3 Differences to Previous Works
To the author’s knowledge, the works of Green and Clutter (2000) and Green and Clutter (2002) are the only two works which have been dedicated to com- posite estimation of stand tables previously to this work. The methodology sug- gested in this thesis differs from these previous works by the weighting scheme which the composite estimator is based on and by the approach for constructing the synthetic estimator:
• The main difference of the weighting scheme suggested in this thesis to the weighting scheme of the previous works is that the mean squared errors instead of only the variances of the component estimators are taken into account. More concretely, the weights of the component estimates are determined by estimating Schaible’s approximation (Schaible, 1978b) to the optimal constant weighting scheme from the observed data.
• The proposed synthetic estimator incorporates external information into the estimation process by kernel smoothing instead of by the traditional partitioning approach (see Drew et al. (1982)). The corresponding linking model is defined at the plot level and not at the stand level as in the previous works. A linking model which is defined at the plot level is helpful especially in heterogenous stands (see Section 7.4.1).
As in the previous works, Monte-Carlo experiments are used for assessing the performance of the different competing estimators. However, in contrast to the previous works, real data sets are used instead of simulated data sets. That is, Monte Carlo experiments are based on real test stands that comprise a wide variety of stand types (see Appendix A.1.1) and sample plot data from the Swiss National Forest Inventory (see Appendix A.1.2) is used as external data source.
CHAPTER 1. INTRODUCTION 16
1.4 Outline of the Work
The thesis is structured into the following seven chapters:
• Chapter 1: In the first chapter, the problem of estimating stand tables in stands with only a few sample plots is formulated. Moreover, differences between this thesis and previous works are emphasized and an overview of the work is given.
• Chapter 2: In the second chapter, the basic terminology of sampling theory, the formalism for characterizing the problem and the selection criteria for competing estimators are introduced.
• Chapter 3: In the third chapter, the principle of composite estimation is described and the optimal (constant) weights for composite estimation are derived. Moreover, Schaible’s approximation to the optimal weights and the properties of the corresponding composite estimator are outlined.
• Chapter 4: In the forth chapter, the design-unbiased direct estimator is outlined under the classical finite population approach as well as under the more modern infinite population approach.
• Chapter 5: In the fifth chapter, a synthetic estimator for estimating stand tables is introduced. It is shown that this estimator has a smaller variance then the direct estimator when sample size is small. However, it is also shown that this estimator can have a considerable bias and that it is not recommendable to use it without correction by a direct estimator when a direct estimator is available.
• Chapter 6: In the sixth chapter, two different approaches for estimating Schaible’s approximation to the optimal constant weighting scheme are considered: the approximation approach and the ordinary approach. It is shown that the composite estimator based on the ordinary approach should be preferred.
• Chapter 7: In the seventh chapter, conclusions from the results of thesis are drawn and future fields of investigation are pointed out.
Finally, the appendix contains some additional material such as the description of the used data sets, additional figures to the Monte-Carlo experiments and technical details to the Mahalanobis transformation.
1.5 Notational Conventions
In this thesis, the following notational conventions, which can differ slightly from some other common conventions, have been adopted:
• Scalars: Constant or variable scalars are denoted by upper- or lower-case letters that are not bold. If a scalar is denoted by a single letter, the letter is italic, such as for exampleE(ˆθ(k)). Otherwise, the letters are not italic, such as for example var(ˆθ(k))or MSE(ˆθ(k)). If the scalar is followed by parenthesis with an argument (or several arguments), either a single value (for a certain argument value) or the corresponding function is meant.
• Vectors: Vectors are denoted by letters that are lower-case and bold. In this thesis, vectors are always column vectors, that is
x=
x(1), x(2), ..., x(p)>
It should be noted that the dimensions are denoted by superscripts in brackets because vectors often have additional subscripts in this thesis.
• Matrices: In order to differ between vectors and matrices, matrices are denoted by letters that are upper-case and bold:
X=
x(11) · · · x(1p) ... . .. ... x(p1) · · · x(pp)
• Random Variables: A common convention is to denote random vari- ables by upper-case letters and their values by lower-case letters. If this convention would be adopted, it would follow thatX denotes both a ran- dom vector and a matrix, and this would in turn become a source of confusion. Therefore, it seems preferable to use the lower-case letter also as the symbol for a multivariate random variable because it is usually clear from the context if the random variables or their values are described (see Krzanowski (2000)).
• Sets: Sets are denoted either by calligraphic symbols or by specifying their elements within curly brackets. For instance,
E ={e1, e2, e3, ..., ene}
denotes the set of population elements, i.e. the trees in the forest stand of interest. It should be noted that in fact every function (such as var(ˆθ(k))) can be regarded as a set of pairs of numbers (see Laugwitz (1996)). How- ever, this view is not adopted in the scope of this thesis, i.e. var(ˆθ(k)) is not written with calligraphic symbols.
Apart from the notation for sets, the adopted conventions correspond to the conventions used in Krzanowski (2000).
Chapter 2
Sample Designs
2.1 Basic Terminology
Two steps have to be conducted in order to estimate the stand table of a subject stand: in the first step data has to be collected in the subject stand (based on a sampling plan) and in the second step this data has to be used as input of an estimation procedure which produces a stand table estimate, where in the scope of this thesis additionally input data from outside of the subject stand is processed by the estimation procedure. The term design stage is used to designate the period during which the sample selection procedure is decided and the sample selected (see Section 2.4). By contrast, the estimation stage (to be optimized in the scope of this thesis) refers to the period when the data is already collected and the required estimates are calculated (see Section 2.5).
These two steps are summarized in the sample design, i.e. the sample design consists of the sampling plan and the method of estimation (see Hansen et al.
(1953) and Dalenius (1988)). A synonymous term is sampling strategy (see Eriksson (1995)). By contrast, the similar term ”sampling design” characterizes only the sampling plan, i.e. how sample units are selected. Given a random sampling plan, it is possible (although not always simple) to state the probability P(s)of selecting a specified samplesfrom the target population. The function P(·)is then usually called the sampling design. It plays a central role because it determines the essential statistical properties (sampling distribution, expected value, variance) of random quantities calculated from a sample, such as the sample mean, the sample median, and the sample variance (Särndal et al., 2003).
2.2 Target Population
In the context of this thesis, the target population consists of a finite set ofne trees (elements) E ={e1, e2, ..., ene} in a forest stand that is given by its area A ⊂ R2 (a closed and bounded subset ofR2). Each tree (element) is defined by the triple (i,pi,yi), i.e. by its labeli∈ {1,2, ..., ne}, by its position vector
pi=
p(1)i , p(2)i >
∈ R2
18
(the set{p1,p2, ...,pn
e} is a subset ofA) and by its attribute vector yi=
yi(1), y(2)i , ..., y(ni a)>
which contains the values of the na variables of interest (such as diameter at breast height, tree height or volume). Within the scope of this thesis only the diameter at breast height is of interest, so that the vector yi reduces to the scalar yi, that is, yi means always the diameter at breast height of tree i in this thesis. The population parameter to be estimated is the stand table θ = (θ(1), ..., θ(c))> ∈ Rc which is given by the number of stems per unit area by diameter classk
θ(k)= Pne
i=1I(yi ∈ I(k))
λ(A) ,
whereλ(A)denotes the surface area ofAand I(yi∈ I(k)) =
1ifyi is in classk 0ifyi is not in classk .
In the scope of this thesis, the most common stand table specification of Switzer- land is used which is the classification by intervalsI(k)of lengthh=4 cm starting at 12 cm (see for example Mahrer (1988), Schmid-Haas et al. (1993), Brassel and Brändli (1999), Brassel and Lischke (2001) and Keller (2005)):
I(k) =
[(k+ 2)h,(k+ 3)h) fork∈ {1,2, ...,19}
[(k+ 2)h,∞) fork= 20
2.3 Sampling Frames
Carrying out a stand inventory calls for having some device of getting observa- tional access to the trees. In a sample survey, such a device is referred to as a sampling frame (Dalenius, 1988) and is given by the list of all units from which the sample is drawn. The entities that make up the population are called ele- ments and the entities of the frame are called the units or sampling units. If it were possible to disregard the cost of collecting data, it would be natural to pre- pare a list of all trees of the subject stand as the sampling frame and to sample directly from that list. Such a procedure is called direct element sampling be- cause the frame directly identifies the population elements (Särndal et al., 2003).
However, direct element sampling would be prohibitively expensive because of the costs of preparing the list, for this reason (random) selection of elements from the population has to be replaced by (random) selection of sampling units from a sampling frame, and some association rule linking population elements with the sampling units has to be specified (Dalenius, 1988). In forestry, the sampling units are typically the points inAfrom which observational access to the population elements (trees) can be achieved by corresponding association rules (see next section).
CHAPTER 2. SAMPLE DESIGNS 20
2.4 Design Stage
It is assumed throughout this thesis that sample trees from the subject stand are selected via ns circular sample plots of 300 m2 during design stage1. The centers of the sample plots are called sample points and are realizations of independent and identically distributed (i.i.d.) random variables zj with a uniform distribution. Such a sampling design can be formalized in the following way:
• The probability density functionf(zj)ofzj is given by f(zj(1), zj(2)) =
0 for zj ∈ A/ 1/λ(A) for zj ∈ A
• The probability that a sample pointzj falls in a certain subareaAs⊂ A of the stand, is given by
P(zj ∈ As) = Z Z
(z(1)j ,zj(2))∈As
f(z(1)j , zj(2))dzj(1)dz(2)j
= Z Z
(z(1)j ,zj(2))∈As
1
λ(A)dz(1)j dzj(2)
• In order to characterize tree selection, the following indicator function ι(pi,zj) = I((pi−zj)>(pi−zj)≤r2)
= I(kpi−zj)k ≤r)
is used (see Eriksson (1995)), where the position vector of thejth sample point (center of the sample plot) is given by zj = (z(1)j , zj(2))>. In the scope of this thesis, circular sample plots of 300 m2 are assumed so that the fixed radius of the sample plots is given byr≈9.77meters.
• The tree with label i (see Section 2.2) is selected at sample point zj if ι(pi,zj) = 1, i.e. ifzj is element of the inclusion area
Ai=C(pi, r)∩ A={q:q∈ A ∧ kpi−qk ≤r}
of the tree, whereC(pi, r)is the circular area around the tree with radius r:
C(pi, r) ={q:kpi−qk ≤r}
• The setEzj of all sample trees included at sample pointzj is given by Ezj ={ei|ei∈ E ∧ι(pi,zj) = 1}
and the setEsof all sample trees of the inventory is given by Es=
ns
[
j=1
Ezj =
ns
[
j=1
{ei|ei∈ E ∧ι(pi,zj) = 1} ⊂ E assuming a stand inventory withnssample pointszj.
1Such plots are standard sample plots in Switzerland, see Schmid-Haas et al. (1993).
• The values of the variable of interest (diameter at breast height) are ob- served on all trees which are element ofEs. The observed values together with the surface areas λ(C(pi, r)∩ A)of the inclusion areas C(pi, r)∩ A will be used during the estimation stage to estimate the unknown stand table (see Chapter 4 for details).
• In Chapter 5, additional sample plots from outside of the subject stand are brought into the estimation process. These additional sample plots originate from the Swiss National Forest Inventory which is based on con- centric plot sampling (see Appendix A.1.2), i.e. the radii of the above mentioned inclusion areas C(pi, r) around the trees ei depend on their diametersyi so that the corresponding indicator function is given by
ι(pi, yi,zj) =
I((pi−zj)>(pi−zj)≤r21)foryi< d I((pi−zj)>(pi−zj)≤r22)foryi≥d ,
where d=36cm,r1=7.98m (surface area of 200m2) andr2=12.62m (surface area of500m2). Furthermore, the sample pointszj are not realizations of a uniformly distributed random variable but build a regular grid of 1.41 km x 1.41 km over Switzerland (see Brassel and Lischke (2001)).
2.5 Estimation Stage
At the estimation stage, the data collected in the subject stand and eventually some additional data from other stands (that are similar to the subject stand) are used for producing an estimate of the unknown stand table. More precisely, these data are used as argument of a function, which is called an estimator θ,ˆ for calculating an estimate of the unknown stand tableθ. If the sampling design is a random design, i.e. when the stand inventory is based on sample plots that are allocated randomly in the subject stand, the estimator is a random variable as well because its argument, i.e. the collected data, is random. That is, the collected data and the resulting estimates differ from inventory to inventory if the inventory is repeated several times. It should be noted that only the de- sign variance is considered, i.e. the variance which is induced by the random sampling design. Apart from the design variance, there are at least two other concepts of variability in survey sampling (see Wolter (1985)): the population parameter to be estimated could be regarded as a realization of a random vari- able (and not as a constant as in this thesis) and the measured diameters could be assumed to have a random measurement error. However, these concepts are not considered in this thesis.
The properties of the three competing estimatorsθˆdir(Chapter 4),θˆsyn(Chap- ter 5) and θˆcom (Chapter 6) will be considered within the scope of this thesis.
Usually, the estimator whose distribution is most narrowly concentrated around the unknown population parameterθis wanted to be used because an estimator that varies little around the unknown value θ is on intuitive grounds ”better”
than one that varies a great deal (Särndal et al., 2003). This suggests using the
CHAPTER 2. SAMPLE DESIGNS 22
criterion of ”small mean squared error” because if MSE(θ)ˆ = E((θˆ−θ)>(θˆ−θ)) =Ekθˆ−θk2
= E
c
X
k=1
(ˆθ(k)−θ(k))2
!
=
c
X
k=1
E
(ˆθ(k)−θ(k))2
=
c
X
k=1
MSE(ˆθ(k)) is small, there is strong reason to believe that an actual realizationθˆis near the true population parameter (Särndal et al., 2003). However, even if the sampling distribution is tightly concentrated aroundθthere is always a small possibility that the particular sample was ”bad”, so that the estimate falls in one of the tails of the distribution, rather far removed from θ(Särndal et al., 2003).
It is important to note that the mean squared error of an estimator can be decomposed in bias and variance in the following way:
MSE(θ)ˆ = Var(θ) + (Bias(ˆ θ))ˆ 2
= E
kθˆ−E(θ)kˆ 2
+kE(θ)ˆ −θk2
= E
c
X
k=1
(ˆθ(k)−Eθˆ(k))2+
c
X
k=1
(E(ˆθ(k))−θ(k))2
=
c
X
k=1
E(ˆθ(k)−Eθˆ(k))2+
c
X
k=1
(E(ˆθ(k))−θ(k))2
=
c
X
k=1
Var(ˆθ(k)) +
c
X
k=1
(Bias(ˆθ(k)))2
Direct estimators are usually unbiased. However, they can have a high variance (can be highly imprecise), in particular when the class width of the stand table is chosen very small or in stands with only a few sample plots. Under such conditions, indirect estimators often perform better. However, indirect estima- tors can have a considerable bias (can be highly inaccurate) if they are based exclusively on data from outside of the subject stand. In this context, it should be noted that the statement that an estimator θˆis biased is a statement of average performance, namely, over all possible samples. However, to say that an estimate is biased is strictly speaking incorrect (Särndal et al., 2003). An estimate is a constant value obtained for a particular realization. This value can be off the mark, in the sense of deviating from the unknown parameter value θ. Because an estimate is a number, it has no variation and no bias (Särndal et al., 2003).
It should be noted that the mean squared relative error MSRE(θ)ˆ = E
kθˆ−θk kθk
!2
=E kθˆ−θk2 kθk2
!
= E
Pc
k=1(ˆθ(k)−θ(k))2 Pc
k=1(θ(k))2
!
is used instead of the mean squared absolute error in Section 5.5 and 6.3 be- cause relative errors reveal the inherent error characteristics of an estimator
better than absolute errors. For example, it is usually reasonable to expect that the relative error of an estimator is less variant than the absolute error as the magnitude of the parameterθ varies (Li and Zhao, 2005). The decomposition in bias and variance can be made in the same way by using the corresponding relative quantities:
MSRE(θ)ˆ = E
kθˆ−E(θ)kˆ kθk
!2
+ kE(θ)ˆ −θk kθk
!2
= E
Pc
k=1(ˆθ(k)−θˆ(k))2 Pc
k=1(θ(k))2
! +
Pc
k=1(E(ˆθ(k))−θ(k))2 Pc
k=1(θ(k))2
!
Apart from the purpose of bias-variance decomposition, the root mean squared relative error
RMSRE(θ) =ˆ E kθˆ−θk2 kθk2
!!12
= E
Pc
k=1(ˆθ(k)−θ(k))2 Pc
k=1(θ(k))2
!!12
and the mean relative error MRE(θ) =ˆ E kθˆ−θk
kθk
!
=E (Pc
k=1(ˆθ(k)−θ(k))2)12 (Pc
k=1(θ(k))2)12
!
are more common measures for evaluating estimators. In this thesis, the mean relative error is evaluated instead of the root mean squared relative error because of its better interpretability (see Li and Zhao (2001) and Li and Zhao (2005)).
Chapter 3
Principle of Composite Estimation
3.1 Component Estimators
In small areas, direct estimators are not precise enough to be of use by it- self. However, they can be used to reduce the bias (and, ultimately, the mean squared error) of indirect estimators. This can be achieved by forming a com- posite estimator. In general, a composite estimator combines two estimators, called component estimators, together resulting in an estimator which may be more accurate than either of its component estimators (Schaible, 1996). In this chapter, the principle of composite estimation is outlined. This principle is not new and has frequently been applied in small area estimation (see for instance Rao (2003)). To the author’s knowledge, it was mentioned for the first time by the U.S. National Center for Health Statistics (1968).
The two component estimators that are used in this thesis for forming a com- posite estimator can be characterized in the following way:
• Direct Component Estimator: The direct component estimator is an estimator which is based exclusively on data from the subject stand. An attractive property of the direct component estimator is that it can be con- structed in such a way that it is design-unbiased (see Chapter 4). However, the direct component estimator has usually a high variance when sample size is small (see (Särndal et al., 2003)).
• Synthetic Component Estimator: The synthetic component estima- tor is constructed as a (weighted) mean of values of the study variable observed outside of the subject stand in similar stands. These values are brought into the estimation process through a linking model that provides a link to similar sample plots (see Chapter 5). In small areas, the variance of the synthetic component estimator can be expected to be lower than that of the direct component estimator because the mean is built over many values, however, it can have a considerable bias (see Särndal et al.
(2003)).
24
It is evident that at some point, as the sample size in a small area increases, a direct estimator becomes more desirable than a synthetic one. This is true whether or not the sample was designed to produce estimates for small areas (Schaible, 1978b). Gonzalez and Wakesberg (1973) and Schaible et al. (1977a) compared errors of synthetic and direct estimates for Standard Metropolitan Statistical Areas and counties. The authors of both papers observed that the direct estimator was outperformed by the synthetic estimator when small area sample sizes were relatively small. By contrast, the synthetic estimator was outperformed by the direct estimator when small area sample sizes were large.
3.2 Combining Estimators
Royall (1973) in a discussion of papers by Gonzalez (1973) and Ericksen (1973), suggested that a choice between direct and synthetic approaches need not to be made but that ”... a combination of the two is better than either taken alone.”
(Schaible, 1978b). A natural way to balance the bias of a synthetic estimator θˆsynagainst the instability of a direct estimatorθˆdiris to take a weighted average ofθˆdirand θˆsyn:
θˆcom=φθˆdir+ (1−φ)θˆsyn
for a suitably chosen weightφ(0≤φ≤1). Such estimators are called composite estimators (see for example Schaible (1978b) and Rao (2003)). Within the scope of this thesis the approximation of Schaible (1978b) to the minimum mean squared error weighting scheme is suggested. An interesting property of composite estimators is that they can have a smaller mean squared error than either component estimator (see Schaible (1978b) and Section 3.4). Another characteristic of composite estimators with weights restricted to 0 ≤ φ ≤ 1 is that the mean squared error of the composite estimator is smaller than the larger of the two mean squared errors of the component estimators regardless of the weighting scheme used (see Schaible (1978b)). Typically, considerable weight is given to the synthetic estimates when sample size within the small area is small, and, as sample size is increased, the weight is gradually shifted onto the design-unbiased direct estimates (see Särndal et al. (2003)).
3.3 Optimal Composite Estimator
The optimal composite estimator θˆcom is based on a weighting scheme which minimizes the mean squared error of the composite estimator (Schaible, 1978b).
The mean squared error is given by
MSE(θˆcom) = φ2MSE(θˆdir) + (1−φ)2MSE(θˆsyn) + 2φ(1−φ)E((θˆdir−θ)>(θˆsyn−θ))
CHAPTER 3. PRINCIPLE OF COMPOSITE ESTIMATION 26
because the above term can be written as
MSE(θˆcom) = φ2E(kθˆdir−θk2) + (1−φ)2E(kθˆsyn−θk2) + 2φ(1−φ)E((θˆdir−θ)>(θˆsyn−θ))
= φ2E
c
X
k=1
(ˆθ(k)dir−θ(k))2
!
+ (1−φ)2E
c
X
k=1
(ˆθsyn(k)−θ(k))2
! +
2φ(1−φ)E
c
X
k=1
(ˆθ(k)dir−θ(k))(ˆθsyn(k)−θ(k))
!
= φ2
c
X
k=1
E
(ˆθ(k)dir)2−2ˆθ(k)dirθ(k)+ (θ(k))2 +
(1−φ)2
c
X
k=1
E
(ˆθ(k)syn)2−2ˆθ(k)synθ(k)+ (θ(k))2 +
2φ(1−φ)
c
X
k=1
E
(ˆθ(k)dirθˆ(k)syn−θˆdir(k)θ(k)−θ(k)θˆsyn(k)+ (θ(k))2
= φ2
c
X
k=1
(E((ˆθdir(k))2)−2φ2
c
X
k=1
θ(k)E(ˆθ(k)dir) +φ2
c
X
k=1
(θ(k))2+
(1−φ)2
c
X
k=1
E((ˆθsyn(k))2)−2(1−φ)2
c
X
k=1
θ(k)E(ˆθsyn(k)) +
(1−φ)2
c
X
k=1
(θ(k))2+ 2φ(1−φ)
c
X
k=1
E(ˆθdir(k)θˆsyn(k))−
2φ(1−φ)
c
X
k=1
θ(k)E(ˆθdir(k))−2φ(1−φ)
c
X
k=1
θ(k)E(ˆθ(k)syn) +
2φ(1−φ)
c
X
k=1
(θ(k))2
= φ2
c
X
k=1
E((ˆθ(k)dir)2)−2φ
c
X
k=1
θ(k)E(ˆθ(k)dir) +
c
X
k=1
(θ(k))2+
(1−φ)2
c
X
k=1
E((ˆθsyn(k))2)−2(1−φ)
c
X
k=1
θ(k)E(ˆθ(k)syn) +
2φ(1−φ)
c
X
k=1
E(ˆθdir(k)θˆsyn(k))
which equals
MSE(θˆcom) = E
kθˆcom−θk2
= E
c
X
k=1
φθˆ(k)dir+ (1−φ)ˆθ(k)syn−θ(k)2
=
c
X
k=1
E(φ2(ˆθ(k)dir)2+ (1−φ)2(ˆθsyn(k))2+ (θ(k))2+ 2φ(1−φ)ˆθ(k)dirθˆsyn(k)−2φθˆdirθ(k)−2(1−φ)ˆθ(k)synθ(k))
=
c
X
k=1
(E(φ2(ˆθ(k)dir)2) +E((1−φ)2(ˆθ(k)syn)2) + (θ(k))2+
E(2φ(1−φ)ˆθ(k)dirθˆ(k)syn)−E(2φθˆdir(k)θ(k))−E(2(1−φ)ˆθsyn(k)θ(k)) ) and therefore it follows that
MSE(θˆcom) = φ2MSE(θˆdir) + (1−φ)2MSE(θˆsyn) + 2φ(1−φ)E((θˆdir−θ)>(θˆsyn−θ))
For minimizing with respect toφ, the terms are rewritten in the following way:
MSE(θˆcom) = φ2MSE(θˆdir) +MSE(θˆsyn)−2φMSE(θˆsyn) +φ2MSE(θˆsyn) + 2φE((θˆdir−θ)>(θˆsyn−θ))−2φ2E((θˆdir−θ)>(θˆsyn−θ))
= φ2(MSE(θˆdir) +MSE(θˆsyn)−2E((θˆdir−θ)>(θˆsyn−θ))) + 2φ(E((θˆdir−θ)>(θˆsyn−θ))−MSE(θˆsyn)) +MSE(θˆsyn) and the first derivative to φis set to zero:
dMSE(θˆcom)
dφ = 2φ(MSE(θˆdir) +MSE(θˆsyn)−2E((θˆdir−θ)>(θˆsyn−θ))) + 2(E((θˆdir−θ)>(θˆsyn−θ))−MSE(θˆsyn)) = 0
It follows that the optimal weightφ∗is given by
φ∗= MSE(θˆsyn)−E((θˆdir−θ)>(θˆsyn−θ)) MSE(θˆdir) +MSE(θˆsyn)−2E((θˆdir−θ)>(θˆsyn−θ))
3.4 Schaible’s Approximation
The optimal weight φ∗ from Section 3.3 becomes more manageable when the term E((θˆdir−θ)>(θˆsyn−θ)) is assumed to be small relative to the mean squared error MSE(θˆsyn) of the synthetic estimator (Schaible, 1978b). Under this condition, Schaible’s approximation (see Schaible (1978b))
φ∗app= MSE(θˆsyn) MSE(θˆdir) +MSE(θˆsyn)
CHAPTER 3. PRINCIPLE OF COMPOSITE ESTIMATION 28
to the optimal weight φ can be used (see also Schaible (1979) and Schaible (1996)). Schaible’s weighting scheme can be viewed as one in which each com- ponent estimator is first weighted by the inverse of its mean squared error, and then the two component weights are normalized so that they sum to unity.
Schaible (1978a) found that the use of the approximate optimal weights rather than the optimal weights produced negligible increases in mean squared errors for selected models and variables. The approximate optimal weight φ∗app de- pends only on the ratio of the mean squared errors (Rao, 2003):
φ∗app = 1
1 +F, whereF = MSE(θˆdir) MSE(θˆsyn)
Under Schaible’s approximation, the mean squared error of the corresponding composite estimator can be approximated by
MSE(θˆcom) = (φ∗app)2MSE(θˆdir) + (1−φ∗app)2MSE(θˆsyn)
= (MSE(θˆsyn))2MSE(θˆdir)
(MSE(θˆdir) +MSE(θˆsyn))2 + (MSE(θˆdir))2MSE(θˆsyn) (MSE(θˆdir) +MSE(θˆsyn))2
= MSE(θˆsyn)MSE(θˆdir)(MSE(θˆsyn) +MSE(θˆdir)) (MSE(θˆdir) +MSE(θˆsyn))2
= MSE(θˆsyn)MSE(θˆdir) MSE(θˆdir) +MSE(θˆsyn)
= φ∗appMSE(θˆdir) = (1−φ∗app)MSE(θˆsyn)
That is, the reduction in MSE achieved by the optimal estimator relative to the smaller of the MSE of the component estimators is given byφ∗app if0≤φ∗app≤ 0.5(because this implies that MSE(θˆsyn)≤MSE(θˆdir)so that the corresponding reduction is given by MSE(θˆsyn)−(1−φ∗app)MSE(θˆsyn) =φ∗appMSE(θˆsyn)) and it equals 1−φ∗app if 0.5 ≤ φ∗app ≤ 1. Thus the maximum reduction of 50 percent is achieved whenφ∗app= 0.5(Rao, 2003). The mean squared error of the composite estimator based on Schaible’s approximation can be approximated in multiples of the mean squared error of the second component estimator, i.e. of the synthetic component estimator, in the following way:
R = MSE(θˆcom) MSE(θˆsyn)
= φ2MSE(θˆdir)
MSE(θˆsyn) +(1−φ)2MSE(θˆsyn) MSE(θˆsyn)
= φ2F+ 1−2φ+φ2
= (F+ 1)φ2−2φ+ 1
Schaible (1978b) studied the behavior of this MSE ratio as a function of φfor selected values ofF (=1,2,6). His results suggest that sizable deviations from the optimal weightφ∗appdo not produce a significant increase in the MSE of the composite estimator, that is, the curve is fairly flat in the neighborhood of the optimal weight. Figures 3.1-3.3 show the reduction in MSE and the range ofφ for which the composite estimator has a smaller MSE than either component estimator for different values ofF:
0 0.2 0.4 0.6 0.8 1 1.2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized MSE
Weight of the Direct Estimator
Composite Estimator Direct/Synthetic Estimator
Figure 3.1: Normalized mean squared errors (with respect to the mean squared error of the synthetic estimator) as a function ofφwhenF = 1 (after Schaible (1978b)).
0 0.5 1 1.5 2 2.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized MSE
Weight of the Direct Estimator
Composite Estimator Direct Estimator Synthetic Estimator
Figure 3.2: Normalized mean squared errors (with respect to the mean squared error of the synthetic estimator) as a function ofφwhenF = 2 (after Schaible (1978b)).
0 1 2 3 4 5 6 7
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized MSE
Weight of the Direct Estimator
Composite Estimator Direct Estimator Synthetic Estimator
Figure 3.3: Normalized mean squared errors (with respect to the mean squared error of the synthetic estimator) as a function ofφwhenF = 6 (after Schaible (1978b)).
Chapter 4
Direct Component Estimator
4.1 Unbiased Estimation
In contrast to small area (or indirect) estimators, direct estimators are based exclusively on data from the subject stand. In the context of this thesis only the design-unbiased approach of direct estimation is considered. Under the design- unbiased approach to sampling, randomness enters the problem only through the deliberately imposed design by which the sample plots are selected (see Thompson (2002)). Furthermore, the property
E(θˆdir) =θ
is fulfilled by making assumptions only about the design and without making as- sumptions about the target population. More concretely, the assumption about the design will be that the sample points are randomly distributed in the subject stand (based on a random variable with an uniform distribution, see Section 2.4).
Under the model-based approach, by contrast, randomness enters the problem by regarding (also) the population as random. Assuming, for instance, a process whereby each tree is dropped onto the stand areaAin such a way that the tree has the same chance of falling on any location. Under such a scenario, an unbi- ased estimate of the mean diameter of the trees could be produced by selecting the n nearest trees to a sampling point (distance sampling, see for instance Lynch and Rusydi (1999) and Lynch and Wittwer (2003)) and by building the mean value of the measured diameters of the sampled trees. However, whether or not the sample point is selected by a uniformly distributed random variable, the unbiasedness of this estimator depends on the assumptions about the pop- ulation, i.e. it is a model-unbiased estimator.
In this thesis, the design-unbiased approach has been chosen for direct esti- mation because it allows to enforce unbiasedness without making assumptions about the population but only by choosing a simple random sampling design.
The design-based approach can be considered under the classical sampling the- ory (Section 4.2) as well as under the more modern infinite population approach
30
(Section 4.3). Under the finite population approach, the population is viewed as the set of all trees. Under the infinite population approach, by contrast, the population is seen as the set of all points of the stand area together with the attribute vectors attached to these points.
4.2 Horvitz-Thompson Estimator
Under the classical sampling theory the selection of sample trees can be consid- ered as nsindependent replications of sample selection, that is at each sample point zj an entire sample of trees is collected (see Eriksson (1995)). It follows that at each sample pointzj an unbiased estimate of the number of trees per hectare in diameter classkcan be produced by the Horvitz-Thompson estimator (Horvitz and Thompson, 1952) which is given by
θˆht(k)(zj) = 1 λ(A)
ne
X
i=1
I(yi∈ I(k))ι(pi, yi,zj)
πi ,
whereπi denotes the inclusion probability P(ι(pi, yi,zj) = 1) =
Z Z
zj∈C(pi,r(yi))∩A
1
λ(A)dz(1)j dzj(2)
= λ(C(pi, r(yi))∩ A) λ(A)
of the tree with label i. In the sampling design for the subject stand, the iota functionι(pi, yi,zj)reduces toι(pi,zj)andr(yi)tor(see Section 2.4 for details). From the above formula it can be seen that the stand areaλ(A)has not to be known in order to estimate the number of stems per hectare (in contrast to the estimation of the total number of stems per hectare in the subject stand):
θˆht(k)(zj) = 1 λ(A)
ne
X
i=1
I(yi∈ I(k))ι(pi, yi,zj)λ(A) λ(C(pi, r(yi))∩ A)
=
ne
X
i=1
I(yi ∈ I(k))ι(pi, yi,zj) λ(C(pi, r(yi))∩ A)
Finally, the Horvitz-Thompson estimatorθˆht for the stand tableθis given by θˆht(zj) =
θˆht(1)(zj),θˆht(2)(zj), ...,θˆht(c)(zj)>
The only random variable in the formula of the Horvitz-Thompson estimator is the position vectorzj of the sample point, so that the expectation value of the Horvitz-Thompson estimator can be written as
E
θˆ(k)ht (zj)
= 1
λ(A)
ne
X
i=1
I(yi∈ I(k))E(ι(pi, yi,zj)) πi
, where the expectation value of the iota functionι(pi, yi,zj)is given by
E(ι(pi, yi,zj)) = Z Z
zj∈A
ι(pi, yi,zj)f(z(1)j , zj(2))dzj(1)dz(2)j
= Z Z
zj∈A
ι(pi, yi,zj) 1
λ(A)dz(1)j dzj(2)
