Modelling Spatial Dependence in an Irregular Natural Forest

The spatial dependence present in a natural stand of Eucalyptus pilularis (Smith) dominated mixed species forest was characterised and modelled. Two wildfires imposed a significant spatial dependence on the post disturbance stand. It was hypothesised that spatial variation in the intensity of the wildfires generated the observed structures. The influence of patch formation, micro-site variability and competitive influences were also noted in the residuals of a distance-dependent individual-tree growth model. A methodology capable of modelling these complicated patterns of observed dependence was sought, and candidates included the spatial interaction, direct specification and Papadakis methods. The spatial interaction method with a moving average autoregression was identified as the most appropriate method for explicitly modelling spatial dependence. Both the direct specification and Papadakis methods failed to capture the influence of competition. This study highlights the possibility that stand disturbances such as natural and artificial fires, insect and fungal attacks, and wind and snow damage are capable of imposing powerful spatial dependencies on the post disturbance stand. These dependencies need to be considered if individual tree growth models are to provide valid predictions in disturbed stands.


Introduction
There is growing agreement among forest growth modellers that the individual trees of forests can no longer be considered independent units (Tomppo 1986, Fox et al. 2001, Garcia 2006).Independence has traditionally been assumed to facilitate the use of classical statistical methods such as ordinary least squares.Several recent studies in even-aged stands (Fox et al. 2007a, Fajardo andMcIntire 2007) have indicated that individual trees are spatially dependent, and in Fox et al. (2007b), this spatial dependence was incorporated in an individual tree growth model.In this study we depart from the relatively simple spatial scenario of evenly-aged and evenly-spaced stands to an irregular naturally occurring mixed species stand dominated by Eucalytpus pilularis.
Irregularly structured E. pilularis is charactersised by complex spatial patterns (Florence 1996) and in this instance patterns are further confounded by spatially heterogeneous disturbance due to wildfire and harvesting.The challenge is to characterise and model these complicated spatial structures so they can be incorporated in an individual-tree growth model.
The study of spatial dependence in forests has diverged into two distinct literatures in the last 20 years; studies that examine spatial pattern in even-aged stands (Samra et al. 1990, Magnussen 1994, Fox et al. 2007a), and studies in irregularly structured stands (Sakai and Oden 1983, Biondi et al. 1994, Zhang et al. 2004).The two scenarios are quite distinct; spatial patterns in even-aged stands are predictable through stand development (Fox et al. 2001, Fajardo andMcIntire 2007), whilst those in irregular stands are rarely predictable and are subject to multiple confounded spatial influences (Kenkel et al. 1997).This contribution examines the later scenario, and intends to advance this literature by explicitly modelling complicated patterns of spatial dependence in an individualtree growth model.
Characterisation of the spatial dependance affecting forest growth and yield models is rarely done (Dennis et al. 1985, Gregoire et al. 1995), even though it may be having a significant influence on the estimation and prediction properties of individual-tree growth models (Garcia 1992, Fox et al. 2001).In this study residuals from a distance-dependent, individual-tree model developed for irregular E. pilularis are examined for the presence of spatial dependence.An observed significant residual spatial dependence would make it desirable to model the dependence for inclusion in growth models.Possible methodologies include the spatial interaction methodology (Cressie and Hartfield 1996), the direct specification methodology (Zimmerman and Harville 1991), and the Papadakis method (Bartlett 1978).These methods were applied to even-aged forests in Fox et al. (2007b), but have not been applied to irregularly structured forests.
This study examines spatial growth dynamics over 50 years of measurement of the Burrawan native forest growth plot.The Burrawan plot is historically significant, representing the first effort to collect growth information on the coastal eucalypt resource of New South Wales, Australia.Coastal E. pilularis does not form reliable annual growth rings, therefore long-term measurements of individual tree growth were required to estimate volume tables.In 1948, the Burrawan plot was measured for individual tree diameter and height.In 1951, measurements from the Burrawan plot were used to compile one of the first scientifically credible volume tables for the New South Wales coastal eucalypt resource, and an optimal rotation age of 80 years was identified (Henry 1955).The Burrawan growth plot was measured 16 times between 1948 and 1998, representing a rare insight into the growth and spatial dynamics of a naturally occurring stand.There is a paucity of long-term, spatially-mapped growth experiments (Biondi et al. 1994), and the Burrawan growth plot is a rare exception.This contribution represents the first published account of the spatial growth dynamics in this historical plot.
Studies in irregular forests have predominantly observed positive spatial dependence among individual-tree attributes (e.g.Sakai and Oden 1983, Geburek and Tripp-Knowles 1994, Biondi et al. 1994, Kuuluvainen et al. 1996, Kenkel et al. 1997, Kuuluvainen et al. 1998, Zhang et al. 2004).Positive dependence has been associated with structural irregularities such as patches of similarly sized trees (e.g.Sakai and Oden 1983, Geburek and Tripp-Knowles 1994, Biondi et al. 1994, Kenkel et al. 1997).In addition to struc-tural irregularities, micro-site variability may be amplifying positive spatial dependence among individual trees (Matern 1960).Spatial variation is further complicated by competition that induces a confounding negative dependence over short inter-tree distances (Matern 1960).This confounding influence has been observed in a mature mixed scots pine (Pinus sylvestris L.) and silver birch (Betula pendula Roth) stand in Northern Finland by Penttinen et al. (1992).Further, the existence of a persistent sub-canopy in an irregular forest may strengthen negative spatial dependence (Kenkel et al. 1997).The presence of multiple spatial influences can generate complicated patterns of spatial dependence in natural forests.If these complicated patterns can be characterised and modelled they can be incorporated in forest growth models for improved estimation and prediction (Fox et al. 2001).
Previous studies in natural forests have quantified spatial dependence for a number of individual tree attributes including diameter (Sakai and Oden 1983, Penttinen et al. 1992, Geburek and Tripp-Knowles 1994, Biondi et al. 1994, Kuuluvainen et al. 1996, Kenkel et al. 1997, Kuuluvainen et al. 1998), diameter increment (e.g.Biondi et al. 1994), height (e.g.Penttinen et al. 1992, Geburek and Tripp-Knowles 1994, Kuuluvainen et al. 1996, Kuuluvainen et al. 1998) and diameter-height relationships (Zhang et al. 2004).Nearly all these studies observed positive spatial dependence in the attributes of individual trees.Biondi et al. (1994) observed consistent positive spatial dependence in diameter through 70 years of measurement for an old growth forest, which they attributed to a constant patch size.Biondi et al. (1994) also observed that spatial dependence in diameter increment was more variable across the 70 years of measurement.Spatial dependence in diameter increment is more variable because it is more sensitive to seasonal climatic influences and the confounding influence of competition.
Competition, micro-site influences, and structural irregularities such as patches have been recognised as key factors that influence the spatial structure of natural forests.However, they are not the only factors affecting spatial structure, and other less recognised factors include genetic architecture and a variety of natural and artificial disturbance events.Spatial genetic architecture, which is only prevalent in natural forests, has been found to produce negligible spatial dependence with respect to the presence of particular genetic characteristics (e.g.Epperson and Allard 1989, Xie and Knowles 1991, Knowles 1991, Geburek 1993, and Leonardi et al. 1996).It is therefore unlikely that it will exert a significant influence on spatial patterns of tree size and growth in natural forests (Sakai and Oden 1983).
Disturbance events affecting forests include artificial and natural fire, insect attacks, fungal attacks, and the impacts of environmental factors such as wind and snow damage.Most of these disturbances are spatially heterogeneous with respect to their intensity.When the intensity of such disturbances is spatially structured it could be hypothesised that the size and growth of trees in the post-disturbance forest will follow this structure (Miller and Urban 1999).Disturbance from fire is an example of an event that is strongly spatially heterogeneous; the intensity of a fire at any point in the stand is dependent on topography, fuel load, vegetation mosaic, moisture, and wind gusts (Williams et al. 1994).Disturbances from insect attack have also been shown to be spatially heterogeneous (Preisler 1993, Preisler et al. 1997).Disturbances from wind, snow and some pathogens would be expected to exhibit similar heterogeneity over space.Disturbance events should impose strong spatial dependencies on the post disturbance stand suggesting that models explicitly incorporating spatial dependence may be required in disturbed stands.The Burrawan plot is subject to two wildfire events and two harvesting events, and examination of post disturbance spatial patterns can test the hypothesised strengthening of spatial dependence, and possibilities for its explicit modelling.
The objectives of this study were to characterise and model the spatial dependence among individual trees in an irregular E. pilularis forest, and examine whether spatially heterogeneous disturbance events induce spatial dependence in the post disturbance stand.

The Burrawan Native Forest Growth Plot
The ), and Turpentine (Syncarpia glomulifera (Smith) Niedenzu).All species were included in the following analysis, and complex spatial structures may evolve in these forests because of this mix of species (Florence 1996).The growth plot has a long measurement history commencing in 1948, and was originally established to satisfy a need for growth data in coastal forests (Henry 1955).
Between 1948 and 1998 the Burrawan plot was measured for individual tree diameters 16 times.
The spatial positions of individual stems were also recorded in 1948, 1960, and 1983.In 1998 the diameter, height, and spatial positions of all remaining stems and stumps were measured by the primary author.
The Burrawan plot was subject to two important disturbance events.The plot was severely burnt in 1952 by a wildfire that affected a large area of coastal eucalypt forest around Wauchope.Some areas of the plot were severely affected and trees from these areas were subsequently removed by a salvage operation in 1960.Wildfire again affected the plot in 1978.Although this fire was less intense, in some areas of the plot trees were scorched to a height of 20 metres.These two wildfires would be expected to exhibit spatially heterogeneous intensity across the large Burrawan plot.We could therefore expect an expression of this spatial heterogeneity in individual tree growth in the post disturbance stand.

Individual-Tree Growth Model
A distance-dependent, individual-tree growth model was developed for irregular E. pilularis (Fox 2000): where ADI i is the annual diameter increment, D i the diameter, APASWC i the squared weighted and constrained area potentially available index (Nance et al. 1987), and TRVRB i the change in the RB competition index as a result of a thinning event for tree i. RB is the distance-independent index of Glover and Hool (1979) and ε i is the residual term.
Model (Eq. 1) was fitted to each individual measurement of the Burrawan plot and a set of residuals estimated.Boundary trees (trees within 15 metres of the plot boundary) were excluded, thus preventing the errors associated with a boundary effect (Monserud and Ek 1974).Residuals were then subject to analysis for the presence of spatial dependence.

Characterising Spatial Dependence in Model Residuals
Spatial dependence in model residuals was characterised using correlograms based on standardised Moran's I (Moran 1950, Cliff andOrd 1981).Moran's I was converted into a standard normal deviate using the expected value and variance under the normality assumption, as defined by Cliff and Ord (1972).The application of correlograms is dependent on the assumption of spatial stationarity, i.e., directional and locational invariance.This was confirmed for measurements of the Burrawan plot using median polish and directional semivariograms (Cressie 1986).
For investigating residual dependence for the candidate methodologies, Moran's I was used to construct unstandardised correlograms.The candidate methods alter the variance-covariance matrix to incorporate spatial dependence.With off-diagonal terms in the matrix, standardised Moran's I is not longer valid (Anselin and Kelejian 1997), and instead an unstandardised correlogram is used.

Modelling Spatial Dependence
Several candidate methodologies exist including the spatial interaction method, the direct specification method, and the Papadakis method (Fox et al. 2007b).The intention of these methods is to ameliorate residual spatial dependence.Because Equation 1 already includes deterministic structures for tree size, competition and thinning, residuals will include adjustments to these structures.Caution needs to be practiced when interpreting these confounded residual spatial structures.A brief description of their application on the Burrawan plot is included below.

Spatial Interaction Models
Spatial interaction models employ probability distributions to model the spatial dependence prevalent among a set of points in space.The simultaneous autoregression (SAR) uses a joint probability distribution to model the autoregressive process in the response (Whittle 1954) while the conditional autoregression (CAR) uses a conditional probability distribution to model the autoregressive process in the response (Besag 1974).The moving average autoregression (MA) uses a joint probability distribution to model the autoregressive process in the error term (Haining 1978).
Implementation of the three autoregressive models requires a correctly designed spatial dependence matrix.A design capable of simultaneously modelling spatial dependence induced by competitive and micro-site influences (Fox et al. 2007b) was adapted for the Burrawan plot.This design used two matrices and facilitated a separate spatial dependence parameter for the two influences.Each element (s ij ) of the final spatial dependence matrix S for the SAR specification was generated using the sum of elements from the two matrices WC and WM: where θ c and θ m are the spatial dependence parameters and wc ij and wm ij are elements of the spatial weights matrices characterising competitive and micro-site influences respectively.Similar definitions of the spatial dependence matrix for the MA and CAR specifications were also used.
Residual correlograms were examined to identify the scale of spatial dependence, and each element of the competitive matrix (wc ij ) and the micro-site matrix (wm ij ) for the various measurements of the Burrawan plot were calculated as in Table 1.The 1952 measurement did not exhibit any negative dependence over short inter-tree distances, therefore the competitive matrix was omitted.Residual correlograms for the 1952 and 1978 measurements showed a slow decline in positive spatial dependence for larger inter-tree distances.To effectively model this decline, the distance over which trees were included was expanded to 20 metres, and the inverse square root of inter-tree distance was used as the weight, thus providing a slower decline in wm ij for larger inter-tree distances.
Estimation of spatial interaction models was carried out by maximum likelihood methods (Mardia 1990, Cressie 1993).The S-plus function slm was used to estimate the spatial interaction models (Kaluzny et al. 1997).
h ij is the physical distance in metres separating trees i and j

Direct Specification Methods
Direct specification methods use a mathematical function to directly specify the variance-covariance matrix (Zimmerman and Harville 1991) thus spatial dependence is explicitly incorporated.Possible functions are restricted in that they must preserve the symmetry and positive definite form of the variance-covariance matrix (Ripley 1981).
Permissible functions applied to the Burrawan plot included the spherical (Eq.3), exponential (Eq.4), and Gaussian (Eq.5) models: where C(h) is the covariance at distance h, σ n is the nugget variance which captures random spatial variation, and ρ is the spatial dependence parameter.Direct specification models are estimated using maximum likelihood methods (Mardia andMarshall 1984, Zimmerman andHarville 1991), and in this instance The MIXED (SAS Institute Inc. 1996) procedure in SAS was used to estimate parameters.

Papadakis Method
The Papadakis method (Papadakis 1937(Papadakis , 1984) ) is perhaps the simplest means of incorporating spatial dependence in model estimation.The Papadakis method uses the average model residual for trees in spatial proximity as a covariate in the model, and has been confirmed as being approximately valid (Bartlett 1978, Taye andNjuho 2007).Residual correlograms were examined to identify the scale of spatial dependence, and based on this the average residual was generated for all trees less than 10 metres from the subject.For the growth periods during which fire disturbances occurred (1952,1978), the distance over which residuals were averaged was expanded to 25 metres.This was done in response to the larger scale of spatial dependence observed for these measurements.
When applied to the earlier developed distance-dependent, individual-tree model (Eq. 1) the Papadakis method can be described as: where p i is the Papadakis covariate, N(j) denotes the number of neighbours of the subject i, ε j is the residual from the fitted regression (excluding term γ p i ) for tree j, which is a spatial neighbour of tree i.

Characterising Spatial Dependence
Spatial correlograms for diameter for all measurements of the Burrawan plot are provided in Fig. 1.Patterns of spatial dependence for individual tree diameter are relatively consistent across the 15 measurement periods of the Burrawan plot.A significant positive spatial dependence is observed for inter-tree distances less than 4 metres.This initially significant dependence declines to zero for increasing inter-tree distance.Structured spatial variation to inter-tree distances of 40 metres is an indication that patches of similarly sized trees are present on the Burrawan plot.The salvage logging of 1960 altered patterns of spatial dependence for diameter.Most noticeably the positive spatial dependence prevalent over short inter-tree distances became more significant.
Spatial correlograms for diameter increment are provided in Fig. 2. Patterns of spatial dependence for diameter increment are far more erratic across the 15 measurement periods than those for diameter.Early measurements (1948, and 1950) exhibit a negative spatial dependence for intertree distances less than 4 metres that may reflect competition between trees in spatial proximity.Later measurements exhibit a significant positive spatial dependence for the smallest distance class.The influence of the salvage logging in 1960 can be observed with altered patterns of spatial dependence after 1960.The salvage logging of 1960 may have reduced competitive pressures with positive dependence emerging over short inter-tree distances after this event.The emerging positive dependence may also be caused by the strengthened positive spatial dependence in diameter after 1960, as can be observed in Fig. 1.The influence of a thinning event in 1990 can be noted with strong positive dependence over small inter-tree distances for 1990.Spatial correlograms of standardised residual Moran's I against inter-tree distance for residuals from model (Eq. 1) are provided in Fig. 3.The influence of wildfires in 1952 and 1979 can be observed, with strong patterns of positive spatial dependence for the 1952 and 1978 measurements.For several measurements there is a significant negative spatial dependence for inter-tree distances less than 4 metres, possibly an expression of competition between trees in close spatial proximity.Several measurements exhibit a significant positive spatial dependence for small inter-tree distances.The presence of significant spatial dependence in residuals indicates a need for explicit models of spatial dependence.

Modelling Spatial Dependence
Possible spatial interaction models included the simultaneous, moving average, and conditional autoregressions.The distance-dependent, individual-tree model (Eq. 1) constituted the deterministic structure.Examination of likelihood ratio tests (LRT) for the significance of spatial dependence parameters and Akaike's information criterion (AIC) identified the moving average autoregression as best.Spatial dependence parameters for the moving average autoregression are shown in Table 2.The spatial dependence parameters     2 show that θ c was either negative or almost zero.This is consistent with this component modelling dependence induced by competition, although no estimate was significantly different from zero.θ m was positive for all but one of the measurement periods, which indicates that this component modelled dependence induced by factors other than competition such as the presence of patches, micro-site influences or wild fire disturbance.
Possible direct specification models include the spherical, exponential, and Gaussian models.A graphical comparison of each model against the empirical semivariogram, as well as examination of Akaike's information criterion identified the spherical model as best.
The Papadakis method is the third possible methodology for modelling observed spatial dependence on the Burrawan plot.Table 3 provides information on the change in model fit when the Papadakis covariate is included in the model.Information on the parameter for the Papadakis covariate as well as its significance are also provided.From Table 3 it can be noted that the Papadakis method is providing significant improvements in model fit for several measurements of the Burrawan plot.In particular, the 1952 and 1978 measurements, which were affected by wildfires, exhibit marked improvement in model R-square, and reductions in mean squared error of 25% and 37% respectively.The parameter for the Papadakis covariate is predominantly positive, indicating that positive spatial dependence is being modelled by the covariate.
Residual correlograms for the original deterministic model (Eq. 1) and the optimal model of each candidate methodology (moving average, Papadakis, and spherical model)

Discussion
Individual tree diameters on the Burrawan plot exhibited consistently structured spatial variation over 50 years of measurement.A significant positive spatial dependence was observed for small inter-tree distances that declined and became negligible for trees separated by more than 40 metres.This structure was indicative of the presence of patches of similarly sized trees (Sakai and Oden 1983, Geburek and Tripp-Knowles 1994, Biondi et al. 1994, Kenkel et al. 1997).Because spatial dependence extended to inter-tree distances of 40 metres it could be inferred that the average patch size on the Burrawan plot is approximately 40 metres in diameter (Biondi et al. 1994, Geburek andTripp-Knowles 1994).This average patch size is more likely a composite of several irregular patches and should not be taken as an indication that patch size is regular or consistent across the Burrawan plot.Other studies in natural forests have observed similarly structured spatial variation in tree diameter (e.g.Sakai and Oden 1983, Biondi et al. 1994, and Geburek and Tripp-Knowles 1994 observed positive spatial dependence to 50, 30, and 20 metres respectively).Salvage logging in 1960 appeared to strengthen positive spatial dependence in diameter.This is a likely result of the removal of suppressed trees or the smaller tolerant species in the logging of 1960.Removal of suppressed trees should extinguish any negative dependence prevalent over short inter-tree distances, and should produce patches of similarly sized trees structured at the scale of micro-site influences.The influence of the 1960 logging was also expressed in diameter increments with significant positive dependence emerging after the event.
Spatial variation in diameter increment was far more erratic across the 50 years of measurement of the Burrawan plot.Several measurements exhibited negative spatial dependence over short inter-tree distances, which is a likely consequence of competition (Penttinen et al. 1992) and the presence of a persistent sub-canopy of tolerant species (Kenkel et al. 1989).For other measurements a significant positive spatial dependence existed over short inter-tree distances, which is an expression of micro-site influences (Matern 1960) and patches of similarly sized trees (Biondi et al. 1994).For several measurements a complicated pattern of spatial dependence emerged because processes inducing positive spatial dependence were confounded by processes inducing negative spatial dependence.This complicated pattern was characterised by increasing positive spatial dependence over small inter-tree distances which reached a maximum in the vicinity of 10 metres before declining for larger inter-tree distances.Those measurements affected by disturbance events, in this instance wildfire, exhibited a strong positive spatial dependence to inter-tree distances up to 80 metres.This observation confirms that a spatially heterogeneous disturbance event such as a wildfire is capable of inducing a strong spatial dependence in the growth of the post disturbance stand.
Characterisation of spatial dependence for growth model residuals confirmed that significant dependence were present in the residuals.Spatial variation in model residuals was similar to that exhibited by diameter increment, indicating that the deterministic model structure is failing to effectively model the spatial variation present on the Burrawan plot.This is occurring despite the inclusion of a spatially-explicit competition index.
All the methods for modelling spatial depend-ence effectively modelled positive spatial dependence, however, only the moving average autoregression was capable of modelling negative spatial dependence.This capability was introduced by manipulating the spatial dependence matrix into two components; one component modelled the positive dependence and the other component modelled the negative dependence.The moving average autoregression was the most appropriate methodology for modelling the complicated patterns of spatial dependence observed on the Burrawan plot.This study confirms the applicability of the moving average autoregression to natural forests of irregular structure.For the Burrawan plot the competitive component (θ c ) was not significantly different to zero.Despite it being not significant, the competitive component prevented the expression of negative dependence in model residuals.Studies in more competition tolerant species, such as the temperate hardwoods of the United States, are required to identify the wider applicability of the identified model.Disturbance events on the Burrawan plot induced a significant spatial dependence among the growth rates of individual trees.Two wildfires in 1952 and 1979 generated very similar patterns of spatial dependence in the post disturbance stand.A strong positive spatial dependence was observed to inter-tree distances of 80 metres.This indicates that individual trees responded to the disturbance more similarly when in spatial proximity.This is an expression of strong spatial heterogeneity in the intensity of the disturbance across the Burrawan plot.Although resultant patterns of spatial dependence were similar, the actual spatial heterogeneity of fire intensity may have been quite different for the two fires.The structural complexity of the mixed species Burrawan plot may have exaggerated the spatial heterogeneity of the two disturbance events.
Disturbance events continually affect forests.Examples include artificial and natural fires, insect and fungal attack, and wind and snow damage.The findings of this study alert us to the possibility that disturbance events impose strong spatial dependence on the stand.This dependence can threaten the validity of individual-tree growth model estimation for disturbed stands.Examination of other disturbance events is warranted to expose the generality of this result.

Fig. 1 .
Fig.1.Spatial correlograms for diameter for the Burrawan plot.The upper two-sided 95% confidence limit is also included.

Fig. 3 .
Fig. 3. Spatial correlograms for growth model residuals for the Burrawan plot.The two-sided 95% confidence limits are also included.

Fig. 4 .
Fig. 4. Spatial correlograms for residuals from the original model (Eq. 1) and the model including a moving average autoregressive (MA) component for the 1966-1990 measurements of the Burrawan plot.

Fig. 5 .
Fig. 5. Spatial correlograms for residuals from the model including the Papadakis covariate (PAP), and the model including a variance-covariance matrix based on the spherical function (SPH) for the 1966-1990 measurements of the Burrawan plot.

Table 1 .
Definition of the competitive and micro-site matrices.
The wildfires of 1952 and 1978 (both fires

Table 2 .
Spatial dependence parameters for the moving average autoregression.

Table 3 .
Assessment of the Papadakis method.
* F-test significant at the 0.05 level; ** F-test significant at the 0.0001 level PRMSE is the Percentage Reduction in Mean Squared Error displayed in Table