A Comparative Study of Spatial Analysis Methods for Forestry Nelder Experiments by David L . R . Affleck B.Sc. (Forest Science), University of British Columbia, 1998 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F Master of Science in THE FACULTY OF GRADUATE STUDIES T H E F A C U L T Y O F F O R E S T R Y Department of Forest Resources Management We accept this thesis as conforming to the required standard The University of British Columbia July 2001 © David L . R . Affleck, 2001 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstrac t The Nelder (1962) series of systematic spacing wheel designs define compact and spatially explicit layouts for forestry spacing experiments. One difficulty in the analysis of Nelder exper-iments is that the compact arrangements of trees may result in significant correlations among neighbouring variable values when these are geographically ordered. When this is the case, clas-sical analysis methods, based on random sampling models and ordinary least squares estimators, provide inefficient estimates of treatment effects and biased estimates of variance parameters. How-ever, the formal geometric structure of Nelder experiments presents an opportunity to evaluate the properties of alternative analysis methods based on spatial correlated error models. The statistical and practical properties of four spatial analysis methods were evaluated in the context of the interpretation of data from forestry Nelder (1962) experiments. The spatial analysis methods considered were based on either regionalized variables or nearest neighbour models and restricted maximum likelihood estimators. The validity and relative efficiency of the spatial analysis methods were assessed under different magnitudes of spatial autocorrelation through a simulation study. The same spatial analysis methods were also applied to data from a Nelder experiment in the University of British Columbia Malcolm Knapp Research Forest in Haney, British Columbia. The spatial analysis methods based on regionalized variables models were valid and more efficient than the classical analysis method when data were strongly spatially correlated. The spatial analysis methods based on nearest neighbour models were not valid and were less efficient than the classical analysis method. Application of the spatial analysis methods required more intensive estimation procedures. ii Contents Abstract ii Contents iii List of Tables v List of Figures vi Acknowledgements vii 1 Introduction 1 2 Nelder Experiments 3 2.1 Spacing Experiments 3 2.2 Nelder Spacing Wheel Designs 6 2.3 The M K N Experiment 9 2.4 Randomization and Statistical Inference 11 3 Analysis Methods for Nelder Experiments 15 3.1 Analysis Methods 15 3.1.1 The General Linear Model 15 3.1.2 Statistical Properties 17 3.2 Classical Analysis Methods 19 3.2.1 The A N O V A Model 19 3.2.2 OLS Estimation 20 3.2.3 Robustness 21 3.3 Spatial Analysis Methods 22 3.3.1 Spatial Correlated Error Models . . 22 iii 3.3.2 R E M L Estimation 29 3.4 Model Selection 32 4 Methods 34 4.1 Simulation Study 34 4.2 Analysis of the M K N Experiment 38 5 Results and Discussion 40 5.1 Results 40 5.1.1 Simulation Study 40 5.1.2 M K N Analysis 45 5.2 Discussion 47 5.2.1 Simulation Study 48 5.2.2 M K N Analysis 53 6 Conclusion 55 Literature Cited 57 Appendix A R E M L Estimation 63 Appendix B On the Equivalence of OLS and G L S 67 iv List of Tables 2.1 Tree density levels in the M K N experiment, by sector. . . . • 11 4.1 Simulation models, analysis methods, and parameter values used in the simulation study 35 4.2 Simulated treatments effects by tree density class 35 4.3 Analysis of variance table for the simulated Nelder experiments 38 4.4 Analysis of variance table for the M K N experiment DBH15 data set 38 5.1 Mean and standard deviations of OLS and R E M L variance parameter estimators at each level of spatial autocorrelation, by simulation trial 42 5.2 ANOVA/OLS analysis of variance table for the M K N experiment DBH15 data set. . . 45 5.3 Estimates of variance parameters, A I C , and restricted L R statistics for the M K N DBH15 data set, by analysis method 46 5.4 Rankings of estimated tree density effects for the M K N DBH15 dataset, by analysis method 47 v List of Figures 2.1 Nelder (1962) systematic spacing designs 6 2.2 A Nelder (1962) type la spacing wheel 7 2.3 Schematic of the Nelder (1962) type la spacing wheel 7 2.4 Layout of the M K N experiment 10 2.5 Aerial photograph of the M K N experiment taken in 1987 12 3.1 Rate of decay of the pairwise covariance between observations as a function of their separation distance under the exponential and Gaussian covariance functions 28 5.1 Relative variance bias and relative empirical variance of analysis methods by simu-lation trial 41 5.2 Empirical rejection rates of among- and within-class treatment F-tests by analysis method, level of spatial autocorrelation, and simulation trial 43 5.3 Empirical rejection rates of the restricted L R test, and the rates at which the AIC for ANOVA/OLS exceeded the AIC for the spatial analysis method, in each simulation trial as a function of the spatial autocorrelation 44 5.4 Frequency distribution and residual plot for the standardized ANOVA/OLS M K N DBH15 residuals 46 vi Acknowledgements I am grateful to Dr. Valerie LeMay, my supervisor, for her assistance and patience, and for many helpful comments concerning the ideas embodied in the present work and in previous drafts. I am also grateful to the other members of my committee, Dr. Antal Kozak, Dr. Peter Marshall, and Dr. David Tait, for their advice and support. Albert Nussbaum of the B . C . Ministry of Forests, Research Branch, generously provided the data for the Malcolm Knapp Research Forest Nelder experiments. Funding for this research was provided by a series of University Graduate Fellowships admin-istered by the University of British Columbia and by fellowships endowed by the Sopron Alumni and the L i Tze Fong family. A Natural Sciences and Engineering Research Council of Canada research assistantship within the Department of Forest Resources Management also provided fi-nancial support. Guillaume Therien and J.S. Thrower and Associates L td . provided me the time and support I needed during the final months of thesis preparation. Finally, I would like to thank Willow, for her support and encouragement throughout. D A V I D L . R . A F F L E C K The University of British Columbia July 2001 vii Chapter 1 Introduction Field experiments are used in forest research to evaluate the effects of different treatments or stand management regimes. Many forestry field experiments consist of large experimental units and span large areas. Spacing experiments based on the Nelder (1962) series of spacing wheel designs are notable exceptions. Individual trees are the basic units of study in these experiments, and are distributed in carefully controlled and relatively compact arrangements. The spatial dimensions of Nelder spacing wheel designs make them suitable candidates for evaluating the statistical and practical properties of spatial analysis methods. Classical analysis methods are based on the principles of random sampling and mutual independence. Data from forestry Nelder experiments may not satisfy the condition of independence because some trees grow in more similar environments than others. Given the relatively compact arrangement of trees in Nelder designs, this may result in appreciable correlations in the growth habits of neighbouring trees. Classical analysis methods cannot account for spatial autocorrelation among observations directly. When data in field experiments are spatially autocorrelated, treatment effect estimators for classical analysis methods are not efficient and variance estimators are biased (Cressie 1993). The device of experimental randomization can be relied upon to neutralize the effects of spatial autocorrelation over repeated trials (Harville 1975). However, it is generally not feasible to repeat experiments. Also, in Nelder spacing wheel designs, spatial randomization of treatments is not possible. In this thesis, spatial analysis methods based on correlated error (CE) models are exam-ined. Spatial C E models allow for spatially explicit autocorrelation among geo-referenced data. The potential advantages of spatial analysis methods are increased precision in treatment effect estimation and unbiased variance estimation (Cressie 1993). Estimates of spatial autocorrelation 1 parameters may also provide useful insights into the form and pattern of uncontrolled sources of variation. However, spatial analysis methods require more involved estimation and computation techniques. The specific objectives of this thesis were: 1. to review the effects of spatial autocorrelation among data in field experiments on the perfor-mance of classical analysis methods consisting of analysis of variance ( A N O V A ) models and ordinary least squares (OLS) estimators; 2. to examine the efficiency and validity of four spatial analysis methods consisting of spatial C E models and restricted maximum likelihood ( R E M L ) estimators, for the interpretation of forestry Nelder experiments; and 3. to evaluate the practicality of spatial analysis methods for forestry Nelder experiments. A description of the Nelder (1962) series of systematic spacing designs and a review of different perspectives on experimental inference are provided in Chapter 2. In Chapter 3, the effects of spatial autocorrelation on the OLS estimators of treatment effects and variance parameters are examined. Spatial C E models and R E M L estimators for spatial analysis methods are also introduced in that chapter. The second and third objectives of this thesis were addressed by means of the application of classical and spatial analysis methods to simulated and actual Nelder experiments, as described in Chapter 4. The results of these two applications are presented and discussed in Chapter 5. Chapter 6 is included as a summary. 2 Chapter 2 Nelder Experiments 2.1 Spacing Experiments In 1943, Toumey and Korstian (p. 18) wrote that one of the chief aims of silviculture was to bring about a stand of "crop trees of best form and most rapid growth consistent with fullness of stand and quality of product." Although somewhat dated, the preceding statement is of particular value because of its explicit recognition of the relationships among tree form, growth rate, and density (i.e., stems per hectare). Silviculturalists may strive to improve site productivity and to ensure the establishment of desired crop species, but managing the distribution of productivity among and within individual crop trees is frequently of equal importance. Often this latter objective is achieved through the regulation of tree density. Silviculturalists may regulate tree density at any age through thinnings. However, because controlling tree characteristics such as lateral branching early in stand development can greatly influence the future economic value of a crop (Nyland 1996), attending to seedling or sapling density can be a better investment. Other economic incentives and legal responsibilities may also render rapid initial tree growth and appropriate re-stocking levels expedient objectives. A detailed knowledge of the effects of initial density on the growth of forest tree crops can therefore be of considerable value. Evaluating the effects of tree density on tree growth has been the focus of many spacing experiments in forestry. Smith (1959) discussed some of the early approaches to studying the problem and Evert (1971, 1984) provided thorough reviews and bibliographies of the subject. The archetypal forestry spacing experiment consists of a number of trees that are grown at different densities by arranging them into a collection of rectangular experimental units. After a period of time, some attribute of the trees is measured (e.g., diameter or survivorship) and trees grown at 3 different densities are compared. These experiments are largely based on the same design principles as agricultural spacing experiments, but differ in two fundamental ways. First, agricultural crops often take only a sin-gle season or year to fully mature, making agricultural spacing experiments relatively short-lived. Forestry spacing experiments that persist for years or decades may become increasingly compli-cated by data compilation issues and tree mortality. Second, agricultural spacing experiments are typically designed to study the effects of both the size and rectangularity1 of the growing space on plant growth. Few spacing experiments in forestry have examined the effects of the spatial arrangement of trees independently of the effects of tree density; the study by Smith (1978) is a notable exception. It has generally been held that tree growth wil l not be unduly influenced by rectangularity, except at extremes, since tree crops have more time to adjust to the dimensions of their growing space (Namkoong 1966). At the same time, the degree of mechanization and control applied in the establishment of agricultural crops is generally not feasible for forest tree crops under operational conditions. The design of forestry spacing experiments is greatly complicated by the nature of the treat-ment employed, namely, tree density. The design of field experiments typically proceeds from the identification of discrete plots of ground as the experimental units. Appreciable systematic het-erogeneity can be found across and within experimental units, so replication and interspersion of treatments are desired. In many field experiments, it is possible to allocate treatments to exper-imental units subsequent to their demarkation. Then experimental randomization (i.e., random allocation of treatments to experimental units) can easily be accomplished. However, in forestry spacing experiments, the factor of interest is inter-tree competition, represented as a set of treat-ment levels by a range of tree densities. As such, the delineation of experimental units and the allocation of treatments must be done simultaneously. Under these circumstances, two options exist for allocating different levels of tree density to experimental units to permit random allocation or interspersion of treatments. One option is to fix the dimensions of the experimental units and plant (or space) each to the desired density. This option permits the experimenter to employ experimental randomization. However, the use of experimental units of a fixed size precludes the use of a single tree per experimental unit for all but the lowest density treatments. At higher densities, the experimental units must contain more trees. Therefore, this strategy entails a less efficient use of materials. Single-tree experimental units are economical in the use of both trees and space, and can lead to more efficient designs than 1The ratio of the dimensions of the growing space (i.e., length : width). 4 multiple-tree experimental units (Loo-Dinkins and Tauer 1987). Multiple-tree experimental units do allow for estimation of within-experimental-unit variability, but different numbers of trees per experimental unit results in differential precision of the estimates. Also, at high densities there can be far more trees than necessary for this purpose (Namkoong 1966, Smith 1978). This results in the collection of redundant information and a potentially wasteful use of trees. The alternative option for treatment arrangement in spacing experiments is to use a constant number of trees in each experimental unit. Experimental units with different tree densities will then have different areas. Unfortunately, this can result in considerable variability in the size of the experimental units, especially when multiple-tree experimental units are used or a broad range of tree densities is studied. In principle, variable-sized experimental units could be randomly allocated over the experimental field (or within blocks). In practice, this is impossible; it would not allow for close packing of plots and would result in a very inefficient use of space. Instead, the experimental field must be partitioned (systematically) into variable-sized experimental units. Often this results in the use of a single experimental unit for each density class, with sub-sampling supplying an estimate of experimental error (e.g., the study described by Reukema 1979). However, replication and interspersion of tree density classes is possible; Smith (1978) described a spacing experiment using multiple-tree, variable-sized experimental units where replicates were interspersed over the field. Ultimately, the most detracting element of spacing experiments employing experimental units with either a constant area or a constant number of trees, is the need for a large number of buffer trees. Maintaining specific tree densities within each experimental unit, while allowing for some degree of randomization (or at least interspersion) of treatments, necessarily prevents the use of spatially contiguous plots. Buffer trees are required around the perimeter of all experimental units, but do not contribute any information for the estimation of treatment effects. Moreover, the proportion of the experimental area (and trees) lost to buffer elements increases as the size of the experimental units is decreased. Thus, in one sense, any efforts to trade off experimental unit size for increased replication results in even less economical use of space and trees. In summary, either constant area or constant number of tree experimental units can be used in spacing experiments in order to adhere to some of the conventions of classic field experiment methodology, such as experimental randomization. However, neither option makes very economical use of trees or space. As a result, researchers have had to consider using alternative designs for forestry spacing experiments. Among such alternatives are Nelder (1962) systematic spacing designs. 5 2.2 Nelder Spacing Wheel Designs In 1962, Nelder published a description of a series of systematic experimental designs specif-ically for spacing trials. Two essential features of this series of designs are that the experimental units consist of single trees and that tree density, rectangularity, or both, are varied in a continuous and systematic manner across the field. Nelder described five designs: type la , in which only tree density is varied; type lb, in which only rectangularity is varied; and types Ic, Id, and II, in which these two factors are varied together but in different directions (Figure 2.1). The first two designs (types l a and lb) can be formulated so as to run through any angle, including 360°. In this case > • • > • "S. • • "s>. • . (a) (b) (c) (d) (e) Figure 2.1: Five systematic designs for spacing experiments: (a) type la : density decreases with increasing distance from the origin (rectangularity fixed); (b) type lb: rectangularity decreases with distance from the origin (density fixed); (c) type Ic: density changes along the horizontal while rectangularity changes on the vertical; (d) type Id: density changes on the vertical while rectangularity changes on the horizontal; (e) type II: density and rectangularity change along both axes but in different directions. Adapted from Nelder (1962). 6 they take on the appearance of a spoked wheel (Figure 2.2). As a result, these designs are often referred to as spacing wheels. In the following text, the term "Nelder experiment" will refer to spacing experiments consisting exclusively of one or more of these spacing wheels. Figure 2.2: A Nelder (1962) type la spacing wheel spanning 360°. Of the five designs proposed by Nelder (1962), the type la has been most widely applied in forestry research (e.g., Freeman 1964; Zedaker 1982; Cole and Newton 1987; Galinski et al. 1994). In the type la spacing wheel, trees are positioned along spokes radiating from the centre at a constant angle (vo) and concentric arcs evolve at increments (r{) along these spokes (Figure 2.3). Tree density decreases (i.e., growing space per tree increases) along the spokes while the arcs form Figure 2.3: Schematic of the Nelder (1962) type la spacing wheel showing the angle (m) between spokes and the arc radii (r,). iso-density contours. Rectangularity is constant throughout. Because trees are located at the intersections of arcs and spokes, the design is most easily described in terms of polar co-ordinates. 7 The radii (in metres) of the concentric arcs form a geometric progression given by U = r j _ i x k = r 0 x kl (2.1) where rj is the radius of the i th arc ( i = 1, 2, . . . , £ , (£ + 1)), ro is the radius of the innermost arc, k is a rate constant, and (t + 1) is the index of the outermost arc. The innermost and outermost arcs of a spacing wheel act as buffers so there are only t experimental levels of tree density. If the design spans less than 360°, then trees on two edge spokes must also stand as guards. The area per tree (Af, in m 2 ) and the tree density (Di\ in stems/ha) anywhere on the i th arc in a type la spacing wheel are given by Ai = rf zu (k - k~1)/2 = r% k2i zu (k - k~l)/2 (2.2) and Di = lOOOO(^)" 1 (2-3) for i = 1,2,... ,t. The rectangularity (u) anywhere in a type l a is v = w/ (k* - k~*y . (2.4) The type l a spacing wheel can be completely described by specifying t, ro, A;, and zu. Alternatively, the latter three parameters can be derived given t, v, and the range of tree densities of interest. Denoting the areas per tree at the endpoints of the density range by A\ and At, the system of equations: ' logfc = l o g ( £ ) / ( 2 - 2 t ) < w = v(kv-k-?) (2.5) r 0 = v/2Ai /o7(* : 3 -* : ) can be solved to yield k, zu, and ro. In his original paper, Nelder (1962) noted that as the rate constant k is increased, the deviation in plant position from the centroid of the quadrangle formed by its four nearest neighbours can become substantial. The distances to the two nearest trees on neighbouring spokes are the same, but the distances to the two nearest trees on the same spoke are not. Nelder defined the non-centrality (C 0 ) of a type la spacing wheel as the deviation in the position of a plant from the centroid of its growing space, taken a percentage of the length of the growing space along the spoke. This can be expressed as 8 Nelder advised that C 0 should be maintained below 5% (i.e., k < 1.1) in agricultural spacing experiments to prevent bias in the estimation of treatment effects. However, Namkoong (1966) pointed out that, just as with rectangularity, because trees will have more time to adapt to the shape of their growing spaces, the effects of non-centrality are likely to be less significant in forestry spacing experiments. Many forestry spacing experiments employing the type l a design have used k > 1.1 (e.g., Zedaker 1982; Paltzat 1984; Cole and Newton 1987; Galinski et al. 1994). Several variations on the original Nelder (1962) spacing wheels have been proposed. In 1966, Namkoong described a number of designs based on the Nelder type l a spacing wheel, in which rectangularity was varied over a relatively narrow range (approximately 1:1 to 1:4). The rationale for these modified designs was that the range of tree densities to be studied in the type la spacing wheel is spanned by a geometric series and most of the density levels fall towards the upper end of the range. Namkoong recognized that spacing wheels with constant density increments could be formulated if rectangularity were varied. More recently, Rustagi (1984) modified the Nelder design to allow for square or rectangular tree spacing. The continuous change in plant density, characteristic of Nelder designs, also forms the basis of several other systematic spacing designs used in intercropping research such as the replacement series design (Radosevich 1987). This latter design closely resembles the rectangular type II design (Figure 2.1e). The type la spacing wheel is an extremely efficient design for studying tree density-growth relationships. When a spacing wheel spans 360° degrees, only the inner and outer arcs (i.e., a total of ^£ trees) are lost as guards. Therefore, greater proportions of area and trees are allowed to contribute information to the estimation of treatment effects. Since the design is so compact, it is often possible to replicate it several times within the same area required by an experiment with multiple-tree, rectangular experimental units. The spacing wheel design is also very flexible; designs running through angles less than 360° can be used in various combinations to conform to the shape of the experimental field (e.g., Freeman 1964; Namkoong 1966). 2.3 The M K N Experiment In 1977, the staff of the University of British Columbia Malcolm Knapp Research Forest in Maple Ridge, British Columbia, initiated a Nelder experiment to investigate the effects of spacing in red alder (Alnus rubra Bong.) plantations ( M K N Experiment). At the time, the prevalence of red alder regeneration on cutovers was sufficient for some to consider managing the species on short rotations (Paltzat 1984). However, up to that point, the species had been widely regarded 9 as a weed species and very little was known about the effects of initial density on its growth and development. The M K N experiment was established at an elevation of 200 m, on a site logged the previous year (Paltzat 1984). Red alder seedlings were planted following scarification. The experiment is made up of two 0.35 ha spacing wheels, as shown in Figure 2.4. Both spacing wheels consist of Figure 2.4: Layout of the M K N experiment. Each of the two spacing wheels consist of three sectors (A, B, and C) with different design parameters. Buffer trees are shown in grey. three sectors, and are oriented so that each wheel is the mirror image of the other. The three sectors are defined by different design parameters and span different tree density levels (Table 2.1). Sectors A and B in Figure 2.4 are both type l a spacing wheels (Nelder 1962). Sector A was designed with ro = 3.93 m, w = 19.59°, and k = 1.3797, while sector B was laid out using ro = 3.93 m, w = 9.795°, and k = 1.1746 (Paltzat 1984). Sector C is a modified type l a design based on a set of design parameters given by Namkoong (1966). Trees on the innermost and outermost arcs, and on spokes that frame different sectors, are used to buffer the experimental trees and maintain the densities given in Table 2.1. These buffer trees are shown in red in Figure 2.4. In all there are 19 spokes per spacing wheel, 6 of which are for buffer trees, and a total of 113 experimental (i.e., non-buffer) trees grown at 25 tree densities per spacing wheel. Data were first collected from the M K N experiment in 1983 when the trees were 7 years old. 10 Table 2.1: Tree density levels in the M K N experi-ment, by sector (see Figure 2.4). Experimental Tree density (stems/ha) Arc Sector A Sector B Sector C 1 3,042 17,001 3,089 2 1,597 12,321 2,348 3 850 8,931 1,607 4 450 6,475 875 5 242 4,693 267 6 131 3,401 134 7 8 2,464 1,787 9 1,295 10 949 11 690 12 504 13 368 The trees have been re-measured twice since then, first at age 15, and then at age 19. A n aerial photograph taken in 1987, when the trees were 11 years old, is shown in Figure 2.5. Mortality has been an issue in this experiment. Paltzat (1984) reported that during the first year, 11% of the trees in one of the spacing wheels, and 61% of the trees in the other, died. The dead seedlings were culled and their positions were replanted in the second year. 2.4 Randomization and Statistical Inference The primary objective of any spacing experiment is to obtain the most accurate estimates of the effects of tree density possible, given the resources at hand. Nelder spacing wheels are very economical designs, but have not found widespread use due, in large part, to historical biases for experimental randomization and design-based statistical inference. The systematic arrangement of treatments makes these designs appealing from a practical standpoint but also renders them invalid as a basis for statistical inference. Since the Nelder spacing wheel designs do not allow for randomization of treatments to experimental units, these designs do not constitute a basis for assessing the bias or precision of an estimator of the effects of tree density. In any field experiment employing a systematic design, it is possible that a set of experimental units assigned the same treatment may also be exposed to a common level of another uncontrolled factor. As an extreme example, consider a spacing wheel 11 Figure 2.5: A e r i a l photograph of the M K N experiment taken in 1987, when the trees were 11 years o l d (scale is approximately 1:1,250). la id out such that the centre of the wheel coincides w i t h a part icular ly elevated microsite and the arcs of the wheel follow successive contours in elevation. In this case, the effects of tree density are completely confounded w i t h the effects of any site factors that vary w i t h the microtopography (e.g., soil moisture). F r o m a design-based approach to statist ical inference, confounding of treatment and uncontrolled factor effects is a major flaw of systematic experimental designs. T h e fact that it is unlikely that the concentric patterns of tree density in a spacing wheel are going to overlap w i t h s imilar expressions of site variables is not important . Nor is it important that part icular layouts of randomized experimental designs may lead to a confounding of treatment and other effects; the resultant bias is zero when averaged over the set of a l l possible randomizations. A second cr i t i sm of the use systematic designs i n field experiments stems from the properties of the experimental material . It has often been observed (e.g., Student 1914; Papadakis 1937; M i l n e 1959; Magnussen 1993) that data from field experiments tend to exhibit systematic spat ia l hetero-geneity such that there can be appreciable correlation among neighbouring observations. T h i s poses a problem for statist ical inference. Spat ial autocorrelation violates the conditions under lying many of the common normal-theory variance estimators and hypothesis tests. T h e device of experimental randomizat ion guarantees the existence of a design-unbiased estimate of experimental error and the 12 validity of inferences on treatment effects. In randomized designs, randomization of treatments to experimental units neutralizes the effects of spatial dependence because it ensures that all spatial arrangements of treatments are equally likely and will occur with the same frequency over the set of all possible randomizations. From the design-based approach to statistical inference, Nelder experiments do not secure a basis for inference. Yet while there is no basis for statistical inference imparted by the design, it is still possible to draw inferences on the effects of tree density in a Nelder experiment on the basis of a statistical model. Model-based statistical inference relies on an explicit definition of a stochastic framework under which the data arise. Inferences are made on the basis of the model and are subject to the adequacy of the model. Under a model-based approach, randomization renders inference more robust, but is not necessary for valid statistical inference (Harville 1975; Thornett 1982; Hooper 1989). Inference can be based on the statistical model alone, the assumptions of which can be examined independently. For example, Galinski et al. (1994) studied tree height, differentiation in a Nelder experiment and investigated the tenability of their model through a series of antecedent soil and tree sampling efforts. It is also noteworthy that even in systematic designs, some sources of variation may be controlled through randomization. Chacko (1965) pointed out that variation due to differences in tree growth potential (i.e., genetic differences) in spacing experiments is often of greater magnitude than variation due to differences in site factors. The effects of this source of variation can be controlled in Nelder spacing wheels through the randomization of seedlings to tree density levels. Much of the controversy over the use of Nelder spacing wheels and other systematic designs in field experiments stems from a misunderstanding of the role of experimental randomization in statistical inference. The literature on the design and analysis of experiments extends back to the early 1900's, yet few authors have articulated, or acknowledged the implications of, their particular preference for design- or model-based inference.2 In the context of design-based inference, the qualifier "over the set of all possible randomizations" has not always been adequately stressed, nor have all its implications been adequately addressed. As Harville (1975) pointed out, experimental randomization serves as a useful basis for statistical inference only when we ignore the fact that any analysis is conditional upon the particular randomization observed. We must also ignore the fact that the properties of the estimators derived from viewing the observed layout as one amongst an infinite set of layouts do not apply to the estimates based on the single layout observed. The design-based approach to statistical inference is associated more with an emphasis on significance 2These topics have, however, received more attention in the sampling literature (e.g., Hansen et al. 1983). 13 testing, and, as Student (1937) pointed out, does not generally lead to the most accurate estimates of treatment effects. Model-based inference can be a pragmatic perspective for the analysis of field experiments, and systematic experimental designs, such as Nelder spacing wheels can provide tremendous material efficiencies. Combined with a suitable analysis method, these two approaches to experimental inference may provide more accurate estimates of treatment effects. 14 Chapter 3 Analysis Methods for Nelder Experiments 3.1 Analysis Methods The interpretation of experimental data from a model-based approach to statistical inference proceeds from the identification and application of a suitable analysis method. A n analysis method is here defined as consisting of a parametric stochastic model and a set of estimators for the model parameters. In this section, a general linear model for Nelder experiments is developed and the statistical properties of an analysis method are defined. 3.1.1 The General Linear Model Consider a Nelder experiment consisting of b spacing wheels, t levels of tree density, and m trees per spacing wheel per density level. The observations on the variable of interest can then be modelled as a realization of the n x 1 random vector where n — htm, is the total number of observations. In addition, the relative spatial positions of the trees, denoted by the n x 2 matrix S = { s i , S 2 , . . . ,sn}T of bivariate spatial co-ordinates (i.e., Si G S 2 ; i = 1,2,...,n), can be derived from the formal geometric structures of the Nelder spacing wheels. Thus, the vector of observations, indexed by their relative spatial positions can be given 15 by the realization ys = {Vsi, 2/s2) • • • > Vsn}T of the random vector Yg. The general model developed below is based on the assumption that the elements of Ys can be additively decomposed as b t b t YSi = fi + Yl /(M>* + E /(*» + E E K)(A x T ) J * + £ S l (3.i) j=l fc=l j=l k=l where p is the overall mean, ctj is the effect of the jth spacing wheel, Tfc is the effect of the fcth level of tree density, (a x r)jk is the interaction effect for the jth. spacing wheel and the kth level of tree density, eSi is an unobservable zero-mean random variable (indexed by its spatial position) associated with the observation taken on the ith tree, and { 1 if the i th observation is assigned to the qth factor level 0 otherwise. Removing the linear dependencies among the mean structure parameters, Equation 3.1 can be stated more concisely as YSi = Xi/3 + eSi (3.2) where X j is a 1 x p design vector relating observations to the main effects and their interactions, (3 is a p x 1 vector of linearly independent parameters (i.e., p — bt), and eSi is as defined above. The tree density effects in Equation 3.1 are then obtained through a linear transformation of the elements of f3, that is = A?> (3.3) where Xt is a p x t matrix of ones and zeros that defines the appropriate transformation of (3. Estimation of r (treatment effects) and of the differences among the elements of r (treatment contrasts) are of central importance. In Equations 3.1 and 3.2, tree density is implicitly defined as a class variable. Researchers studying the effects of tree spacing are generally interested in modelling the effects of tree density over a range of tree densities, not only at the t levels used in an experiment. However, the effects of tree density are likely to be different for different variables of interest, tree species, and site qualities. Since the emphasis in this study is on comparing different estimators of tree density effects, the discontinuous form in Equation 3.2 is used. This avoids complicating these comparisons with possible effects of mis-specifying the density-growth response. Notwithstanding, the results of this study could be extended to the estimation of more complex density-growth 16 T = n Tt response functions. Continuous and/or non-linear functions could be fitted by generalizing Xj/3 to a more flexible mean structure (e.g., / (x, ; /3)) . Spacing wheel effects and interactions among spacing wheels and tree densities can also be obtained from a transformation of j3. These factors could be modelled as either random or fixed effects. However, the focus of this study is the estimation of the effects of individual levels of tree density (or of contrasts among them) not the estimation of the expected growth response. As a result, variation due to differences among spacing wheels, and variation due to interactions among spacing wheels and tree density levels, can be treated as fixed effects without loss of generality (Casella and Berger 1990). Throughout this thesis, it was assumed that the joint probability density function (pdf) of the random errors (eS i) in Equation 3.2 could be approximated by a multivariate normal (i.e., Gaussian) distribution. The Central Limit Theorem (CLT) states that the sum of a sequence of independent random variables converges in distribution to a normal random variable (Chung 1974).1 The random errors in Equation 3.2 consist of the effects of many sources of variation (e.g., measurement error, differences among trees used in the experiment, variability in soil properties within spacing wheels, etc.). If the effects of these sources of variation can be taken to be additive, the C L T can be relied upon to provide a measure of assurance that their cumulative effects will be normally distributed. A multivariate normal distribution may not be suitable when the tree variable of interest is discrete (e.g., tree vitality) or has a censored domain (e.g., diameter increment cannot be negative) but may be a reasonable approximation for many continuous tree growth attributes: At this point, Equation 3.2 can be restated in matrix notation as Y s = X /3 + es (3.4) where X = {xi , X 2 , . . . , x n } T is a n x p matrix of full column rank, and es = {eSl,sS2,..., e S n } T is a re x 1 vector of random errors. In addition, es ~ N (0, V ) where 0 is an n x 1 vector of zeros, and V is an n x n positive definite matrix. Finally, let V — t r 2 f l where o2 is a general scale parameter (cr2 > 0) , and Q is an n x n positive definite matrix. This implies that Y s ~ N (X/3, a2£l). 3.1.2 Statistical Properties Each analysis method discussed in the following sections is based on the general linear model given by Equation 3.4. Evaluation of these analysis methods will be made on the basis of both practical and statistical considerations. Practical considerations including parsimony, biological 1 Chung also details several formulations of the CLT in which the condition of independence is relaxed. 17 relevance, and feasibility will be discussed in Chapter 5. Statistical considerations will be discussed here. The statistical properties of an analysis method considered here are bias, efficiency, and validity. A n analysis method is unbiased if the estimators of treatment effects are unbiased, that is E ( f ) = E ( A t r ) 9 ) = T (3.5) where T = A;T/3 is an estimator of r , /3 is an estimator of f3, Xt is the appropriate linear transfor-mation of /3, and the expectation is taken over the pdf defined by the model. A n analysis method is efficient if the estimators of treatment effects are efficient, that is Var ( f ) = Var ( r 0 ) - A (3.6) where f0 is any other unbiased estimator of r and A is a t x t positive semidefinite matrix. In other words, the estimators of treatment effects should be as precise as possible. Finally, if an analysis method is unbiased and E .k=i = £ V a r ( f f c ) (3.7) k-i where Var (T^) is an estimator of Var (r^) (A; = 1 ,2 , . . . , t), then it is valid (Besag and Kempton 1986). Validity is an important property of an analysis method because, in practice, analysis meth-ods are often selected on the basis of the estimated precision of T. If these estimates are biased, the comparison is meaningless. Lack of validity is also dangerous, because it can vitiate the results of hypothesis tests. For example, a testable hypothesis in the context of the general linear model developed above is the absence of a treatment effect (i.e., Ho • T = 0). Assuming all other factors in the model (e.g., spacing wheels) are considered fixed, this hypothesis can be tested using an F-statistic of the form F«r*H»-rt ~ M W ( 3 - 8 ) where MSTRT is the mean squared treatment effect; M S ^ ; ^ is the mean squared error; and (t — 1) and (n — p) are the treatment and error degrees of freedom, respectively. In matrix notation, Equation 3.8 can be stated as tTLT (hXf ( x ^ V - i x ) XtLT) _ 1 L T ^(t-i),(n-p) = ^ — (3-9) 18 where V is an estimator of V , and L is a (t — 1) x t contrast matrix L = 1 - 1 0 1 0 0 0 0 0 0 0 0 0 0 (3.10) 0 0 0 0 ••• 1 - 1 (Beckers 1997). In a randomization context, validity is defined as the equality of the expectations of the treatment and error mean squares in Equation 3.8 under Tio, which implies Equation 3.7 (Azais et al. 1998). In the context of the general linear model given by Equation 3.4, the equality of the expectations of the treatment and error mean squares also implies that the F-statistic in Equation 3.9 will be distributed according to a Fisher (F) distribution with t—1 and n — p degrees of freedom. If an analysis method is not valid, it is not legitimate to test Ho using an F-statistic. 3.2 Classical Analysis Methods The classical analysis method, A N O V A / O L S , 2 for field experiments is characterized by the assumption that the errors in Equation 3.4 are normal, independent, and identically distributed (i.i.d.) random variables. Estimation of model parameters proceeds using the ordinary least squares (OLS) estimators. If the assumption of normal i . i .d . random errors is justified, the OLS estimators of treatment effects are unbiased and efficient, and unbiased estimates of their variances are readily available. 3 .2 .1 T h e A N O V A M o d e l The analysis model for ANOVA/OLS is a special case of Equation 3.4. It can be stated as Y = X/3 + e (3.11) where e ~ AT (0, a2 I „ ) . The A N O V A model does not incorporate the information in S on the spatial arrangement of observations, so the subscript is dropped from Y and e. Also, the matrix Cl is reduced to I n , a n x rc identity matrix. This implies that a2 i f i = Z CoV(Yi,Yl) = Cov(ei,el) = for i, I = 1 , 2 , . . . , rc, and that Y ~ AT (X/3, a2ln). 0 otherwise 2The naming convention for analysis methods is [analysis model\/[estimation procedure]. 19 3.2.2 OLS Estimation The simple form of the variance matrix of the random errors, £, assumed in the A N O V A model establishes a set of conditions under which the OLS estimators of r are unbiased and efficient. The OLS estimator of r is given by TO1S = \? ( X T X ) _ 1 X T Y . (3.12) Equation 3.12 defines a linear combination otY~N (X/3, a2In), so it follows that T0IS also has a multivariate normal distribution. Its mean and variance matrix are given by E ( t o Z s ) = \[ ( X T X ) _ 1 X r E ( Y ) = AjT/3 = T (3.13) and Var (r0is) = E ( f o i s - r ) (T0IS - r)1 = E ( X T X ) 1 X r e e r X (X T X) \ = o2\[(xTx) X A t . (3.14) Equation 3.13 states that T0IS is an unbiased estimator of r . It can also be shown that Var(r 0 / S ) achieves the Cramer-Rao lower bound (Judge et al. 1985). Therefore, T0IS is an efficient estimator of r . In addition, the mean squared error 2 _ ( y - x ^ ) r ( y - x ^ s ) aois — = (3.15) n—p n—p where f3ols is the OLS estimator of (3 and e 0z s is the n x 1 vector of ANOVA/OLS residuals, provides an unbiased estimator of a2 (Judge et al. 1985). Substituting a2ls for a2 in Equation 3.14 Var ( f ols) = a2olsXf ( x T x ) _ 1 A* (3.16) provides an unbiased estimator of Var ( r 0 / S ) . Under the assumption of multivariate normal i . i .d . random errors, ANOVA/OLS is unbiased, efficient, and valid. Under the null hypothesis (i.e., r = 0), the treatment F-statistic T J , L t (LAf ( X T X ) _ 1 A , ! / ) " ' L f „„ F < ' - 1 « " - " = 5 «L(«-0 ( 3 ' 1 7 ) is distributed according to an F distribution with t — 1 and n — p degrees of freedom. These are very powerful results. 20 3.2.3 Robustness One drawback of the classical analysis method is that the OLS estimators lose their well-defined and attractive properties when the condition of normal i . i .d . errors cannot be justified. As discussed in Section 2.4, data from field experiments can exhibit spatial dependence, such that appreciable correlation can exist among observations located close together. In such cases, the form for V taken in Equation 3.11 is oversimplified (i.e., V = a2ft ^ cr2In). When the random errors are not independent, TQIS is normally distributed and unbiased but has variance Var (Tots) = E (tols - r ) ( r o / s - r ) T = E Xj (xTx) ^ e ^ X ^ x ) 1 At a2x[ ( x r x ) _ 1 xTn x ( x T x ) _ 1 At (3.18) Using Equation 3.18 it can be shown that the OLS estimators are not efficient under spatial autocorrelation. The generalised least squares (GLS) estimator of r given by - l fgu = Xf ( x r V " 1 x ) X T V - X Y (3.19) is normally distributed (under the assumption of multivariate normality for e ) with mean and variance E (fgU) = XT ( X T V - 2 X ) _ 1 X T V - X E ( Y ) = T (3.20) and Var (Tgis) = E (tgis - T) (fgis - r E xf ( x ^ - ^ x ) " 1 x T v - 1 £ £ r v - 1 x ( x r v - 1 x ) " 1 A4 = Xf (X^-'X)'1 XTV-'E ( e e r ) V " X X ( x ^ " ^ ) " 1 A = xj ( x ^ - ' x y 1 Xt = a2Xf ( X r n - 1 x ) _ 1 At . (3.21) From Equation 3.20, fgis is an unbiased estimator of r . More important, however, is that V a r ( r 0 / S ) — V a r ( f s ^ ) = A is a positive semidefmite matrix (Judge et al. 1985). As a result, when the error terms have a non-diagonal variance matrix, the OLS estimator of r is, in general, not efficient. 21 -E -E - E In addition, when ^ In the OLS estimator of a2 in Equation 3.15 is biased \r(eT ( l n - X ( x T x ) _ 1 X T ) e~ tr ( - x ( X T X ) _ 1 x T ) £ £ l " ( i n - x ( x r x ) _ 1 x r ) v n — p 1 n — p 1 n — p 1 n — p 1 -tr n — p ,2 = ^ ^ t r ( l n - X (X T X) 1 XT^) 11 (3.22) where tr(-) is the trace operator. The resulting magnitude and sign of the bias in a2ls wil l depend on the structure of the design and correlation matrices (Judge et al. 1985). But as a result, the treatment F-statistic in Equation 3.17 wil l also be biased. In summary, when the multivariate normal i . i .d. condition of the A N O V A model is justified, ANOVA/OLS is unbiased, efficient, and valid. However, if the errors are autocorrelated, the OLS estimators of r are not efficient and the variance estimators are biased. ANOVA/OLS will still be unbiased but, in general, will be neither efficient nor valid. 3 . 3 S p a t i a l A n a l y s i s M e t h o d s The gains in efficiency of treatment estimation provided by G L S estimators are of limited value since V is usually unknown. In addition, because there are more elements i n V than there are observations, it is impossible to estimate V directly. The basis of the spatial approach is to reduce the dimensions of V to an estimable number of covariance parameters by making assumptions about the correlation structure (i.e., f i ) . The resultant models are generally nonlinear with respect to the covariance parameters and require likelihood-based estimation. 3.3.1 Spatial Correlated Error Models Spatial C E models are derived from the principle known as Tobler's Law (e.g., Antle and Marshall 1996), that "everything is related to everything else, but near things are more related than distant things" (quoted in Gould 1970). The spatial C E models considered in this thesis are based on the parameterization V = <r2fi where fl = fi(S;0) is a matrix-valued function of a vector of unknown parameters 6 and the matrix of spatial coordinates S. In these models, the off-22 d i a g o n a l e l e m e n t s o f Cl (i.e., the set of p a i r w i s e c o r r e l a t i o n s a m o n g o b s e r v a t i o n s ) are m o d e l l e d as a f u n c t i o n of t h e s p a t i a l d i s t r i b u t i o n o f t h e d a t a . A s i n m a n y t i m e - s e r i e s m o d e l s , t h e a u t o r e g r e s s i v e f o r m u l a t i o n of a u t o c o r r e l a t i o n p l a y s a c e n t r a l role . H o w e v e r , d i s t i n c t a u t o r e g r e s s i v e f o r m u l a t i o n s c a n b e u s e d i n a s p a t i a l s e t t i n g a n d s p a t i a l a u t o c o r r e l a t i o n c a n b e m o d e l l e d as e i t h e r a c o n t i n u o u s o r d i s c r e t e p r o c e s s . I n a t e m p o r a l s e t t i n g , t h e a u t o r e g r e s s i v e m o d e l o f a u t o c o r r e l a t i o n c a n b e s p e c i f i e d i n a n y o f t h r e e e q u i v a l e n t w a y s . D e n o t i n g t h e o b s e r v a t i o n t o b e t a k e n at t i m e q (q = . . . , — 1 , 0 , 1 , . . . ) b y Yq, t h e f i r s t - o r d e r a u t o r e g r e s s i v e m o d e l c a n b e s p e c i f i e d : 1. b y d e f i n i n g Yq as a f u n c t i o n of Yq-\ Yq = 9Yq-! + eq (3.23) w h e r e E(Yq) = 0 a n d eq ~ N ( 0 , a2) is a r a n d o m e r r o r u n c o r r e l a t e d w i t h Yq^a a n d £ 9 _ a for a > 0; 2. t h r o u g h t h e c o n d i t i o n a l d i s t r i b u t i o n of Yq g i v e n Y * = {Yi : i ^ q) E(Yq | Yq*) = (3.24) V a r ( y , | Y*) = a2 (3.25) w h e r e yq-\ is a r e a l i z a t i o n o f Y g _ i ; o r 3. t h r o u g h t h e u n c o n d i t i o n a l d i s t r i b u t i o n o f Yq E(Yq) = 0 • (3.26) C o v ( y „ y , _ a ) = r ^ 2 ^ a i (3.27) ( C l i f f a n d O r d 1981). T h e s e d e f i n i t i o n s are r e f e r r e d t o as t h e s i m u l t a n e o u s , c o n d i t i o n a l , a n d u n c o n -d i t i o n a l f o r m u l a t i o n s , r e s p e c t i v e l y . A s C r e s s i e (1993) p o i n t e d o u t , t h e e q u i v a l e n c e of these t h r e e f o r m u l a t i o n s o f t h e t e m p o r a l a u t o r e g r e s s i v e m o d e l relies o n t h e u n i d i r e c t i o n a l n a t u r e o f d e p e n d e n c e t h r o u g h t i m e . I n a s p a t i a l s e t t i n g t h e r e is n o p a s t , p r e s e n t , o r f u t u r e - s p a t i a l d e p e n d e n c e e x t e n d s i n a l l d i r e c t i o n s s i m u l t a n e o u s l y . A s a resul t , s p a t i a l C E m o d e l s t h a t d e f i n e s p a t i a l a u t o c o r r e l a t i o n a c c o r d i n g to s i m u l t a n e o u s , c o n d i t i o n a l , o r u n c o n d i t i o n a l f o r m u l a t i o n s are d i s t i n c t . A s e c o n d d i s t i n c t i o n a m o n g s p a t i a l C E m o d e l s c a n b e m a d e o n t h e b a s i s o f t h e s p a t i a l e x t e n t of t h e a u t o c o r r e l a t i o n . T h e c o l l e c t i o n of responses , y s , f r o m a N e l d e r e x p e r i m e n t c a n b e c o n s i d e r e d as a r e a l i z a t i o n of a s p a t i a l s t o c h a s t i c p r o c e s s Y s = {YSi =Y(Si): Si EDcn2,i = 1,2,..., n} (3.28) 23 where s, ranges over a 2-dimensional index set D. A dichotomy between spatial C E models can be drawn on the basis of whether the index set, D, is discrete or continuous. One class of spatial C E models that have been widely applied in the agricultural and social sciences are the nearest neighbour (NN) models. In N N models, the n rows of S (the locations of the n trees) define the extent of the stochastic process Y s (i.e., D = S) . Under this frame-work, spatial autocorrelation is defined using spatial analogues of the simultaneous and conditional formulations of temporal autoregression. Consequently, analyses based on N N models account for spatial autocorrelation through differencing or using responses from neighbouring observations as covariates. A second class of spatial C E models applicable to the analysis of field experiments originated from geological spatial interpolation methods developed in the 1940s (Cressie 1990). Regionalized variables (RV) models have been widely used in geological and mining applications but were only introduced as potential models for the analysis of field experiments by Zimmerman and Harville in 1991. These models characterize spatial autocorrelation as a continuous process (i.e., D C TZ2) and incorporate spatial autocorrelation through explicit specification of the covariance structure. This approach is analoguous to the unconditional formulation of the temporal autoregressive model. Before examining N N and R V models in detail, it is important to note some of the assump-tions involved in modelling spatial autocorrelation. In general, autocorrelation among spatial data poses a major problem for statistical inference. Unless assumptions are made about the joint pdf of the spatial process, y s can only be viewed as a single (and possibly incomplete) realization of the stochastic process Y s - To avoid this interpretation, Y s is generally modelled as spatially stationary stochastic process. A stochastic process is spatially stationary if the joint pdf of Y s depends on S only through the relative spatial arrangement of the observations (Wackernagel 1998). Under the assumption of multivariate normality, this is equivalent to the conditions of covariance (or second order) stationarity E ( Y S ) = X /3 (3.29) and Cov(YSi,YSl) = *2g(hu\0) (3.30) where the parameters /3 do not depend on S, and g(-) is a function of 8 andhu — (s^ — s>), for all Si and S/ in D (i, I = 1 ,2 , . . . , n). Under the condition of stationarity, the pdf of a spatial stochastic process depends only on the relative spatial distribution of the observations, not on their absolute positions. In the absence of fixed large-scale spatial trends, covariance stationarity is generally a 24 reasonable assumption. A second condition that is often imposed on the pdf of spatially autocorrelated processes is isotropy. A stationary process Y s is isotropic if its joint pdf depends on S only through rf«=||hii 11= hu (3-31) for all Sj and S; in D. Under covariance stationarity, this implies that CoV(YSi,YSl)=a2f(du\0) (3.32) for all Si and Si in D. Hence the correlation between any two observations depends on their spatial positions only through their separation distance and is direction invariant. Isotropy is generally adopted as a simplifying assumption and is not a necessary condition for statistical inference. Anisotropy (i.e., pairwise correlations are a function of the directional distances between obser-vations) may be more suitable when, for example, the observations are located on a topographic gradient. In such cases, there may be much stronger pairwise correlations extending downslope than across the slope. N e a r e s t N e i g h b o u r M o d e l s In N N models, spatial autocorrelation is accounted for indirectly through simultaneous or conditional autoregressive formulations. Central to the use of N N models are the definitions of neighbourhoods and spatial weight matrices. The simultaneous spatial autoregressive (SAR) model was introduced by Whittle (1954). The error correlation structure in the S A R model is given by e s = 0 W e s + T) = (In-dWy'r, (3.33) where 6 is an unknown spatial autocorrelation parameter, 77 ~ A^(0, <r2In) is a n x 1 random vector, and W = {wu} is a n x n spatial weight matrix with elements { l/dii if es, is in the neighbourhood of es- (3.34) 0 otherwise. For a Nelder experiment, it is useful to define the neighbourhood of a tree as the set of nearest inter-spoke and intra-spoke trees. The spatial weight matrix in Equation 3.34 then reflects the 25 assymetrical distribution of these neighbours by weighting their influence as a diminishing function of their distance from tree i. Combining Equations 3.4 and 3.33, the S A R model can be stated as Y s = X / 3 + (In - eW)-1^ . • (3.35) The mean and variance of Y s are E ( Y S ) = X / 3 + ( I n - ^ W ) - 1 E ( r 7 ) = X / 3 (3.36) and V a r ( Y s ) = E ( Y s - X/3) ( Y s - Xpf = E(i n - 0 W ) " W ( i n - ew7')-1 = ff2(In-0W)-1(In-0Wr)-1 = a2(In-29W + 62WW)-1 . (3.37) Under the assumption of multivariate normality, the S A R model implies that Y s ~ A r ( x / 3 , C T 2 ( I n - 2 0 W + 0 2 W W ) - 1 ) . (3.38) A disadvantage of the S A R model is that the autoregressive process, Y s , is not independent of the underlying i . i .d . random process, rj. That is C o v ( Y s , ri) = E ( ( I n - 0 W ) ~ V ) = <r2(In - 0 W ) - 1 (3.39) is not a diagonal matrix. Consequently, application of OLS or iterative least-squares methods (e.g., Mora 1996) will not result in consistent estimators of 8 (Anselin 1988). In addition, this result has prompted some authors (e.g., Cressie 1993) to suggest that the conditional spatial autoregressive (CAR) model (Besag 1974) is a more suitable model for spatial autocorrelation. The C A R model can be specified through the conditional means and variances of the ele-ments of Y s E ( F S i I Y * J = Xi/3 + 9J2 wu(yai - xi/3) (3.40) and Var(Y S i | Y * . ) = a2 (3.41) where Y * . = {YSl : I ^ i}, a2 is the conditional variance of YSi, W = {u>ij} is a spatial weight matrix, and 8 is an unknown autocorrelation parameter. Besag (1974) showed that under the assumption of multivariate normality the C A R model can be expressed unconditionally as Ys~N(X/3,o-2(In-8W)-1) . (3.42) 26 The spatial weight matrix used for the S A R model (Equation 3.34) can also be used for the C A R model. For Equations 3.38 and 3.42 to define valid pdfs, the variance matrices V = ff2(In-2W + ^ 2 W W ) - 1 (3.43) and V = a2(In-6W)-1 (3.44) must be positive definite. For the C A R model this means that 0 < 6 < l/u)(\) where u>^ is the largest eigenvalue of W (Wall 2000). For the S A R model, this requirement can be met by the less restrictive condition 9 ^ l/^>(i) where ui^ is the i th largest eigenvalue of W (i = 1,2,... ,n). However, in order to maintain the interpretation of 9 as a spatial autocorrelation parameter, the restriction 0 < 9 < l/w(i) should also be applied in the S A R model (Wall 2000). A drawback of the N N class of models is that spatial autocorrelation is modelled indirectly. Examination of Equations 3.38 and 3.42 reveals that the autocorrelation structure is modelled through the inverse of the covariance matrix. As Wall (2000) pointed out, this means that the actual spatial autocorrelation structure cannot be understood simply by examining W . In addition, except on infinite regular lattices, S A R and C A R models are not stationary, and the diagonal elements of Q. = ( I n - 20W + ^ W W ) ' 1 and ft = (I„ - 9W)'1 are not all equal to 1. Most disturbingly, Wal l showed that the actual pairwise correlations implied by modeling the inverse covariance matrices can be contrary to what is intended: pairwise correlations between non-neighbours can exceed those between neighbours. R e g i o n a l i z e d Var iab les M o d e l s The regionalized variables approach to spatial modelling was developed primarily for spa-tial interpolation and spatial prediction applications. A l l R V models are based on unconditional autoregressive formulations and direct modelling of the variance matrix. A n R V model requires an explicit definition of the joint pdf of Ys. The elements of are modeled directly using a covariance function C{hu;o2,6) = Cov(e 8 i ,e B l ) (3.45) where C(-) is a real-valued function of hu = (SJ — s/), the scale parameter (cr2), and a vector of covariance parameters 6. The covariance function should be continuous, non-negative, and 27 monotone decreasing. Wackernagel (1998) provided an extensive list of permissable covariance functions. Two of the simplest, and most commonly used, isotropic covariance functions are: 1. the exponential covariance function C(du) = 2. the Gaussian covariance function C(dit) = { a28d» for du > 0 0 otherwise; a29dl for du > 0 0 otherwise; (3.46) (3.47) where dj* =|| hu ||, a2 > 0, and 9 E [0,1). These covariance functions specify different rates of decay in the magnitude of the correlation between two observations as a function of their separation distance (Figure 3.1). 1.0 0.8 i 0.6 CO 8 0-4 0.2 0.0 —Exp. \ \ \\ Gau. \ \ \ \ v -1 2 3 4 Separation Distance (d) Figure 3.1: Rate of decay of the pairwise covariance between observations as a function of their separation distance under the exponential and Gaussian covariance functions (cr2 = 1; 9 = 0.5). A n advantage of using R V models is that it is possible to interpret the form of the spatial correlation directly from the covariance function. On the other hand, the selection of an appropriate covariance function for a data set is not straightforward. One approach is to base selection on the empirical covariance function 0M=MJ-F5>A (3-48) 1 y V d 1 Md 28 where iSi is an estimator of sSi, Nd is the set of pairs of observations separated by a distance (d), and | Nd | is the number of pairs in Nd (the summation is over all pairs in Nd)- If C(d) is plotted against d, it may reveal an appropriate form for the covariance function; unfortunately, C(d) is a biased estimator of C(d) (Cressie 1993). Another approach is to base the selection of the covariance function on a cross-validation procedure. One observation at a time is excluded and subsequently predicted from the rest of the set, and the covariance function that minimizes the mean squared prediction error is taken as the model covariance function. Both approaches are complicated by heteroskedasticity and correlations among model residuals (Ripley 1981; Zimmerman and Harville 1991). In this thesis, only the exponential and Gaussian covariance functions were considered. The R V models corresponding to these two covariance functions will be referred to as the E X P and G A U models, respectively. 3.3.2 R E M L Estimation In N N and R V models, the vector of observations is modelled as a nonlinear function of the autocorrelation parameter, 9. Under these models, OLS estimators of j3 are not efficient and the OLS estimators of the variance parameters are not consistent. On the other hand, the restricted maximum likelihood ( R E M L ) estimators described by Patterson and Thompson (1971, 1974) have desirable asymptotic and finite sample properties in a spatial setting. Under the assumption of multivariate normality, the log-likelihood of Equation 3.4 is L(f3,a2,9\Ys) = c o - ^ l n l V I - ^ Y s - X / ^ f V - ^ Y s - X ^ ) = Cl-?lincT2-±in\n\-^(Ys-X(3)TCl-l(Ys-Xf3) (3.49) where CQ and c\ are constant terms. The maximum likelihood (ML) estimators of (3 and a2 are Pml = ( x ^ ^ x ) _ 1 X ^ - 1 Y S (3.50) and a2ml = 1(YS - X & r f f i W i Y s - X ^ ) (3.51) where Cl = Cl(8mi) is an estimator of Cl (under the parameterization Cl(0)), and 8mi is the M L estimator of 9. No analytic solution for 9mi exists, so it must be obtained from a (constrained) maximization of Equation 3.49. ' • A disadvantage of M L estimation is that the estimator of a2 (Equation 3.51) is biased in small samples (Judge et al. 1985). The M L estimator tends to underestimate cr2, especially when 29 the number of fixed-effect parameters p is large relative to n. This bias can be attributed to the fact that this estimator does not account for the loss of degrees of freedom incurred in estimating /3. To rectify this problem, Patterson and Thompson (1971, 1974) proposed the R E M L approach. The basis of the R E M L approach is to: (1) estimate the variance parameters from a restricted portion of the likelihood that is confined to the n — p set of linearly independent error contrasts (i.e., contrasts with expectation zero); and then (2) estimate /3 using the estimated generalised least squares (EGLS) formula. Harville (1974) showed that the restricted likelihood required for variance parameter estimation is L*(a2,9 | Y s ) = c 3 - ^ l n a 2 - | In | n | -1 In | X ^ ^ X | - ^ Y ^ Q Y s (3.52) where Q = O - 1 — f i - 1 X ^ X T f i - 1 X J X r f i _ 1 and C3 is a constant. The R E M L estimators of the variance parameters are those values of a2 and 6 that maximize Equation 3.52. There is no closed form solution for the R E M L estimator of 9 but the value of a2, say a2, for which Equation 3.52 is a maximum is given by o2 = —-—YgQY s (3.53) n — p (Harville 1977). Substituting Equation 3.53 into Equation 3.52 yields L*(9\ Y s) = c 4 - In (Y^QYs) - | In | O | - ^ l n I X ^ Y r ^ l (3.54) where L*c(9) is the concentrated restricted log-likelihood function (up to a constant C4) . The R E M L estimator of 9 is then obtained by maximizing Equation 3.54. Maximization of Equation 3.54 is a constrained optimization problem (9 € © C [0,1)). However, L*(-) is a function of 9 alone. This greatly facilitates the implementation of numerical search algorithms. In general, search algorithms requiring a user-specified gradient vector and Hessian matrix are preferable (the first and second partial derivatives of Equation 3.54 are given in Appendix A ) . One impediment to the successful application of search algorithms is the fact that Equation 3.54 can have multiple modes. Under the S A R model, modes will exist for L*c(9) at all values of 9 corresponding to the eigenvalues of W when 9 is not restricted to [0, 1/CJ(1)]. Dietrich (1991) showed that the multimodality of L*c(9) increased up to approximately 10% with increasing sample size under the G A U model. His study also found that L*c(9) was consistently unimodal under the E X P model. In practice, because the concentrated log-likelihood is a univariate function and 9 has a relatively narrow range (0 C [0,1)), visual inspection of the modality of L*(9) is 30 feasible. This would ensure that the global maximum is found. A second option is to repeat the optimization procedure using different starting points. Once Oreml is found, the R E M L estimator of a2 is obtained by replacing ft in Equation 3.53 with ft = ft(9remi). This substitution yields 1 -2 n — p Y s Q Y s (3.55) where Q = ft"1 - fl^X ( x T « - 1 x ) 1 XTft~x. The R E M L estimator of (3 is subsequently ob-tained from the E G L S formula Preml = ( X ^ ^ X ) ' X ^ ^ Y g (3.56) Since the R E M L estimators of a2 and 8 are based on a transformation of the log-likelihood, they have several properties in common with M L estimators. In particular, the R E M L estimators are consistent and asymptotically normally distributed (Cressie and Lahiri 1996). The asymptotic variance matrix is given by the inverse of the matrix J(a\9) = - E ^-eL*(a\9) §,L*(a2,9) d2 n—p 27* tr ( Q § ) (3.57) (Cressie and Lahiri 1996; Appendix A) . In a spatial setting these asymptotic results are conditional on certain assumptions about the form of the joint pdf of Y s and the relative spatial distribution of the data as n —• co. The assumptions concerning the joint pdf of Y s are the usual regularity conditions required for M L estimation (see Cassella and Berger 1990). But Cressie and Lahiri (1996) showed that the consistency and asymptotic normality of the R E M L estimators in a spatial setting also require that the minimum separation distance between observations does not approach 0 as the sample size grows. Finally, note that a 2 e m l = ^ C T 2 ^ for fixed 9, and that a 2 e m l also reduces to the OLS estimator given by Equation 3.15 when ft = In (i.e., when 9 = 0). Although °remi * s asymptotically less precise than a^, it has superior finite-sample bias properties. It is difficult to specify the exact distribution of (3reml fr°m Equation 3.56 because Y s is correlated with ft. However, Breusch (1980) showed that under the assumption of multivariate normality, f3remi is unbiased (provided the expectation exists) regardless of whether the parame-terization of ft is correct. The asymptotic variance of / 3 r e m / is nlim VM(preml) = a2 ( x T f i - 1 x ) _ 1 (3.58) 31 (Judge et al. 1985) and a common estimator of the variance is V^iPreml) = tfend ( X T f T J X ) . (3.59) However, K e n w a r d and Roger (1997) noted that E q u a t i o n 3.59 yields a biased estimator of Var(/3 r e mj) because (1) E q u a t i o n 3.58 does not account for the effect of est imating 9 and a2 on the variance of /3reml; and (2) Var(/3rem;) is a biased estimator of l i m n ^ o o Var(/3 r e m j). In general, E q u a t i o n 3.59 w i l l tend to overestimate the precision of 0remi-Fina l ly , when a l l other factors in the model are fixed, a treatment F - s ta t i s t i c can be com-puted as TTEMLLT (L\T ( x r f W x ) _ 1 A* 1 7 ) _ 1 L f r e m i F(t-l),(n-p) = 77 7^2 (3-60) \Z L)(7reml where tt = fl(9remi). However, the treatment F -s ta t i s t i c given i n E q u a t i o n 3.60 w i l l not s t r ic t ly follow an F(t_i-)(n_p) d is t r ibut ion because this statistic does not account for est imation of 9. K e n -ward and Roger (1997) discuss several approximations that could be used i n place of E q u a t i o n 3.60, the properties of which depend on several factors inc luding the model covariance structure. 3.4 Model Selection Spat ia l analysis methods based on the C E models out l ined i n Section 3.3 offer potent ia l gains i n efficiency and val idi ty relative to A N O V A / O L S . However, these gains are only realized i f the data are spat ial ly autocorrelated and if this autocorrelation is appropriately model led. Legi t -imate problems are ascertaining when spatial autocorrelation is significant, and choosing among compet ing spat ia l C E models. Since spatial autocorrelation is defined differently by different models, measures of spat ia l autocorrelat ion i n a data set are contingent on the identif icat ion of a suitable model . T h e most commonly used tests for spatial autocorrelation are based on N N models. For example, the M o r a n I statistic (Cl i f f and O r d 1981) is a ratio of cross products reflecting the degree of association defined by a spat ia l weight m a t r i x relative to the overall sum of squares: I = n e \ W e ^ (3.61) z eols eols where W = {wu} is the spat ia l weight m a t r i x (wu = 0), e 0 j s is the n x 1 vector of A N O V A / O L S residuals, and z = J2i=i Sl=i wu is a constant. T h i s statistic can be used to test the n u l l hy-pothesis of independence against the alternative hypothesis of non-zero (i.e., posit ive or negative) 32 spatial autocorrelation; if the magnitude of J is sufficiently large, the null hypothesis is rejected. However, the exact distribution of I under the null hypothesis is complicated by autocorrelation among the ANOVA/OLS residuals. Also, the I statistic is conditional on the spatial weight matrix providing a suitable reflection of the pairwise dependence among observations but does not account for estimation of W. A more broadly applicable measure for spatial autocorrelation is the likelihood ratio (LR) statistic. The L R statistic (A) is defined as A = 2 L(J3,a2,9 | Y s ) - L(/3,a2 \YS,9 = 0) (3.62) where L(f3,a2,9 | Y s ) is the unconstrained log-likelihood of Y s and L(/3, a2 | Y s , 8 = 0) is the log-likelihood constrained by the condition 9 = 0. A n L R statistic based on the ratio of the restricted log-likelihoods (A*) can be defined similarly A * = 2 [L*(G2, 9 | Y s ) - L*(a2 | Y s , 9 = 0)] (3.63) where L*(-) is the restricted log-likelihood of Y s - The restricted L R statistic given by Equation 3.63 is easier to compute than A when R E M L estimation is used; the unconstrained function L* (a2,9 \ Y s ) is readily available from the R E M L estimation procedure and L*(a2 | Y s ) is simply a function of a2ls. The restricted L R statistic converges in distribution to a x2 random variable with a single degree of freedom when 9 = 0 (Beckers 1997). This statistic can therefore be used to test the null hypothesis of no spatial autocorrelation (9 = 0) against an alternative hypothesis, the interpretation of which depends upon the spatial C E model under consideration. A drawback of the restricted L R test is that it can only be used as a model selection tool when comparing nested models (i.e., a spatial C E model vs. A N O V A ) . This test cannot be used as a basis for model selection when several different spatial C E models are being considered. Model selection should incorporate the principle of parsimony both in terms of the model itself and in terms of the data given the model. The Akaike Information Criterion (AIC) klC = -L*(a2,8\ Ys) + q (3.64) where q is the number of variance parameters in the analysis model (Akaike 1974), is often used for model selection in time-series settings. The idea is to select the analysis method that yields the smallest A I C . In the context of the analysis methods discussed in this chapter, the A I C penalizes the spatial analysis methods for fitting the additional covariance parameter 9. 33 Chapter 4 Methods 4.1 Simulation Study The second objective of this thesis was to evaluate the efficiency and validity of spatial analysis methods for forestry Nelder experiments. When data are spatially autocorrelated and ANOVA/OLS is used, treatment effects are not estimated efficiently and variance estimators are biased. Analysis methods based on spatial C E models and R E M L estimation may be asymptotically valid and efficient but finite-sample properties can not be established analytically. Therefore, the distributions of the R E M L estimators, and the efficiency and validity of spatial analysis methods, were evaluated by means of a simulation study. The simulation study consisted of four trials. In each trial, a different spatial C E model was used to simulate an autocorrelated series at one of several levels of spatial autocorrelation. The models used were the E X P , G A U , SAR, and C A R models discussed in Section 3.3. Different magnitudes of 9 were simulated in each trial because the level of the spatial autocorrelation implied by 9 is model-specific. Also, this parameter has a restricted range in the two N N models. In each trial, at every level of 9, 500 Nelder experiments were simulated. Since estimation of 9 is invariant to both cr2 and /x (Breusch 1980) these parameters were fixed throughout the study. Every simulated Nelder experiment was analysed using both ANOVA/OLS and a spatial analysis method consisting of the spatial model used in the simulation with R E M L estimation (Table 4.1). The simulated spatial lattice and treatment structures were taken directly from the spacing wheel design of the M K N experiment (Section 2.3). Only a single spacing wheel was simulated, as differences among spacing wheels can be modelled by the mean structure. The spatial weight matrix given by Equation 3.34 was used in the S A R and C A R simulations. In each simulation, 113 normally distributed random deviates were generated using the multiplicative congruential random 34 Table 4.1: Simulation models and analysis methods used in the simulation study. The simulated mean, variance, and spatial autocorrelation are given by /J,, cr2, and 9, respectively. Trial Simulation model a'2 9 Number of simulations Analysis methods 0 500 1 E X P 0 1 0.2 500 ANOVA/OLS, 0.4 500 E X P / R E M L 0.6 500 0 500 2 G A U 0 1 0.2 500 ANOVA/OLS, 0.4 500 G A U / R E M L 0.6 500 0 500 3 S A R 0 1 0.1 500 ANOVA/OLS, 0.2 500 S A R / R E M L 0 500 4 C A R 0 1 0.1 500 ANOVA/OLS, 0.2 500 C A R / R E M L number generator in S A S / I M L (SAS Institute Inc. 1999). These independent normal deviates were then assigned to 25 levels of tree density and transformed through the root of the simulated variance matrix to produce an autocorrelated series. For the S A R model trial, the root of the variance matrix was obtained directly. For the other three simulation models, the Cholesky decomposition (Cressie 1993) was applied to the simulated variance matrix to obtain the necessary transformation. The effects of tree density ( r ) were assigned by dividing the 25 levels of tree density into three classes. The 8 highest density treatments were assigned a value of —0.2, the 8 lowest density treatments a value of +0.2, and the 9 intermediate treatments a value of zero (Table 4.2). This provided the Table 4.2: Simulated treatments effects by tree density class. Class Tree density Simulated treatment (stems/ha) effect (r) 1 > 2500 - 0 . 2 2 800-2500 0.0 3 < 800 +0.2 means to evaluate the effects of spatial autocorrelation on both the power and size of treatment F-tests for different analysis methods. 35 R E M L estimation of covariance parameters was performed using the Newton-Raphson line search algorithm in S A S / I M L (SAS Institute Inc. 1999). The range of §remi under both N N models was restricted to the interval [0, 0.275], as 0.275 corresponded to the reciprocal of the largest eigenvalue of the spatial weight matrix. For the R V models, estimates of 6 were restricted to the range [0, 0.99] since the variance matrix is singular at 9 = 1. The efficiency and validity of each spatial analysis method relative to ANOVA/OLS was assessed following the procedure outlined by Besag and Kempton (1986). Under this scheme, the average empirical variance (EMP) of the set of pairwise treatment contrasts is used to compare the efficiency of each pair of analysis methods in a simulation trial at each level of 6. The empirical variance (Vare) of an estimated treatment contrast (o^^f = f j — TJ) is V a r ^ ^ f ) = — £ [ ( ^ f H ) - (4,,^)] (4.1) fe=i where is a t x 1 vector with a +1 and a —1 in the ith and j t h rows (i,j = 1 ,2, . . . , *) , respectively, and T[fc] is the estimate of r from the kth simulation. Equation 4.1 is the average squared deviation of around its expected value, ^Jij)T ( T is unbiased). E M P is obtained from Equation 4.1 by averaging over all pairwise treatment contrasts — ^ £ E V a r £ ( C > , r ) . (4.2) ^ ' i=l j>i E M P = tAt. - 1 ) —^' . J> Smaller values of E M P correspond to greater precision. Therefore, in each trial, the relative effi-ciency of the spatial analysis method (relative to ANOVA/OLS) could be evaluated at each level of 6 through the E M P statistics. The relative empirical variance (REV) of an analysis method was calculated as E M P . M R E V = JVl x 100% (4.3) i^MfANOVA/OLS where ¥MPM is the E M P of the analysis method and EMPANOVA/OLS is the E M P of ANOVA/OLS. For example, in trial 1, the R E V of E X P / R E M L would be R E V = ™Z EXP,REML x 100% (4.4) E M F ' ANOVA/OLS for a given 0. B y this convention, the R E V of ANOVA/OLS is always 100%. The validity of an analysis method can be assessed by comparing E M P to its estimated value, the average predicted variance (PRE) of the set of all pairwise treatment contrasts. The average predicted variance of a treatment contrast is 500 500 -, V a i <C> f ) = o^o E 4 ] (6h) Xt ( x T S f x ) " X T 6{id)) (4.5) fc=i 36 P R E 2 where fi[fc] = fl(O^); a2^ and 9^ are the estimates of a2 and 9 from the kth simulation, respectively; and Xt defines the linear transformation such that r = Xf/3. P R E is the average of Vai ' (c^^f-) over the t(t — l ) / 2 pairwise treatment contrasts = ^ E » L > - > • (4.6) If an analysis method is valid, estimates of precision are unbiased and P R E should be approximately equal to E M P . If E M P - P R E is positive (negative), it is an indication that the variance of treatment estimates is underestimated (overestimated) (Brownie and Gumpertz 1997). In each trial, the relative variance bias (RVB) of an analysis method at a given level of 9 is R V B = | M ^ - P R E M x l o o % ANOVA/OLS where EMP^vi and P R E ^ f are the E M P and P R E of the analysis method and EMP ' ANOVA/OLS 1 5 the E M P of ANOVA/OLS. The R V B of ANOVA/OLS corresponds to the percent difference between E M P and P R E since E M P ANOVA/OLS is the denominator of Equation 4.7. The distributions of 9, a2, and the treatment F-statistics for both the spatial analysis method and ANOVA/OLS were also of interest. In particular, the distributions of the F-statistics can be used to assess the power and size of resultant F-tests. The treatment effect structure (Table 4.2) provided the opportunity to examine the power and size of different treatment F-tests. The overall treatment effect was partitioned into an among-classes component and a within-classes component (Table 4.3). The significance of each of these components was tested against the experimental error. The empirical rejection rate of the among-classes F-test provided an assessment of its power (i.e., the probability of rejecting the false null hypothesis). Similarly, the empirical rejection rate of the within-classes F-test provided an assessment of its size (i.e., the probability of not rejecting the true null hypothesis). The power of the ANOVA/OLS overall treatment F-test under this treatment effect structure and zero spatial autocorrelation is approximately 0.1 at the 0.05 significance level. In each trial, the empirical rejection rate of the restricted L R test at the 0.1 significance level was calculated for each simulated level of 9. At 9 = 0, this rejection rate indicated the size of this test; for 9^0, the rejection rate indicated the power of the test. For each trial, and each level of 9, the proportion of simulations or rate at which A I C 5 < AlCANOVA/OLS (4.8) where AIQs was the A I C for the spatial analysis method and AIQANOVA/OLS w a s * n e A I C for ANOVA/OLS, was also tracked across the simulated ranges of 9. This provided an indication of the value of this comparison with respect to model identification and selection. 37 Table 4.3: Analysis of variance table for the simulated Nelder experiments. Source of Variation Degrees of freedom Mean Square F-test Treatments 24 Among-classes 2 Within-classes 22 Error 88 Total 112 MSTRT M S A M S B MSERR MSTRT/MSERR MSA/MSERR MSB/MSERR 4.2 Analysis of the MKN Experiment The third objective of this thesis was to assess the practicality of spatial analysis methods for forestry Nelder experiments. ANOVA/OLS, and each of the four spatial analysis methods evaluated in the simulation study (Section 4.1), were therefore applied to data from the M K N experiment (Section 2.3). The analysis of variance table for the M K N experiment age 15 diameter outside bark at 1.3 m above ground (DBH15) data set is given in Table 4.4. The difference between the two Table 4.4: Analysis of variance table for the M K N experiment DBH15 data set. Source of Variation Degrees of freedom Mean Square F-test Spacing wheel 1 Tree density 24 Spacing wheel x tree density 24 Error 163 Total 212 M S S V Y MSTRT M S / J V T MSERR MSSW/MSERR MSTRT/MSERR MSTNT/MSERR spacing wheels was modeled as a fixed effect. As might be expected, the DBH15 data set was incomplete; 3 observations were missing, all from the same spacing wheel. Missing data is a major concern in spacing experiments, because it implies that the original treatment structure has been compromised. It is straightforward to modify an analysis to accomodate missing data when it is known that their absence is due to errors made during collection or recording. But if data may be missing because-trees have died, the analysis must be modified to account for concomitant changes in tree density. In the analyses undertaken here, the assumption was made that changes in tree density resulting from tree mortality would not have had significant effects on tree growth outside the immediate vicinity of the dead trees. Therefore, at each of the three points in the 38 M K N experiment with missing DBH15 values, the nearest inter-spoke and intra-spoke trees were identified and excluded from the analyses. As a result, all five analyses were based on a total of only 213 trees: 113 trees from one spacing wheel, 100 from the other. The spatial analysis methods modeled spatial autocorrelation only within spacing wheels. The S A R / R E M L and C A R / R E M L analysis methods used the spatial weight function given by Equa-tion 3.34. The restricted likelihood functions for the spatial analysis methods were based on a multivariate normal distribution with variance matrix V = <r2n = a 2 « ! 0 0 ft2 (4.9) where Qj is the covariance matrix for the j t h spacing wheel (j = 1,2). R E M L estimation was based on the Newton-Raphson line search algorithm in S A S / I M L (SAS Institute Inc. 1999). The range of possible values for 0remi was restricted to the interval [0, 0.339] for both S A R / R E M L and C A R / R E M L and to [0, 0.99] for E X P / R E M L and G A U / R E M L . 39 Chapter 5 Results and Discussion 5.1 Results The second and third objectives of this thesis were addressed by means of the simulation study and M K N DBH15 analysis described in Chapter 4. The results of these two applications are presented below. 5.1.1 Simulation Study The R E V and R V B of ANOVA/OLS and the spatial analysis method used in each simulation trial are presented in Figure 5.1. Each graph in Figure 5.1 is a plot of R V B vs. R E V for the two analysis methods at the levels of spatial autocorrelation simulated in one of the four trials. The trial numbers and the spatial C E simulation models are indicated above each graph. The simulated level of autocorrelation is referenced by its decimal digit (e.g., 9 = 0.2 is indicated by a 2). The R E V of E X P / R E M L and of G A U / R E M L in trials 1 and 2, respectively, tended to de-crease with the magnitude of the spatial correlation (Figure 5.1). At the highest levels of spatial autocorrelation (9 = 0.6) the E M P statistics for these two spatial analysis methods were approxi-mately 2% smaller than E M P A N O V A / O L S - I n contrast, the R E V of the two analysis methods based on N N models, S A R / R E M L and C A R / R E M L , increased with increasing spatial autocorrelation in trials 3 and 4. ANOVA/OLS was more efficient than these two spatial analysis methods at all levels of e. The R V B of ANOVA/OLS increased with the level of spatial autocorrelation in all four-trials (Figure 5.1). In all trials, PREANOVA/OLS underestimated EMPANOVA/'OLS when 0 was greater than zero. The R V B of E X P / R E M L in trial 1, and the R V B of G A U / R E M L in trial 2, were consistently below zero. The R V B of E X P / R E M L was lowest for 9 = 0.6 (-3.4%) while the R V B of 40 1. EXP 50 40 30 co 20 > or 10 -10 •ANOVA/OLS °EXP/REML 97 oo > or 40 30 20 10 -10 98 99 REV (%) 3. SAR 99 100 101 102 REV(%) i 02 100 101 •ANOVA/OLS "SAR/REML 25 15 m > or •ANOVA/OLS "GAU/REML 97 20 15 10 > or -5 103 104 98 98 2.GAU 99 REV (%) 4. CAR 100 101 •ANOVA/OLS ;;CAR/REML 100 102 104 106 108 REV (%) F i g u r e 5.1: R e l a t i v e v a r i a n c e b i a s ( R V B ) a n d r e l a t i v e e m p i r i c a l v a r i a n c e ( R E V ) o f t h e t w o a n a l y s i s m e t h o d s u s e d i n e a c h s i m u l a t i o n t r i a l . T h e s i m u l a t e d level o f s p a t i a l a u t o c o r r e l a t i o n (8) is i n d i c a t e d b y i ts d e c i m a l d i g i t (e.g., 9 = 0.2 is i n d i c a t e d b y a 2) . E X P / R E M L was lowest for 8 = 0 ( - 1 . 5 % ) . T h e R V B o f SAR/REML a n d t h e R V B of C A R / R E M L i n c r e a s e d w i t h t h e l e v e l o f s p a t i a l a u t o c o r r e l a t i o n i n t r i a l s 3 a n d 4, r e s p e c t i v e l y . F o r these two s p a t i a l a n a l y s i s m e t h o d s , P R E o v e r e s t i m a t e d E M P at 8 = 0 a n d u n d e r e s t i m a t e d E M P at 9 = 0.2. T h e m e a n s a n d s t a n d a r d d e v i a t i o n s of t h e v a r i a n c e p a r a m e t e r e s t i m a t o r s are g i v e n i n T a -b l e 5.1 ( reca l l t h a t cr 2 = 1 i n a l l s i m u l a t i o n s ) . In t r ia ls 1 a n d 2, t h e b i a s i n <5"2js i n c r e a s e d w i t h t h e l e v e l o f s p a t i a l a u t o c o r r e l a t i o n b u t t h e m e a n of t h e R E M L e s t i m a t o r o f a2 w a s c o n s i s t e n t l y close to 1. H o w e v e r , t h e p r e c i s i o n of v2eml w a s less t h a n t h a t o f <r^s a n d d e c r e a s e d as 8 i n c r e a s e d i n these 41 Table 5.1: Mean and standard deviations of the OLS and R E M L estimators of a2 and 9 at each level of spatial autocorrelation (9), by simulation trial (500 replications). Trial 9 &2ols ~ 2 @reml Mean Std.Dev. Mean Std.Dev. •Mean Std.Dev. 1. E X P 0 0.998 0.152 1.011 0.157 0.086 0.123 0.2 0.978 0.147 1.013 0.161 0.207 0.160 0.4 0.922 0.142 1.001 0.175 0.383 0.150 0.6 0.831 0.142 1.013 0.230 0.593 0.114 2. G A U 0 1.000 0.150 1.007 0.152 0.073 0.103 0.2 0.983 0.150 1.002 0.156 0.209 0.129 0.4 0.963 0.152 1.003 0.160 0.409 0.092 0.6 0.923 0.150 1.000 0.165 0.603 0.048 3. S A R 0 0.992 0.155 0.995 0.155 0.034 0.047 0.1 0.989 0.151 0.983 0.148 0.095 0.064 0.2 1.062 0.177 0.982 0.145 0.185 0.052 4. C A R 0 1.002 0.161 1.003 0.162 0.054 0.076 0.1 0.996 0.150 0.991 0.150 0.110 0.094 0.2 0.998 0.148 0.981 0.145 0.164 0.090 two trials. In trials 3 and 4, the mean of the OLS estimator of a2 was close to 1 at all levels of 9. Also, in these last two trials the bias in o 2 e r n l increased, though its precision decreased, with the level of spatial autocorrelation. Overall, the biases in 9remi were greatest at 9 = 0 and, in trials 3 and 4, were also notable at 9 = 0.2. The standard deviation of 9remi was always greater towards the centre of the sampled range of 9 than at the extremeties. The empirical rejection rates for among- and within-class treatment F-tests are summarized in Figure 5.2. Each graph displays the rejection rates for the null hypothesis of no differences among treatment classes on the vertical axis and for the null hypothesis of no differences within treatment classes on the horizontal axis. The trial and simulation model are indicated above each graph, and the results are colour-coded by analysis method. As in Figure 5.1, the level of spatial autocorrelation is indicated by its decimal digit (e.g., 9 = 0 is referenced by the number 0). The rejection rates differed by analysis method within simulation trials (e.g., ANOVA/OLS vs. E X P / R E M L ) because the two analysis methods used different estimates of r and of V. The size of the ANOVA/OLS within-class treatment F-test was badly controlled (Figure 5.2) when 9 > 0 in all four trials. A t the 0.05 significance level, the true null hypothesis was rejected under ANOVA/OLS in over 40% of the simulations in trial 1 at 9 = 0.6. By comparison, the rejection 42 1. EXP 2.GAU 0.6 0.5 t= 0.4 0.3 E < 0.2 0.1 0.55 0.50 0.45 I 0.40 o 0.35 <u in ra 0.30 o 0.25 < 0.20 0.15 0.10 •ANOVA/OLS °EXP/REML 0.0 0.1 0.2 0.3 0.4 Within-class rejection rate 3. SAR 0.5 • 2 • 1 • 0 ' *1 2 •ANOVA/OLS "SAR/REML 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Within-class rejection rate 0.45 0.40 0.35 0.30 o 0.25 < 0.20 0.151 0.00 • 6 • 4 • 2 • 0 2 b O 4 6 •ANOVA/OLS "GAU/REML 0.05 0.10 0.15 0.20 Within-class rejection rate 4. CAR 0.25 0.45 rejection rate o o CJ1 O • 2 • 1 ra 6>°-30 c • 0 1 < "o 0.25 i 2 »ANOVA/OLS °CAR/REML 0.20 1 1 . 0.00 0.05 0.10 0.15 0.20 Within-class rejection rate 0.25 Figure 5.2: Empirical rejection rates (at the 0.05 significance level) of among- and within-class treatment F-tests by analysis method, level of spatial autocorrelation, and simulation trial. The level of spatial autocorrelation (9) is indicated by its decimal digit (e.g., 6 = 0.2 is referenced by the number 2). rates of the same test for the spatial analysis methods were lower and consistently below 0.1. The empirical rejection rates of the false null hypothesis of no differences among treatment classes were higher for ANOVA/OLS in all trials at all levels of spatial autocorrelation (Figure 5.2). This rejection rate for ANOVA/OLS also increased with the level of spatial autocorrelation. The rejection rates of the same test for the spatial analysis methods decreased with the level of spatial autocorrelation. 43 The empirical rejection rate of the restricted L R test, and the proportion of simulations in which AIC^NOVA/OLS exceeded the AIC of the spatial analysis method, are shown in Figure 5.3. The null hypothesis of 8 = 0 (no spatial autocorrelation) and the 0.1 significance level were used for the restricted L R test. 1. EXP 2.GAU e e Figure 5.3: Empirical rejection rates (at the 0.1 significance level) of the restricted L R test for zero spatial autocorrelation, and the rates at which the AIC for ANOVA/OLS exceeded the AIC for the spatial analysis method, in each simulation trial as a function of the spatial autocorrelation (8). The proportion of simulations in which AICANOVA/OLS exceeded the AIC of the spatial analysis method was at least as high as the rejection rate for the restricted L R test (Figure 5.3). Also, when 8 = 0 the restricted L R test consistently rejected the null hypothesis in fewer than 10% 44 of the simulations. Although the rejection rate of the restricted L R test increased with 9 more rapidly in some trials than in others, the scales of the horizontal axes are not strictly comparable. The same value of 6 does not imply the same magnitude of spatial autocorrelation in different models. 5.1.2 M K N A n a l y s i s The ANOVA/OLS analysis of variance table for the DBH15 data set is given in Table 5.2. Note that differences among spacing wheels and differences within tree density treatments between spacing wheels were modeled as fixed effects. Under ANOVA/OLS, interaction and treatment effects were both significant at the 0.05 significance level. The experimental error (i.e., MSERR — aois) corresponded to an estimated coefficient of variation of 22.7% (y = 129.0). Table 5.2: ANOVA/OLS analysis of variance table for the M K N experiment DBH15 data set. Decrees Source of Variation . . & . Mean square F p1 ot freedom Spacing wheel 1 1621.3 L 9 l 0.168 Tree density 24 8673.9 10.24 < 0.001 Spacing wheel x tree density 24 2114.4 2.50 < 0.001 Error 163 846.9 Total 212 1 Probability of observing an F-statistic at least as great under the null hypothesis of no differences among factor effect levels. The frequency distribution of the ANOVA/OLS standardized residuals (i.e., e0is/aois; Fig-ure 5.4) was not significantly different from a normal distribution at the 0.05 significance level (Shapiro-Wilk statistic 0.979; p-value 0.251). The plot of these standarized residuals against the predicted DBH15 value did not reveal any heteroskedastic trends. The estimates of the scale and spatial autocorrelation parameters obtained from each anal-ysis method are given in Table 5.3. Convergence was obtained for the R E M L estimators of 9 for S A R / R E M L and C A R / R E M L . Estimates of 9 for both E X P / R E M L and G A U / R E M L were too close to zero to obtain absolute convergence. The concentrated restricted log-likelihood was examined across the potential ranges of 6 for each spatial analysis method to verify that the R E M L estimates were in the neighbourhood of a global maximum. The restricted L R statistic was used to test the null hypothesis of zero spatial autocorrelation (9 = 0) for each spatial analysis method (Table 5.3). The restricted L R tests indicated that the 45 180 200 220 Class midpoint Predicted DBH 15 (cm) Figure 5.4: Frequency distribution and residual plot for the standardized ANOVA/OLS M K N DBH15 residuals. DBH15 data did not exhibit significant positive spatial autocorrelation, as defined by any of the four simulation models applied. The restricted L R statistic for C A R / R E M L was negative, indicating that the resultant C A R model was a worse fit for the data than the i . i .d . A N O V A model. ANOVA/OLS also had the smallest A I C among the five analysis methods applied (Table 5.3). Table 5.3: Estimates of scale (a 2) and spatial autocorrelation (6) pa-rameters, A I C , and restricted L R (A*) statistics (with the associated p-values) for the M K N DBH15 data set, by analysis method. Analysis method a2 e A I C A* P1 ANOVA/OLS 846.9 - 781.72 - -E X P / R E M L 846.9 1.36 x 1 0 - 1 8 782.72 0.000 1.000 G A U / R E M L 846.9 9.88 x 1 0 - 1 9 782.72 0.000 1.000 S A R / R E M L 855.7 0.045 782.40 0.645 0.422 C A R / R E M L 873.5 0.106 782.99 -0.522 N / A Probability of observing an A* statistic at least as great under the null hypothesis of 0 = 0. Estimates of the effects of tree density differed little among the analysis methods applied. Rankings of the tree density effects as estimated by ANOVA/OLS, E X P / R E M L , G A U / R E M L , and S A R / R E M L were indentical (Table 5.4). However, rankings under the C A R / R E M L analysis method were different. Only the average treatment effect rankings are listed in Table 5.4. Due to the significant interaction effect for tree density by spacing wheel, the within spacing wheel treatment effect rankings may be different. 46 Table 5.4: Rankings of estimated tree density effects for the M K N DBH15 dataset, by analysis method (largest effect is rank 1; ranks different from ANOVA/OLS are in bold). Tree density Estimated effect I ank (stems/ha) ANOVA/OLS E X P / R E M L G A U / R E M L S A R / R E M L C A R / R E M L 17,001 23 23 23 23 23 12,321 16 16 16 16 16 8,931 24 24 24 24 24 6,475 25 25 25 25 25 4,693 18 18 18 18 18 3,401 22 22 22 22 22 3,089 7 7 7 7 7 3,042 14 14 14 14 14 2,464 17 17 17 17 17 2,348 21 21 21 21 21 1,787 20 20 20 20 19 1,607 13 13 13 13 13 1,597 11 11 11 11 11 1,295 15 15 15 15 15 949 19 19 19 19 20 875 10 10 10 10 10 850 8 8 8 8 9 690 12 12 12 12 12 504 9 9 9 9 8 450 4 4 4 4 4 368 5 5 5 5 5 267 6 6 6 6 6 242 3 3 3 3 3 134 2 2 2 2 2 131 1 1 1 1 1 5.2 Discussion The simulation study described in Section 4.1 was undertaken to answer three questions. The first was was whether the spatial analysis methods described in Chapter 3 were more efficient than ANOVA/OLS when applied to spatially autocorrelated Nelder experiments. The second question was whether the estimates of precision supplied by the spatial analysis methods were unbiased. The last question was whether the restricted L R test, or the A I C , could be used to determine the presence or significance of spatial autocorrelation in a data set. The analyses of the M K N DBH15 data set were undertaken to assess the practical properties of spatial analysis methods relative to ANOVA/OLS. 47 5.2.1 Simulation Study The two spatial analysis methods based on R V models, E X P / R E M L and G A U / R E M L , were slightly less efficient than ANOVA/OLS at 9 = 0 (Figure 5.1). E X P / R E M L was also less efficient than ANOVA/OLS at 9 = 0.2. As the spatial autocorrelation was increased, however, the relative efficiency of these two spatial analysis methods increased. This result was consistent with findings from other published studies, but the observed degree of difference in precision was relatively minor, and deserves some attention. Neither the E X P / R E M L nor G A U / R E M L analysis methods offered much more than a 2% gain over ANOVA/OLS in the efficiency of treatment effect estimators over the range of simulated autocorrelation (Figure 5.1). Zimmerman and Harville (1991) reported much larger gains in the relative efficiency of spatial analysis methods based on R V models. Their study was based on agricultural uniformity tr ia l 1 data containing strong large-scale spatial trends. They surmised that some of the apparent efficiency of the spatial analysis methods relative to ANOVA/OLS in their study could be attributed to the potential for these analysis models to "soak up" some of this large-scale spatial variation. No large-scale fixed trends were simulated in this study. Brownie and Gumpertz (1997) also reported more appreciable increases in the efficiency of treatment effect estimators of analysis methods based on RV models relative to ANOVA/OLS with increasing spatial autocorrelation. Their study was simulation-based, but the simulated field experiment was constructed on a regular spatial grid and treatments were randomly allocated to experimental units. Their study showed that E X P / R E M L could be as much as 40% more efficient than ANOVA/OLS. The discrepancy between the results of the present study and those of Brownie and Gumpertz may be due in part to the differences between the irregularly spaced but more symmetrical systematic Nelder experiment simulated here, and the randomized rectangular field experiment simulated by Brownie and Gumpertz. This point is considered further below. The simulation study also indicated that the analysis methods based on N N models are relatively less efficient than ANOVA/OLS (Figure 5.1). Furthermore, the relative efficiency of S A R / R E M L and C A R / R E M L decreased with the magnitude of the spatial autocorrelation. This was a surprising and somewhat disconcerting result given the popularity of these two analysis methods. As Wall (2000) noted, these analysis methods have achieved sufficient popularity that the statistical software package S-Plus© has built-in functions for fitting S A R and C A R models. The efficiency of ANOVA/OLS relative to these two spatial analysis methods in the present study •"•A field experiment conducted with an uniformly applied treatment. 48 may have been due in part to the irregular spatial distribution of the observations in the simu-lated Nelder experiments. These N N models were originally conceived of as models for processes distributed over regular, infinite lattices (Whittle 1954, Besag 1974). The poor performance of these analysis methods may also have been due to the complex nature of the correlation structure implied by the S A R and C A R models. Because the correlation structure is modelled through the inverse variance matrix, the resultant variance matrices do not generally have constant diagonal elements or even similar covariances among pairs of observations separated by the same distance (Wall 2000). This inherent lack of stationarity was built into the S A R and C A R simulations. Regardless of their relative efficiency, however, unless estimates of the precision of treatment contrasts are unbiased, more complex analysis methods are of limited value. The simulation study demonstrated that, as expected, ANOVA/OLS is not a valid analysis method when data are spa-tially autocorrelated. ANOVA/OLS underestimated the average precision of the treatment contrast estimators when 9 was greater than zero, and the bias increased with the spatial autocorrelation. Brownie and Gumpertz (1997) also examined the validity of ANOVA/OLS through the dis-tribution of E M P and P R E statistics. In their study, EMPANOVA/OLS - ^^^ANOVA/OLS WAS negative when the data were spatially autocorrelated, implying that the variance estimators were positively biased. This seemingly contradictory result can again be explained by the different treat-ment structure imposed in their simulations. Brownie and Gumpertz noted that PBEANOVA/OLS was greater than EMPANOVA/OLS because of an induced block effect. This induced block effect resulted from the simulated spatial autocorrelation and a restricted randomization of treatments to experimental units. Since no block effects were included in their analysis models, the induced block effect was allocated to the experimental error. In the present study, spatial trends may have been induced by the simulated spatial autocorrelation. However, because the simulated treatment structure was constant throughout, the spatial trend was divided into variability within and among treatments. ANOVA/OLS was not valid in any of the simulation trials when 9 was greater than zero, but the validity of the spatial analysis methods varied by trial (Figure 5.1). In trials 1 and 2, appli-cation of the spatial analysis method corresponding to the simulation model resulted in relatively minor differences between E M P and P R E . These results were consistent with those of Brownie and Gumpertz (1997). The R E M L estimators of 9 and a2 also performed well in the first two trials (Table 5.1). In trials 3 and 4, the spatial analysis methods were not valid and the R E M L estimators of a2 were biased. In these two trials, a2ls was a better estimator of a2 than o2eml . The simulated spacing wheel is an example of a case in which R E M L estimation of the scale 49 parameter is preferable to M L estimation. The number of fixed effect parameters (p — 25) in the analysis models was large relative to total number of observations (n = 113). In this case, the finite sample bias of the M L estimator of a2 E ( ^ ) - a 2 = - V = -^a 2 (5.1) became relatively large. The R E M L estimator for a2 is unbiased and, in the absence of spatial autocorrelation, reduces to the usual OLS estimator. The R E M L estimators of 9 were consistently biased at 9 = 0 (Table 5.1). In a few of these simulations, the R E M L estimation optimization procedure failed to converge because the Newton-Raphson line search algorithm employed in S A S / I M L (SAS Institute Inc. 1999) did not strictly adhere to the restrictions placed on the range of possible values for 9. As the estimate of 9 approached zero, the algorithm would sometimes attempt to evaluate the gradient or Hessian at values of 9 < 0. Also, in some simulations, the algorithm did not possess sufficient sensitivity to optimize the objective function because the log-likelihood was very flat in the neighbourhood of 0. Notwithstanding, the unit range and unidimensional nature of the objective function permit more direct search operations for any single experimental data set. Although the differences between E M P and P R E were used to determine the validity of the analysis methods, validity is more commonly understood in terms of the distribution of the treatment F-statistics. The empirical rejection rate or size of the ANOVA/OLS within treatment class F-test was badly controlled when spatial autocorrelation was simulated (Figure 5.2). The true null hypothesis of no differences among treatments within classes was rejected by the ANOVA/OLS within-class F-test in an increasing proportion of the simulations as the value of 9 was increased, regardless of the simulation model. This result is directly related to the underestimation of the scale parameter (a2; Table 5.1) and of E M P (Figure 5.1). These biases are also responsible for the increasing empirical rejection rates of the among-classes treatment F-test for ANOVA/OLS as the spatial autocorrelation was increased (Figure 5.2). The false null hypothesis of no differences among treatment classes was rejected more often than it should have been (given the actual treatment and correlation structures) because the variance of treatment contrasts was underestimated. Figure 5.2 should not be taken to imply that the ANOVA/OLS treatment F-tests are more powerful when data are spatially autocorrelated: the rejection rate increased with 9 because the variance estimators were increasingly biased. The rejection rate for the among-classes treatment F-test would increase regardless of the simulation structure (i.e., even if r = 0). Unfortunately, the rejection rate of the within-class treatment F-test for ANOVA/OLS was not 0.05 (the significance level used) when 50 9 = 0 (Figure 5.2). This was an artefact of the simulation study and suggests that the number of simulations was too small. 2 The empirical rejection rate of the within-classes treatment F-test for the spatial analysis methods was relatively well-controlled at around 0.05 (the significance level used; Figure 5.2). This was particularly true in the first two trials. What was more interesting was that the rejection rate of the among-classes treatment F-test for all four spatial analysis methods decreased as 9 increased. The power of the treatment F-test is generally regarded as a function of the variance (cr2); the effect size (5.2) cr where T/J) is the jth largest treatment effect; the sample size (ra); and the number of observations per treatment (Nemec 1991). Clearly, this is an oversimplification; all these factors were held constant in this study. First, the measure of effect size given by Equation 5.2 is too liberal when data are autocorrelated. The variability of an autocorrelated series can not be adequately described by cr 2 alone. Secondly, as Griffiths (1992) pointed out, spatial autocorrelation impacts on the information content or dimensionality of a sample of a given size. A more positively autocorrelated sample can therefore be interpreted as having a smaller effective sample size. 3 This will impact on the degrees of freedom of the treatment F-test. In the simulation study, the efficiency and validity of the spatial analysis methods were evaluated only where the true model for the spatial autocorrelation was known and exploited. In practice, it wil l be necessary to ascertain which model is most appropriate for a given data set. The simulation study indicated that both the restricted L R test and the AICs provide useful measures for spatial autocorrelation and model selection (Figure 5.3). However, in this simulation study, the comparison of AICs is essentially the same as the restricted L R test. The proportion of simulations in which A I C A N O V A / O L S exceeded the A I C of the spatial analysis method was always at least as high as the empirical rejection rate for the restricted L R tests because of the manner in which the simulation study was constructed. When the AICs of nested models are compared, the comparison can be written as a ratio of the restricted likelihoods, that is, if A I C s - A I C A N O V A / O L S < 0 2This discrepancy between nominal and empirical rejection rates may alternatively be ascribed to the imperfections of the random number generator. However, the quality of the random number generator in SAS/IML (SAS Institute Inc. 1999) has been well substantiated by Fishman and Moore (1982). interestingly, the reverse is also true; negatively autocorrelated data of dimension n carry a greater information content and have an effective sample size greater than n (Griffiths 1992). 51 where AIQs is the A I C of one of the spatial analysis methods used in the simulation study, then from Equation 3.64 -L*(a2,9\ Y s ) + 2 + L * ( a 2 | Y s , 8 = 0) - 1 < 0 (5.3) and from Equation 3.63 -0 .5A* + 1 < 0 (5.4) A* > 2 Therefore, an A I C comparison among nested models is equivalent to a restricted L R test with a significance level of approximately 0.16. The value of the A I C lies more in its potential for comparing non-nested models for spatial autocorrelation. Overall, the simulation study indicated that the analysis methods based on R V models were efficient and valid relative to ANOVA/OLS for the simulated Nelder experiment. The study also indicated that S A R / R E M L and C A R / R E M L were neither efficient nor valid for this same application. However, it must be emphasized that these results are conditional on the application of the correct model for spatial autocorrelation. Also, as mentioned above, the fact that the empirical rejection rates of the within treatment class F-tests for ANOVA/OLS deviated from the nominal a level of 0.05 indicated that it may have been useful generate more simulations. Even with such an improvement, however, Ripley (1981; p. 17) has pointed out that "a problem with using simulation to find empirical distributions is that one is usually interested in the tails of the distribution and most of the simulations provide very little relevant information." It can thus be argued that uniformity trial data would have provided a more secure basis for evaluating different analysis methods. For the purposes of this thesis, however, uniformity trial data are unobtainable; it is impossible to conduct a Nelder experiment, or any other spacing experiment, without imposing the effects of differential tree density levels. Simulation experiments are the only option available for evaluating the statistical properties of different analysis methods for Nelder experiments. Of all the considerations that should be weighed in evaluating the results of this study, perhaps the most important is the fact that all simulations were based on a single experimental design and layout. Judge et al. (1985) noted that simulation results are often highly conditional on the design matrix ( X in Equation 3.4). This restriction in the scope of this study is of added concern here for two reasons. First, the effects of spatial autocorrelation will certainly be a function of the spatial scale of the experiment. If the simulation study were repeated with the same values of 8 but different dimensions for the simulated spacing wheel, the effects of spatial autocorrelation on the efficiency and validity of a given analysis method could be quite different. Second, Nelder 52 spacing wheels are relatively symmetrical, and the symmetry of a design can have a major impact on the relative efficiency of ANOVA/OLS. Prior to undertaking this study, a simple type l a spacing wheel (Figure 2.2) was considered as a candidate experimental design for the simulations. However, it was soon discovered that as a consequence of the perfect symmetry in the type l a (or type lb) spacing wheel, the OLS and G L S estimators of /3 were equivalent under all the isotropic spatial C E simulation models used. 4 The equivalence of the OLS and G L S estimators of /3 in symmetrical designs is due to the form of the G L S estimator. This estimator is essentially a weighted least-squares estimator. When the data are spatially autocorrelated under one of the isotropic spatial C E models described in Chapter 3, the G L S estimator weights the observations by their relative isolation; observations that have many close neighbours are given a smaller weight than observations that are relatively remote. But in a spatially symmetrical design like the type la spacing wheel, all observations assigned to a given treatment have exactly the same relative spatial position. Thus, all observations on a treatment level must be weighted equally and the resulting G L S estimator of T wil l be equivalent to the OLS estimator. It follows that under these circumstances the OLS estimator of r is efficient and attains the Cramer-Rao lower bound. The OLS estimators of variance parameters wil l still be biased, however. The spacing wheels used in the M K N experiment were not perfectly symmetrical, but were still highly symmetrical relative to the designs simulated in other published studies. This may partially explain why the efficiencies of E X P / R E M L and G A U / R E M L were not found to be as large, relative to ANOVA/OLS, as those reported by Brownie and Gumpertz (1997). 5 .2 .2 M K N Analysis Considerable variability was observed in the M K N DBH15 data set, but the data did not appear to exhibit spatial autocorrelation. Based on the empirical distribution of the ANOVA/OLS residuals (Figure 5.4), the normal distribution appeared to be a suitable distribution on which to base R E M L estimation of variance parameters for the spatial analysis methods. For the M K N DBH15 data set, the R E M L estimation procedure for the E X P / R E M L and G A U / R E M L analysis methods failed to attain convergence. Visual inspection of the restricted log-likelihood functions confirmed that the estimates of 9 were in the neighbourhood of a global maximum. The estimates of 9 for these two spatial analysis methods were very close to zero, 4 A formal proof of this assertion is given in Appendix B. 53 so the corresponding estimates for a2 were nearly identical to those given by ANOVA/OLS (Ta-ble 5.3). Also, both these spatial analysis methods produced the same treatment effect rankings as ANOVA/OLS (Table 5.4). Since the estimates of spatial autocorrelation were minimal for both E X P / R E M L and G A U / R E M L , all observations received approximately the same weight (i.e., as Oreml ~* 0, f i -» In) and fremi « T0is. The R E M L estimation procedure did converge for the S A R / R E M L and C A R / R E M L analysis methods. Visual inspection of the restricted log-likelihood functions confirmed that the positive estimates for 9 were in the neighbourhood of a global maximum. Based on the restricted L R tests, the estimates of 9 for these two spatial analysis methods were not significantly different from zero (Table 5.3). However, the value of 9Teml for C A R / R E M L was sufficiently large that the resulting estimates of T were different from ANOVA/OLS to an extent that different overall rankings were produced (Table 5.4). The restricted L R tests indicated that none of the spatial analysis methods proved any more useful than ANOVA/OLS in explaining the variation in this dataset (Table 5.3). Based on the minimum A I C rule, ANOVA/OLS provided the most concise stochastic description of the data (Tables 5.3). Application of the four spatial analysis methods was a much more involved procedure than the application of ANOVA/OLS. Although the spatial coordinates of the trees in a Nelder spacing wheel are relatively easy to obtain, the construction of the spatial weight matrix ( W ) for N N models can be a difficult task, especially for large asymmetric designs with missing observations. Also, the spatial analysis methods require the derivation and differentiation of the appropriate restricted log-likelihood functions. Finally, although numerical search procedures are built into many statistical software packages, R E M L estimation will generally require some programming skills. Nonetheless, for the spatial analysis methods examined here, the actual R E M L estimation procedure is quite efficient in practice. There is only a single parameter estimate to search for, and a direct inspection of the result is facilitated by the restricted range for 9remi. 54 Chapter 6 Conclusion Spatial analysis methods are most useful for data that exhibit autocorrelation when variable values are ordered in geographic space. Classical analysis methods rely on the device of experimental randomization to neutralize the effects of autocorrelation over repeated trials. But in any single field experiment, the application of spatial analysis methods to spatially autocorrelated data can provide more efficient estimates of treatment effects. Nelder (1962) spacing wheels are economical designs for studying tree density-growth re-lationships. The formal geometric structure of these designs also lends them to spatial analysis applications. In this thesis, a simulation study was used to assess the efficiency and validity of four spatial analysis methods relative to ANOVA/OLS. In the simulation study, the two spatial analysis methods based on R V models and R E M L estimation tended to be both valid and more efficient than ANOVA/OLS when spatial autocorrelation was simulated. The reverse was true for the two spatial analysis methods based on N N models. The application of the spatial analysis methods to the M K N DBH15 data did not yield any additional information about the variability in the variable of interest or of the experimental material itself. Estimates of spatial autocorrelation parameters for all four spatial analysis methods were low relative to the ranges of spatial autocorrelation used in the simulation trials. However, this application did demonstrate the complexities of spatial modeling and likelihood-based estimation for field experiments with large numbers of irregularly distributed observations. These issues present new challenges, but given the power and flexibility of modern statistical software, it is certainly feasible to implement spatial analytic techniques. Whether spatial autocorrelation will have a significant impact on the precision of treatment effect estimates, or on the accuracy of the estimates of their precision, in actual forestry Nelder experiments could not be determined. However, the levels of spatial autocorrelation estimated for 55 the M K N DBH15 data should not necessarily be taken as typical or representative of the levels found in a more general class of Nelder experiments. Galinski et al. (1994) reported significant spatial dependence in Scots pine (Pinus sylvestris L.) height growth patterns in a Nelder experiment in Sweden. Significant spatial autocorrelation has also been reported in agricultural trials operated at spatial scales in line with the Nelder experiments examined here (e.g., Bal l et al. 1993; Brownie et al. 1993). The presence and significance of spatial autocorrelation in Nelder experiments must be determined on a case-by-case basis. This thesis has provided some indication that L R tests or information criteria may be useful tools for discerning the significance of spatial autocorrelation in a data set and for selecting an appropriate spatial C E model. Most importantly, the results of this study indicate that R V models may be the most suitable class of spatial C E models for forestry field experiments. The spatial analysis methods based on N N models studied in this thesis were neither valid nor efficient, and were more difficult to implement and interpret. 56 Literature Cited Akaike, H . 1974. A new look at the statistical model identification. I E E E Trans. Automatic Control AC-19: 716-723. Anselin, L . 1988. Spatial Econometrics: Methods and Models. Dordrecht: Kluwer. 284 pp. Antle, A . N . and P.L. Marshall. 1996. Fill ing in missing forestry data: exploring autocorrelation techniques. In Proceedings Spatial Accuracy Assessment in Natural Resources and Environ-mental Sciences: Second International Symposium, May 21-23, Fort Collins, C O . U S D A For. Ser. Gen. Tech. Rep. RM-GTR-277. pp. 593-600. Azais, J . - M . , H . Monod, and R . A . Bailey. 1998. The influence of design on validity and efficiency of neighbour methods. Biometrics 54: 1374-1387. Bal l , S.T., D . J . Mulla , and C F . Konzak. 1993. Spatial heterogeneity affects variety trial inter-pretation. Crop Sci. 33: 931-935. Beckers, F . 1997. On the Spatial and Space-Time Analysis of Field Experiments. Ph .D. thesis, Universite Catholique de Louvain, Belgique. 178 pp. Besag, J . 1974. Spatial interaction and the statistical analysis of lattice systems. J . Royal Statist. Soc. B 36: 192-236. Besag, J .E. and R . A . Kempton. 1986. Statistical analysis of field experiments using neighbouring plots. Biometrics 42: 231-251. Breusch, T.S. 1980. Useful invariance results for generalized regression models. J . Econometrics 13: 327-340. Brownie, C , D . T . Bowman, and J .W. Burton. 1993. Estimating spatial variation in analysis of data from yield trials: a comparison of methods. Agron. J . 85: 1244-1253. 57 Brownie, C. and M . L . Gumpertz. 1997. Validity of spatial analyses for large field trials. J . Agiic . Biol . Environ. Stat. 2:1-23. Casella G . and R . L . Berger. 1990. Statistical Inference. Belmont, C A : Duxbury. 650 pp. Chacko, V . J . 1965. Designs for forestry spacing and thinning experiments. In Proceedings of the Conference of Advisory Group of Forest Statisticians of the International Union of Forest Research Organisations Section 25, 27th Sept. to 1st Oct., Stockholm, Sweden. Paper No.5. Chung, K . L . 1974. A Course in Probability Theory. 2nd ed. New York: Academic Press. 365 pp. Cliff, A . D . and J . K . Ord. 1981. Spatial Processes: Models and Applications. London: Pion. 266 pp. Cole, E . C . and M . Newton. 1987. Fifth-year response of Douglas-fir to crowding and nonconiferous competition. Can. J . For. Res. 17: 181-186. Cressie, N . 1990. The origins of kriging. Math. Geol. 22: 239-252. Cressie, N . A . C . 1993. Statistics for Spatial Data. New York: Wiley. 900 pp. Cressie, N . and S.N. Lahiri . 1996. Asymptotics for R E M L estimation of spatial covariance pa-rameters. J . Statist. Plann. Inference 50: 327-341. Dietrich, C R . 1991. Modality of the restricted likelihood for spatial Gaussian random fields. Biometrika 78: 833-839. Evert, F . 1971. Spacing studies - A review. Forest Management Institute Information Report FMR-X-37 . Canadian Forestry Service: Ottawa. 95 pp. Evert, F. 1984. A n update (1970/71-82) of the annotated bibliography on initial tree spacing. Petawawa National Forestry Institute Information Report PI-X-44. Canadian Forestry Ser-vice: Chalk River, Ontario. 157 pp. Fishman, G.S. and L . R . Moore. 1982. A statistical evaluation of multiplicative congruential random number generators with modulus 2 3 1 — 1. J . Amer. Statist. Assoc. 77: 129-136. Freeman, G . H . 1964. The use of a systematic design for a spacing trial with a tropical tree crop. Biometrics 20: 200-203. 58 Galinski, W . , J . Witowski, and M . Zwieniecki. 1994. Non-random height pattern formation in even aged scots pine (Pinus sylvestris L.) Nelder plots as affected by spacing and site quality. Forestry 67: 49-61. Gould, P. 1970. Is statistix inferens the geographical name for a wild goose? Econ. Geog. 46: 439-448. Griffith, D . A . 1992. What is spatial autocorrelation? Reflections on the past 25 years of spatial statistics. L'Espace Geographique 3: 265-280. Hansen, M . H . , W . G . Madow, and B . J . Tepping. 1983. A n evaluation of model-dependent and probability sampling inferences in sample surveys. [With comments and rejoinder.] J . Amer. Statist. Assoc. 78: 776-807. Harville, D . H . 1974. Bayesian inference for variance components using only error constrasts. Biometrika 61: 383-385. Harville, D . H . 1975. Experimental randomization: Who needs it? Amer. Statist. 29: 27-31. Harville, D . H . 1977. Maximum likelihood approaches to variance component estimation and to related problems. J . Amer. Statist. Assoc. 72: 320-338. Hooper, P . M . 1989. Experimental randomization and the validity of normal-theory inference. J . Amer. Statist. Assoc. 84: 576-587. Judge, G . G . , W . E . Griffiths, R . C . H i l l , H . Lutkepohl, T . - C . Lee. 1985. The Theory and Practice of Econometrics. 2nd ed. New York: Wiley. 1019 pp. Kenward, M . G . and J . H . Roger. 1997. Small sample inference for fixed effects from restricted maximum likelihood. Biometrics 53: 983-997. Loo-Dinkins, J .A . and C . G . Tauer. 1987. Statistical efficiency of six progeny test field designs on three loblolly pine (Pinus taeda L.) site types. Can. J . For. Res. 17: 1066-1070. Magnussen, S. 1993. Bias in genetic variance estimates due to spatial autocorrelation. Theoiet. App. Gen. 86: 349-355. Milne, A . 1959. The centric systematic area-sample treated as a random sample. Biometrics 15: 270-297. 59 Mora, A . 1996. Metodo de ajuste de un modelo rendimiento-densidad para observaciones prove-nientes del diserio sistematico de Nelder. [In Spanish.] Rev. Facultad de Agronomia 22: 31-36. Namkoong, G . 1966. Applications of Nelder's designs in tree improvement research. In Pro-ceedings 8th Southern Conference on Forest Tree Improvement, June 16-17, 1965, Savannah, Georgia, pp. 24-37. Nelder, J . A . 1962. New kinds of systematic designs for spacing experiments. Biometrics 18: 283-307. Nemec, A . F . L . 1991. Power analysis handbook for the design and analysis of forestry trials. B . C . Ministry of Forests Biometrics Information Handbook 2. Victoria, B . C . : Crown Publications. 26 pp. Nyland, R . D . 1996. Silviculture: Concepts and Applications. New York: McGraw-Hil l . 633 pp. Paltzat, R. 1984. Performance of red alder in Nelder plot designs at age 7 at the U . B . C . Research Forest. B.S.F. thesis, University of British Columbia, Vancouver, B . C . 45 pp. Papadakis, J.S. 1937. Methode statistique pour des experiences sur champ. [In French.] Bul l , de l'lnstitut d'Amelioration des Plantes a Salonique 23: 13-27. Patterson, H . D . and R. Thompson. 1971. Recovery of inter-block information when block sizes are unequal. Biometrika 58: 545-554. Patterson, H . D . and R. Thompson. 1974. Maximum likelihood estimation of components of vari-ance. In Proceedings of the 8th International Biometrics Conference, Aug. 25-30, Constanta, Romania, pp. 197-207. Radosevich, S.R. 1987. Methods to study interactions among crops and weeds. Weed Tech. 1: 190-198. Reukema, D . L . 1979. Fifty-year development of Douglas-fir stands planted at various spacings. U S D A For. Serv. Pap. PNW-253. 21 pp. Ripley, B . D . 1981. Spatial Statistics. New York: Wiley. 252 pp. Rustagi, K . P . 1984. A modified Nelder design to achieve square or triangular spacing. Indian Forester 110: 521-528. 60 SAS Institute Inc. 1999. SAS OnlineDoc. V.8. SAS Institute Inc.: Gary, N C . Smith, J . H . G . 1959. Comprehensive and economical designs for studies of spacing and thinning. For. Sci. 5: 237-245. Smith, J . H . G . 1978. Design factors from Nelder and other spacing trials to age 20. Com. For. Rev. 57: 109-119. Student. 1914. The Elimination of Spurious Correlation due to Position in Time or Space. Biometrika 10: 179-180. Student. 1937. Comparison between balanced and random arrangements of field plots. Biometrika 29: 363-379. Thornett, M . L . 1982. The role of randomization in model-based inference. Aust. J . Stat. 24: 137-145. Toumey, J .W. and C F . Korstian. 1943. Seeding and Planting in the Practice of Forestry. 3rd ed. Wiley, New York. 520 pp. Wall , M . M . 2000. A close look at the spatial correlation structure implied by the C A R and S A R models. Res. Rep. No. 2000-022, Division of Biostatistics, University of Minnesota, Minneapolis, M N . 16 pp. Wackernagel, H . 1998. Multivariate Geostatistics. 2nd ed. Berlin: Springer. 291 pp. Whittle, P. 1954. On stationary processes in the plane. Biometrika 41: 434-449. Zedaker, S .H. 1982. Growth and development of young Douglas-fir in relation to intra- and inter-specific competition. Ph .D. thesis, Oregon State University, Corvallis, Oregon. 175 pp. Zimmerman, D . L , and D . A . Harville. 1991. A random field approach to the analysis of field-plot experiments and other spatial experiments. Biometrics 47: 223-239. Zyskind, G . 1967. On canonical forms, non-negative covariance matrices and best and simple least squares linear estimators in linear models. Ann. Math. Stat. 38: 1092-1109. 61 APPENDICES 62 Appendix A R E M L Estimation Consider the statistical model for the n x 1 random vector Y Y = X /3 + e ( A . l ) where X is a n x p matrix of full column rank, e is a n x 1 vector of random errors, e ~ N (0, V ) , and V is a n x n positive definite matrix. Further, let V = a2 ft where a2 is a general scale parameter (a2 > 0) and ft = ft(8) is a n x n positive definite matrix known up to a parameter 8. The restricted log-likelihood is L*(*2,8 | Y ) = c - In a 2 - 1 In | ft \ - ± In | X T f i " 1 X | - ^ 2 Y r Q Y (A.2) A B C D where Q = ft-1 — f i _ 1 X ( x T f i _ 1 x ) X T f i - 1 , c is a constant term, and £*(•) is divided into the components A, B, C , and D, for ease of reference below (Harville 1974; Section 3.3.2). Differentiation of the Restricted Log-Likelihood The first two partial derivatives of the restricted log-likelihood with respect to the scale parameter (a2) are and £ l V , » | Y ) - ^ - ^ < J Y . ( A . 4 ) Differentiation of L*(-) with respect to the covariance parameter (8) requires partial differ-entiation of the inverse and determinant of the matrix ft. The subsequent derivations make use of the following results 63 i . j ^ - i = _ n - i (&n) ft-1 n 1 ( * « ) ] In I ft 1= tr 3. ^ Y ^ Q Y = - Y r Q (^n) Q Y the latter being due to Beckers (1997). Making use of these results, the first partial derivative of L*(-) with respect to 9 can be obtained from •2-A = 0 & R - A aeu - A. — ae = i t r i n - 1 = ^ [n- 1 (An) (x^- ix)" i ( x ^ x ) (x^^x)" 1 x^ft-1 (&n) n-ix The first partial derivative of the L*(-) with respect to 9 is then rtr 1 + 2 t r . + i Y T Q ( c l n ) Q Y (xTft-xx) 1 x T n - 1 (^n) n - 1 x (A.5) The second partial derivative of the restricted likelihood with respect 9 can be obtained from tr [n- 1 (&n)] = - t r [ f t - 1 (&n) ft-1 (^n)J + tr [ f t - 1 (£n); a. de tr = tr (xTft- 1x) _ 1xTft- 1 (|ft) ft-xx x ^ x ^ x ^ f t - 1 (An) ft-1x(xrft-1x)-1xTft-1 (An) ft-xx ( x ^ ^ x ) - ^ ^ - 1 (An) n - 1 (An) ft*x "(x^ft^X^X^YT 1 ( J j r f t ) ft-XX - 2tr + tr tr (xrn-!x) 1 x^ft-1 (^ft) ft-1 (j|ft) ft-ix 2tr (xTft-Jx) 1 X^ft- 1 (^ft) ft"1 (^ft) ft-iX 64 + tr - t r + tr (x^ x^) 1 x T n - 1 (jjpri) n-1* (X^TJ^X)" 1 x^ -1 ($n) n-1 (£n) n^x (x^ x^)-1 xTn-1 (J^n) n-ix| A ^ Q ( i n ) Q Y = - ^ Q ( i n ) Q ( i n ) Q Y + ^ Q ( f e n ) Q Y . Combining these three results according to Equation A.5 , the second partial derivative of the restricted log-likelihood function with respect 9 is a2 > | ? £ V , » | Y ) ^ ( I " ) Q ( £ " ) Q Y + 2 ^ ^ I W°) Q Y 1 + 2 t r 1 " 2 t r *»(»n)n-'(ln /a 2 N (A.6) Finally, the second partial derivative of the L*(-) with respect to both a2 and 9 can be obtained from partial differentiation of Equation A.3 with respect to 9 d2 L*(a2,9\Y) = da289 89 2a' Y Q Y (A.7) Asymptotic Covariance Matrix The asymptotic covariance matrix of the R E M L estimators of a2 and 9 is given by the inverse of the information matrix J(a2,9) = -E £*L*(a2,9) ^ L * ( a 2 , 9 ) ^-9L*(a2,9) g,L*(a2,9) (Cressie and Lahiri 1996). The elements in the right-hand side of Equation A.8 are given by (A.8) - E ( £ t L * ( a 2 , 8 ) ) = - ^ + E ( ^ Q Y ) = - ^ + £ E ( t r [ Q Y Y T ] ) = - ^ + £ ( ^ t r [Qfi] + tr [QX/3(3TXT]) 2 a 4 65 -E (s&^V.«0) = E ( ^ Q (|ft) Q Y ) = 2 ^ E ( t r [ Q ( ^ f t ) Q Y Y r ] ) Q ($0) Q (<7 2ft + X / 3 0 T X T ) = oirtr 2<x 1 2? 1 •tr = 27>tr Q ( A « ) Q (A«) Q X / 3 / 3 T X 2 -E ( | U V 2 , )^) = E ( ^ Q (&n) Q (A") Q Y - ^ Y ^ Q ( & n ) Q Y ) tr + i t r Q(£n) Q(An)n^(An)] = r^ E (tr [Q (An) Q (An) Q Y Y R ] ) - ^ E (tr [Q ( * ,n ) Q Y Y ^ ] ) - ^ [Q (An) ft-1 (An)] + i t r [Q ( & n ) ] = £ t r [Q (An) Q (An) Q (*2n + x/3/3Txr)] - 2 > [Q ( f t n ) Q ( a 2 f t + X(3(3TXT)] - [Q (An) n - 1 (An)] + ±tr [Q = tr [Q (An) Q (An) Qn] - £tr [Q (J^ft) Qft] - Jtr [Q (An) ft-1 (An)] + i t r [Q ( J J f t ) ] = tr [Q (An) Q (An)] - i t r [Q (&n) ft-1 (An) Q(An)Q(An) = i t r The information matrix in Equation A.8 can therefore be written as 2pr 2 ^ z r V^S 961 J ( Q Hr (Q (A .9 ) 66 A p p e n d i x B On the Equivalence of OLS and GLS Consider the general linear model of the n x l response vector ( Y s ) for a Nelder experiment consisting of a single type l a spacing wheel Y s = X / 3 + £ s (B . l ) where S is a n x 2 matrix that defines the relative spatial locations of the trees, X is a n x p design matrix, (3 is a p x 1 vector of unknown parameters, E$ ~ N ( 0 n x i , V ) is a n x 1 vector of random errors and V is an n x n positive definite matrix. The relative spatial location of the i t h tree is completely determined by the radius of its arc (r^; k = 1, 2 , . . . , t) and the angle of its spoke [zuv, I = 1 , 2 , . . . , m), that is, rk TZ>1 Izu (B.2) where sf is the i th row of S (i = 1, 2 , . . . , n) and zu is the angle between successive spokes in the design. Theorem 1 If the random vector es for a Nelder experiment consisting of a single complete type la spacing wheel with m spokes and t experimental arcs is defined as an isotropic process, then the ordinary least squares (OLS) and generalized least squares (GLS) estimators of (3, (3ols = ( X T X ) X T Y S (B.3) and PgU = ( X T V ~ 1 X ) X ^ ^ Y g • (B.4) respectively, are equivalent, andf30is is the minimum-variance unbiased estimator of/3. 67 Zyskind (1967) outlined the general conditions under which OLS and G L S estimators of (3 are equivalent. In particular, Zyskind demonstrated that a necessary and sufficient condition for the equivalence of the OLS and G L S estimators is that V can be decomposed as V X = X D (B.5) where X is the n x p design matrix and D is any p x p matrix. It is not required that X be of full column rank. The following two lemmas will be used to show that the identity defined by Equation B.5 applies under the conditions of Theorem 1. Lemma 1: If es is an isotropic stochastic process, then the elements of V = {vij}, arranged lexicographically by spoke and arc, are constrained by the relation E v(h+(a-l)t),j = E Vh,j 3=1 3=1 for all a — 1, 2 , . . . , m; and h = 1, 2 , . . . , t (i.e., (h + (a — l)t) = 1 ,2 , . . . , n). Proof: (B-6) u(fc+(o-i)*)j = E f ( 3=1 S ( M - ( o - l ) t ) - s3 3=1 t ra ( E E / k=l 1=1 t m ZUl = E E f { rh + r l ~ 1rhrk cos (azu - Ivufj fc=l l=i t m E E / {rh + T l - 2rhrk COS (ZU - C U 7 ) ) \c = 1 + I - al fc=l c = l t m = E E / k=lc=l ZUl rk, ZUr = £/ (K- s i l l ) 3=1 n = zZVhJ 3=1 where the first equality follows from the condition of isotropy. Lemma 2: If es is an isotropic stochastic process, then the elements of V = {vij}, arranged lexicographically by spoke and arc, are constrained by the relation E V(h+(a-l)t),(h'+(l-l)t) = E vh,(h'+(l-l)t) 1=1 1=1 (B.7) 68 for all a = 1, 2 , . . . , m; h = 1, 2 , . . . , t; and h! = 1, 2 , . . . , t. Proof: 1=1 rh rw Y V(h+(a-l)t),(h'+(l-l)t) - Ylf{ S(h+(a-l)t) ~ S ( V + ( / - l ) t ) ) ' " 1=1 m = Yf i=i m = Yf (rh+ rh' ~ ^rhrh, cos (azu - Izu)} 1=1 m = Yf(rh + r\> - 2rhrh, cos (zu - czu)} c=\ m = Yf [c = 1 + I - a] c = l ZUl rh' zuc III = Yf { S " ~ S ( A ' + ( c - l ) t ) ) c = l m - Y vh,(h'+{c-i)t) c=l from which Equation B.7 follows. It is now possible to show that the conditions of Theorem 1 imply the identity given by Equation B.5. Since X need not be of full column rank, it is possible to write V X = V [ [ l t x l I t ] ® l m x l ] where ® is the direct (Kronecker) product operator. Therefore, (B.8) V X = =1 VU 2-^r= =1 v2,j 2^r= £"= -1 vn,j 2_JT= ai Ki •• b •• b J2r=l vl,t+(r-l)t T,r=l v2,t+(r-l)t Em r=l un,t+(r-l)t (B.9) i,t %t (B.10) an bn,l • • • bnt 69 From lemmas 1 and 2 , respectively, a h = a h + { l _ l ) t and &(fc+(i-i)t),(fe'+(J-i)t) = bh>(h,+{l_1)t) (h = 1 , 2 , . . . , t;ti = 1 , 2 , . . . , t; I = 1 , 2 , . . . , m). Thus, V X = 0-2 &2,1 at h,i 1(txl) It • • • £>2,t ® 1 ' r raxl) ( B . 1 1 ) • • • h,t 0 0 0 a\ 61,1 02 h,i ^2,t ® l ( m x l ) ( B . 1 2 ) at hti 0 0 ... o ai ••• 6i ,t ® l ( m x l ) 02 &2,1 • • • h,t ( B . 1 3 ) at ^,1 • • • h,t = X D ( B . 1 4 ) Since V X = X D is a necessary and sufficient condition for the equivalence of the OLS and G L S estimators of (3 (Zyskind 1 9 6 7 ) , Theorem 1 is proved. The two lemmas define a type of symmetry that is implied by the condition of isotropy assumed in Theorem 1 ; the elements of the rows of V corresponding to same arc, but different spokes, are simply permutations of one another. This symmetry in V is a consequence of the symmetry in the treatment structure in the type la design. This symmetry, and the equivalence of the OLS and G L S estimators of /3, is lost if trees have died or observations are missing from the analysis. Theorem 1 can also be generalized to Nelder experiments consisting of several type la spacing wheels, or combinations of several type la and type lb spacing wheels, if differences among spacing wheels are modelled as fixed effects. 7 0
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A comparative study of spatial analysis methods for...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A comparative study of spatial analysis methods for forestry Nelder experiments Affleck, David L. R. 2001
pdf
Page Metadata
Item Metadata
Title | A comparative study of spatial analysis methods for forestry Nelder experiments |
Creator |
Affleck, David L. R. |
Date Issued | 2001 |
Description | The Nelder (1962) series of systematic spacing wheel designs define compact and spatially explicit layouts for forestry spacing experiments. One difficulty in the analysis of Nelder experiments is that the compact arrangements of trees may result in significant correlations among neighbouring variable values when these are geographically ordered. When this is the case, classical analysis methods, based on random sampling models and ordinary least squares estimators, provide inefficient estimates of treatment effects and biased estimates of variance parameters. However, the formal geometric structure of Nelder experiments presents an opportunity to evaluate the properties of alternative analysis methods based on spatial correlated error models. The statistical and practical properties of four spatial analysis methods were evaluated in the context of the interpretation of data from forestry Nelder (1962) experiments. The spatial analysis methods considered were based on either regionalized variables or nearest neighbour models and restricted maximum likelihood estimators. The validity and relative efficiency of the spatial analysis methods were assessed under different magnitudes of spatial autocorrelation through a simulation study. The same spatial analysis methods were also applied to data from a Nelder experiment in the University of British Columbia Malcolm Knapp Research Forest in Haney, British Columbia. The spatial analysis methods based on regionalized variables models were valid and more efficient than the classical analysis method when data were strongly spatially correlated. The spatial analysis methods based on nearest neighbour models were not valid and were less efficient than the classical analysis method. Application of the spatial analysis methods required more intensive estimation procedures. |
Extent | 5249598 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-08-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0075393 |
URI | http://hdl.handle.net/2429/11615 |
Degree |
Master of Science - MSc |
Program |
Forestry |
Affiliation |
Forestry, Faculty of |
Degree Grantor | University of British Columbia |
GraduationDate | 2001-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2001-0303.pdf [ 5.01MB ]
- Metadata
- JSON: 831-1.0075393.json
- JSON-LD: 831-1.0075393-ld.json
- RDF/XML (Pretty): 831-1.0075393-rdf.xml
- RDF/JSON: 831-1.0075393-rdf.json
- Turtle: 831-1.0075393-turtle.txt
- N-Triples: 831-1.0075393-rdf-ntriples.txt
- Original Record: 831-1.0075393-source.json
- Full Text
- 831-1.0075393-fulltext.txt
- Citation
- 831-1.0075393.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0075393/manifest