Modeling the spatial structure of topical forests

Forest Ecology and Management 110 (1998) 293±314

Modeling the spatial structure of topical forests JoaÄo L.F. Batistaa,*, Douglas A. Maguireb a

University of SaÄo Paulo, Luiz de Queiroz College of Agriculture, Department of Forest Sciences, CxP 9, 13400-970 Piracicaba - SP, Brazil b Oregon State University, College of Forestry Department of Forest Resources, 209 Peavy Hall, Corvallis OR 97331-5703, USA Received 8 July 1997; accepted 24 March 1998

Abstract The spatial structure of tropical forest stands under different management conditions was modeled as a series of different spatial point processes. Spatial patterns were ®rst assessed by K-function analyses to help choose a point process appropriate for observed patterns. The homogenous Neyman±Scott process accurately described live tree distribution in clear cut areas, where tree patterns tended to be aggregated. Parameters were estimated by minimizing Diggle's modi®ed least squares criterion, and goodness-of-®t was assessed by comparison to con®dence envelopes constructed by Monte Carlo simulation. Parameter estimates can be interpreted to help understand the ecological processes in¯uencing re-colonization of disturbed areas. The inhomogeneous Poisson process was investigated for simulating the spatial pattern of ingrowth trees in lower canopy strata. The intensity function of this process was inversely proportional to variables representing canopy density. As assessed by Monte Carlo generation of con®dence envelopes, the inhomogeneous Poisson process successfully portrayed the in¯uence of canopy structure on understory plant distribution in most stands. Tree mortality was modeled as a thinning process in which the probability of individual tree mortality was conditional on subject tree attributes and competitive environment. The thinning function took the form of a generalized linear model with a binomial error distribution and logit link function. In most stands, tree neighborhood variables were powerful predictors of mortality, but they were not important predictors in all plots. This suggests that the surrounding forest structure of a subject tree has considerable in¯uence on its morality, but competition is not the sole cause of tree morality in tropical forests. # 1998 Elsevier Science B.V. Keywords: Forest management; Spatial pattern modeling; K-function; Point process model; Neyman±Scott process; Inhomogeneous Poisson process

1. Introduction Traditional research on plant spatial pattern has emphasized testing for departures from complete spatial randomness (CSR). In tropical forests, spatial pattern analyses have been used to test the Janzen±

*Corresponding author. Tel: +55 19 4336155; fax: +55 19 4330681; e-mail: [email protected] 0378-1127/98/$19.00 # 1998 Elsevier Science B.V. All rights reserved. PII S0378-1127(98)00296-5

Connell theory, one implication of which is that the young trees should have more aggregated spatial pattern than mature trees. In this theory, predation is the driving mechanism that reduces plant population density and allows the existence of highly diverse plant communities (Janzen, 1970; Connell, 1970). The results of studies on the spatial pattern of several tree and shrub species seems to indicated that the Janzen±Connell theory might explain the spatial pattern of some species but cannot be taken as a rule for

294

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

all tropical forest woody plants (Hubbell, 1979, 1980; Augspurger, 1983; Coley, 1983; Clark and Clark, 1984; Sterner et al., 1986; Wong et al., 1990; Awasthi, 1990; Gibson and Brown, 1991; Bariteau, 1992; Condit et al., 1992). Tests of the CSR hypothesis have also been applied in studies of competition effects on dynamics of evenaged homogeneous forests (West, 1984; Franklin et al., 1985; Kenkel, 1988; Kenkel et al., 1989), mixed species temperate forests (Collins and Klahr, 1991; Szwagrzyk and Czerwczak, 1993), semi-arid shrubdominated stands (Welden et al., 1990), and savanna vegetation (Skarpe, 1991). In general, the competition is assumed to cause a systematic change in the spatial pattern from random or aggregated to a regular pattern as a result of density-dependent mortality. Model of forest dynamics that incorporate spatial interactions among trees and other forest vegetation have been used as an alternative approach to testing hypotheses about the in¯uence of forest structure on spatial pattern. This mechanistic approach has been applied to describe the spatial structure of tropical and temperate forests (Maguire et al., 1993) and as a component of models simulating forest dynamics (e.g., Busing, 1991). The objective of this paper is to identify stochastic point processes that can simulate the non-random spatial patterns of trees in tropical Atlantic forests in Brazil. The point processes were selected to provide realistic reproduction of the observed patterns under a simple set of biologically plausible assumptions. Three stand development processes were studied: 1. Tree establishment after clear cut: Spatial pattern is usually aggregated in clear cut areas during the re-colonization stage. This pattern is usually generated by environmental heterogeneity, spatial variation in the soil seed bank, non-random seed dispersal and predation, and a combination of these mechanisms. As the re-colonization process proceeds, spatial pattern becomes less aggregated and shifts toward a random or regular pattern. This pattern was simulated with a homogeneous Neyamn±Scott process to test if a simple clustering model could adequately ®t the observed pattern and to model the variation in spatial pattern among different clear cut areas through time.

2. Tree establishment under forest canopy: Variation in stand density and canopy structure results in a spatial mosaic on the scale of micro-sites that are very influential on regeneration and tree establishment. Light conditions are a major determinant of the undercanopy environment but tree establishment is affected by several other factors as well. To mimic this process, the spatial pattern in ingrowth trees was simulated as an inhomogeneous Poisson process (IPP). 3. Tree mortality: Tree mortality to a large extent is a direct response to the competitive environment created by the surrounding forest vegetation. In general, prediction of mortality is based on competitor and subject tree attributes, such as age, size, growth rate and crown size. The influence of tree neighborhood is here represented by point-specific basal areas, and mortality is modeled as a function of tree attributes and spatial variables. The spatial pattern of mortality is simulated as a `thinning' process on the existing live trees. 2. Material and methods 2.1. Field site The study area is located in the Atlantic Forest of eastern Brazil. The Atlantic Forest is a vegetation formation that occurs along the eastern Brazilian Costs and is coincident with the `Serra do mar'which produces a climate with high precipitation and almost no dry season. The study site is located at Linhares Reserve, owned by Cia. Vale do Rio Doce (CVRD), a state mining company, 30 km north of the town of Linhares, state of EspõÂrito Santo, on BR-101 Highway (Fig. 1). Its geographical location is between the latitudes 198060 S and 198190 S and longitudes 398450 W and 408190 W. (Jesus, 1987). The reserve has 21 787 ha and represents 25% of the remaining natural vegetation of the state, which is presently reduced to 1.5% of its original area (Jesus et al., 1992). The topography is dominated by ¯at-topped hill that are between 28 and 65 m above wide shallow valleys. The geological substrate is mainly composed of tertiary and quaternary sediments, and deep (2±4 m) podzolic soils of sandy texture and low fertility. The climate is considered `Awi' in the KoÈppen classi-

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

295

Fig. 1. Location and map of Linhares Reserve. Numbers in the map indicate the blocks of sustainable management experiment.

®cation, i.e., warm and humid with wet summers and dry winters. The air temperatures range from a minimum of 8.38C, with a mean of 238C, and precipitation ranges from 817 to 1639 mm (mean of 1230). Relative humidity remains high throughout the year, between 81.0% and 85.5% (mean 83.5) (Jesus et al., 1992). Study plots were located in a uniform dense forest called `Mata de Tabuleiro' that is the typical lowland

Atlantic forest of South Bahia and North EspõÂrito Santo, representing 63.1% of the reserve (Jesus, 1987). This forest is ¯oristically and physiognomically similar to the `terra-®rme' Amazon forest, but is distinct from other Atlantic forests covering the slopes of `Serra do Mar' because of the sparser understory and the abundance of palms (Rizzini, 1979). Although it is considered a rain forest (Veloso and GoÂes-Filho,

296

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

1986), it has low water surplus (47 mm) and can suffer a water de®cit of 66 mm between January and February (summer) (Jesus, 1987). Mean stand density is around 540 stems/ha (DBH>10 cm), but few large trees are present; on average, only 15 trees per hectare have DBH greater than 80 cm and the mean basal area is about 34 m2/ha (Jesus, 1987). At Linhares Reserve, 571 tree species are known, but only 298 tree species ocurried in the study plots (DBH>10 cm). Rosewood (Dalbergia nigra) is probably the most well known tree species, the only remaining populations of which are found in the Mata de Tabuleiro of Bahia and EspõÂrito Santo. 2.2. Field measurements Field data were obtained from a forest management experiment installed in 1980 with the objective of comparing the effect of different partial cuttings on forest growth, yield, and species diversity. The management alternatives are primarily concerned with wood production for lumber manufacture, and secondarily for ®rewood. The experiment was originally installed in a completely randomized block design with ®ve replications and nine treatments. With the exception of control and clear cut treatments, all silvicultural treatments are a type of `selective' management. In this analysis, only ®ve treatments are considered. 1. Control: Control treatment with no silvicultural intervention. 2. BA reduction: 15% of basal area was removed starting with the largest trees (by DBH). 3. Lower and upper cut: All trees with DBH<10 cm and DBH>80 cm were removed, as well as 30% of the remaining basal area, starting from the largest residual trees. 4. Clear cut: All trees were removed. 5. Selective cut: All trees with DBH>50 cm were removed, as well as 25% of the remaining basal area, starting from the largest residual trees. Lianas and climber were also cut in all treatments except in the control. The plots were installed in an area of `Mata de tabuleiro' that was considered homogeneous and free of previous human intervention. Initially, eight sites

were selected, and ®ve locations were randomly drawn among these. These locations were 75 m from the nearest road, to facilitate wood removal and to minimize the impact of roads on treatments response. Each 0.5 ha experimental plot was rectangular (100 m50 m) and plots were separated by buffer zones. Before treatment implementation, all trees with DBH>10 cm were identi®ed and measured for circumference at breast height. Circumferences were later transformed to DBH. Immediately after treatment in 1980, the plots were re-measured and mapped. Three additional remeasurements were made at 3-year intervals from 1983 to 1989. At each remeasurement ingrowth trees were identi®ed, measured and mapped, and dead trees were tallied. The plots were mapped by dividing them into 10 m10 m quadrats and taking the distance of each tree to two sides of the quadrat. 3. Analysis of spatial pattern Tree spatial patterns were investigated with the Kfunction proposed by Ripley (1977). The K-function, also called reduced second-order analysis, is a powerful descriptor of spatial patterns because it takes into account a wide range of scales, not only nearest neighbor distances (Cressie, 1991). This method has been applied in the analysis of tree spatial pattern in a wide variety of different forest types (Diggle, 1983; Cressie, 1991; Penttinen et al., 1992; Moeur, 1993; Batista, 1984; Haase, 1995). The K-function can be intuitively interpreted when applied to a tree map: if l is the number of trees per unit area, then lK(r) is the expected number of trees within the distance r from an arbitrary tree (Ripley, 1977). Under complete spatial randomness (CSR), the Kfunction has a simple but quadratic form K…r† ˆ r2 To linearize the function and stabilize the variance, the K-function is usually transformed to s ^ ^ ˆ K…r† ÿ r (1) L…r†  This transformation also facilitates graphical evaluation of departures from CSR because the random

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

patterns are represented by the horizontal line through ^ ^ zero in the plot of L…r† on r. Positive values of L…r† indicate aggregation and negative values indicate regularity. Edge correction is an essential component of estimating K-functions for trees mapped on ®eld plots. Most forest applications have used an isotropic edge correction suggested by Ripley (1977). In this study, we used a translation edge correction proposed by Osher and Stoyan (1981) because it is computationally simpler and allows a wider range of scales in the analysis:  2 X n n X A I…kxi ÿ xj k < r† ^ (2) K…r† ˆ n w…kxi ÿ xj k† iˆ1 j ˆ 1 i 6ˆ j where A is the size of the ®eld plot, n the number of trees in the tree map, kxi ÿ xj k the distance between tree i and tree j, I equal to 1 if kxi ÿ xj k < r and 0 otherwise and w…kxi ÿ xj k† is the translation edge correction, which for a rectangular region with sides a and b (a
297

with a mean equal to   jAj (jAj is the area of the region). Numerous stochastic point process models that depart from CSR have been developed and explored in the spatial statistics literature. In this study, different point process models are ®tted to observed point patterns representing `tree maps' on permanent plots. 4.1. Homogeneous Neyman±Scott process The spatial pattern of established trees in the clear cut plots was modeled as homogeneous Neyman± Scott process. This process generates clustered patterns and is de®ned by the following stochastic mechanism (Diggle, 1983; Cressie, 1991): 1. Parent events follow a HPP process where the `intensity' of parent events (or number of parent events per unit of area) is de®ned by the parameter . 2. The number of offspring per parent even follows a Poisson distribution. 3. The offspring events are located in relation to the parent events according to a bivariate Gaussian distribution with mean 0 and variance s2I. To ®t this model Diggle (1983) suggested the use of the expected value of K-function. It is frequently used to ®t stochastic point processes because its expected value can be derived from the stochastic mechanism that de®nes the process. For the homogeneous Neyman±Scott process, the expected value of the K-function is K…r; ; † ˆ r 2 ‡

1 ÿ exp…ÿr 2 =42 † ; r > 0 (3) 

where s and r are parameters that control the process (Cressie, 1991). To ®t the homogeneous Neyman±Scott process to an observed point pattern, Diggle (1983) suggested the `modi®ed least square criterion' which consists of minimizing the expression Zr0 D…; † ˆ

c ^ f‰K…r†Š ÿ ‰K…r; ^; ^†Šc g2 dr

(4)

0

^ where K…r† is the estimated K-function for the observed tree map (the estimator proposed by Osher

298

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

and Stoyan (1981) based on translation edge correction was used), K…r; ^; ^† is the expected K-function for the process given the parameter estimates ^ and ^, r0 and c are the `tuning constants' necessary to control the variance of the estimated K-function to assure more stable convergence in the minimization process. The starting values for the minimization process were obtained by ®tting Eq. (3) through nonlinear least squares, and using sˆ1 and rˆ1 as initial values in Eq. (4). In the case that the non-linear regression failed to converge, other arbitrary values were used to start the minimization of D(s,r). The modi®ed least squares criterion also requires choices of the tuning constants: c (power transformation) and r0 (limit of integration). For aggregated patterns, Diggle (1983) suggested using c<0.25, and argues that values of r0 become less relevant for appropriate values of c. Several values of the tuning constants were tested using one plot (plot 1) to verify the effect of the constants on the result of the minimization procedure. The values for the power transformation were 0.10, 0.15, 0.25 and 0.50, while the limits of integration were set at 10, 25 and 45 m. Approximate con®dence envelopes (sensu Guttorp, 1995) for the K-function based on the ®tted models were constructed by Monte Carlo simulation of the Neyman±Scott process with the obtained parameter estimates (Diggle, 1983). 4.2. Inhomogeneous Poisson process The inhomogeneous Poisson process (IPP) was applied to the spatial pattern of ingrowth trees. This process is an expansion of the homogeneous Poisson process in the sense that intensity of the process is not constant but a function l(x) that varies with spatial location (x). Diggle (1983) presents the following stochastic mechanisms as de®nition of this process: 1. The number of events in a region A of observation follows that Poisson distribution with mean Z …x† d…x† A

2. Given the number of events in a region A is equal to n, the n events form an independent random sample from the spatial distribution on A with probability density function proportional to l(x).

The intensity function l(x) determines the IPP and hence, ®tting this process is a matter of ®nding an appropriate intensity function. To model the spatial pattern of ingrowth trees, we de®ned the intensity as inversely proportional to pointspeci®c basal area estimates. Therefore, spatial location in the study plots that had a high surrounding basal area would yield a low probability of having an ingrowth tree, while locations with a low surrounding basal area would have a high probability of having an ingrowth tree. The point-speci®c basal are estimates were consistent with the well-established concept of horizontal point sampling in forestry (Husch et al., 1982). For any point x within the study plots, point basal area PBAF …x†, is de®ned as the basal estimated that would be obtained with a basal area factor of F. Basal area factors (F) of 1, 2 and 3 m2/ha were tested and to avoid bias due to edge effects, the plots were considered to be torus. Since in some plots (clear cut) point-speci®c basal area estimates were 0 m2/ha, a small constant was added yielding the following intensity function: 1 ; >0 (5) …x† ˆ PBAF …x† ‡ 0:5 where is a parameter that must be estimated from the data. PBAF …x† is a non-negative function, so the requirement that l(x) be non-negative is satis®ed. An estimate of was obtained by using the property that the expected number of events of an IPP in a region A is given by Z N…A† ˆ jAj …x† jdxj A

so for each plot could be estimated by n R ^ ˆ jAj A ‰1=PBAF …x† ‡ 0:5Š jdxj where n is the observed number of events and the integral was numerically evaluated using the Composite Simpson's rule (Johnson and Riess, 1982). The performance of the model was assessed by comparison of approximate con®dence envelopes (Guttorp, 1995) for the linear transformation of the K-function (L(r)) with empirical estimates. The approximate con®dence envelopes were constructed in two steps by Monte Carlo simulation of the IPP (Diggle, 1983):

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

1. Simulate a HPP over the region A with intensity M ˆ sup‰…x†Š: x2A

2. The event xi (iˆ1,2,. . .,n) generated in step 1 are retained with probability p…xi † ˆ

…xi † : M

The intensity function de®ned in terms of pointspeci®c basal area is completely dependent on the stand structure within each study plot. To ®nd the maximum value of the intensity function and to integrate it over the entire plot, a grid of points was used. Two different grid spacings were considered, 1 m1 m and 5 m5 m (Fig. 2). Although the overall pattern was visually similar under both resolutions, the 11 gird was more computationally expensive when several simulations of the IPP were required to con-

299

struct the con®dence envelopes; hence, the 5 m5 m grid was used. In general, the effect of basal area factor (F) on the estimation of the intensity function and on the con®dence envelopes was not very strong. An F of 1 tended to produce con®dence envelopes conforming more closely to CSR, while larger F's seemed to enhance the differences between area with high and low basal areas, causing the con®dence envelopes to vary more from plot to plot. Approximate con®dence envelopes were constructed with Fˆ3 m2/ha. 4.3. Thinned process Tree mortality was related to stand structure by representing the process as a `thinned stochastic process'. In a thinned process, some events of an initial point process are delated and the remaining events form a new process. The `thinning function' de®nes

Fig. 2. Basal area point estimates (PBA) for (A) 1 m1 m and (B) 5 m5 m grids in control plot 1, using a basal area factor (F) 3.

300

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

the stochastic mechanism by which events are deleted. In this study, the initial point pattern was the pattern observed at the beginning of the study period (1980) and the thinning function was de®ned as a conditional probability of mortality in which both subject tree characteristics and subject tree neighborhood variables serve as predictors. Tree mortality data for most plots cover the 9-year period after silvicultural invervention; however, for clear cut plots mortality covers the 3-year period between 6 and 9 years after silvicultural intervention. The variables used as predictors were: DBH
ABA
RBA PBAF1 PBAF2 PBAF3 PBASF1 PBALF1
APBAF1

Diameter at breast height. Individual tree basal area growth in terms of annual absolute growth rate. Individual tree basal area growth in terms of relative growth rate (
ABA/BA). Basal area point estimate, basal area factor (F) 1, taking each tree as sampling point. Basal area point estimate, Fˆ2. Basal area point estimate, Fˆ3. Basal area point estimate (Fˆ1) of trees with DBH smaller than the subject trees. Basal area point estimate (Fˆ1) of trees with DBH larger than the subject trees. Basal area growth (annual absolute rate) of the trees included in the basal area point estimate (Fˆ1) around the subject.

The ®rst three variables are individual tree attributes, while the others characterize the tree neighborhood based on: (1) the basal area point estimates with different basal area factors (PBAF1 , PBAF2 and PBAF3 ; (2) the basal area of trees larger and smaller than the subject tree PBASF1 and PBALF1 and (3) the basal ara growth of trees in the neighborhood of the subject tree …
A PBAF1 †. These spatial variables imply that the identi®cation of a given tree as a neighbor of the subject tree depends on its size (DBH) and distance from the subject tree. In other words, the neighbor identi®cation is controlled by the basal area factor. A given tree with a DBH, Di , and at a distance Ri from the subject tree is considered a neighbor if the quantity Fi ˆ

D2i 4R2i

is greater or equal to the speci®ed basal area factor (Fˆ1,2 or 3). To ®t the thinning function as a mortality model, an approach similar to that suggested by Cressie (1991) was used, but the model was ®tted using generalized linear models with binomial distribution error and logit link function: ˆ

exp‰ 0 xŠ 1 ‡ exp‰ 0 xŠ

where is the parameter vector of the model, x is the vector of explanatory variables, as m is the expected mortality probability (McCullagh and Nelder, 1989). Starting with the full model (all predictors included), predictors with little contribution to the goodness-of®t are removed by backward elimination. Because some variables were strongly correlated, the model produced by the step-wise procedure was further reduced to retain only variables whose individual contribution to the model was signi®cant (p<0.05) and to avoid any adverse impact of collinearity. This was possible in all but two plots (lower and upper cut, plots 4 and 5). The goodness-of-®t of the models was analyzed using the c2 test for the null hypothesis that the parameter estimates do not differ from zero. This null hypothesis corresponds to the random mortality hypothesis because it implies that neither the tree attributes nor tree neighborhood are important in de®ning the probalility of mortality for individual trees. Approximate con®dence envelopes (5%) for the empirical K-function were constructed by applying the mortality model as a thinning function over the initial pattern of live trees. The spatial pattern analyzed was not that of surviving trees (outcome of the process), but rather that of mortality events (the `thinned trees'). 5. Results 5.1. Spatial pattern after clear cut 5.1.1. Tuning constants The results obtained for the test with tuning constants were somewhat different from that suggested by Diggle (1983). One possible reason is that the plot used to test the constants have a complex aggregated

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

301

Fig. 3. Empirical (thick line) and expected K-function for Neyman±Scott process with parameter estimates obtained with different values of the tuning constants c and r0. Data from clear cut plot 1.

pattern with clustering at two scales: < 3 m and around 25±30 (Fig. 3). For the limit of integration r0 ˆ 10 m, the clustering at small scales dominated the ®tting process and a large variation of parameter estimates was obtained for the different values of the power transformation constant (Table 1). Increasing the limit of integration to r0 ˆ 25 m resulted in more similar parameter estimates for different values of c. Although the ®tted models were closer to the overall pattern of the empirical K-function, different values of c still resulted in expected K-functions that were very dif-

ferent at distances larger than 20 m. The ®tted models were very similar for r0 ˆ 45 m; i.e., regardless of the c value in the power transformation, little difference was observed for the estimates of the intensity of the parent process (r) and the cluster size (s). Also, faster and better convergence was obtained for values of c equal to 0.15 and 0.25. Larger limits of integration resulted in models that ®tted the larger scale patterns better and the smaller scale patterns more poorly. These results suggest that the choice of limit of integration and power transfor-

302

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

Table 1 Values of parameter estimates (^ , ^), number of iterations (n), different values of power transformation (c) and the limit of integration (r0) for the fitting of the Neyman±Scott model. Data are from clear cut plot number 1 (6 years after treatment implementation) and starting values were sˆ13.25 and rˆ0.0014 Limit of integration r0 (m)

10 25 45

Power transformation (c) ^ ^ n ^ ^ n ^ ^ n

0.10

0.15

0.25

0.50

0.1315 0.1832 49 9.7869 0.0022 18 11.6497 0.0017 13

0.1405 0.1733 47 12.0913 1.0016 17 12.1871 0.0016 12

3.2991 0.0152 37 16.2857 0.0010 18 12.5444 0.0016 13

8.1679 0.0039 40 13.2700 0.0014 2F 12.0793 0.0016 10

F: False convergence.

mation depends on the scale which the model should represent, and that in complex patterns with clustering at different scale it might not be possible to obtain a good ®t for all scales. In the following analyses, the values of r0 ˆ 45 and c ˆ 0.15 were used in the minimization of the modi®ed least squares criterion. 5.1.2. Parameter estimates The initial parameter estimates obtained by the nonlinear regression (s(i) and r(i)) were very close to the ®nal values obtained by the minimization of the least squares criterion (^  and ^ Table 2). The non-linear regression method seems to provide good initial approximations: however, since goodness-of-®t measures are dif®cult to obtain in stochastic point processes, the improvement obtained by the modi®ed least squares was not assessed. In plots 4 and 5, the non-linear regression failed to converge, so arbitrary initial values had to be used in the minimization procedure. The cause of lack of convergence in plot 4 is probably its random pattern for all scales (except distances <2 m), so the K-function model is not appropriate for the empirical values (Fig. 4). But in plot 5, no apparent reason seems to justify the lack of convergence. The modi®ed least squares criterion converged in all plots, even in plot 4. This seems to indicate that this method is more robust than simple application of non-linear regression. The parameter estimates for plot 4 indicates that Neyman±Scott process is not an appropriate model for this plot. The estimates of parent event intensity (r)

suggest that the number of clusters is equal to the number of trees in the plot, and the average size for the cluster (r) is less than 0.05 m, so the model seems to make each individual tree a `one-tree-cluster' (Table 2). Since the parent event is an HPP, this indicates that the tree patterns is completely spatially random. The con®dence envelope for the K-function (Fig. 4), however, shows that the model overemphasized the aggregation at small scales (<3 m) and produced a pattern at large scales that is more random than the empirical K-function. Similar trends were observed in the other plots. In general, the empirical K-function is outside the 95% con®dence envelope at distances less than 0.5 m, but the strong aggregation at larger scales is well represented by the model. 5.1.3. Tree spatial pattern Despite some variation in number and cluster size among the plots 1±3, the Neyman±Scott model indicates that they have similar aggregated patterns because this variation is relatively small. Over time, these plots show decreasing aggregation (Fig. 4), and the model re¯ects these changes by the increase in number of cluster …^ † and lesser and inconsistent trends in the cluster size …^ † (Table 2). Stand density (trees/ha) and basal area also increased which indicates that the process of re-colonization after clear cutting is still an ongoing process and that the ingrowth trees tend to occupy the more open areas. The random pattern and the decrease in the number of trees/ha indicate that plot 4 is under a different

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

303

Table 2 Parameter estimates (^ , ^) of te Neyman±Scott process with Gaussian distribution for different clear cut plots (treatment c) after 6 and 9 years of treatment implementation Years after clear cut

Plots 1

2

3

4

5

6

…i† …i† ^ ^ ^ np Trees/ha BA (m2/ha)

10.1014 0.0019 12.1817 0.0016 8.0 454 4.10

11.1598 0.0011 15.6394 0.0008 4.0 252 1.66

10.9583 0.0016 8.8488 1.0019 9.5 166 1.38

10a 0.001a 0.0434 0.0434 217.0 434 3.59

10.1a 0.0019a 7.1443 0.0082 40.9 334 2.09

9

…i† …i† ^ ^ np ^ Tree/ha BA (m2/ha)

10.2995 0.0025 10.3551 0.0025 12.3 470 4.77

11.2503 0.0012 16.5355 0.0008 4.2 364 2.72

6.6167 0.0043 8.8127 0.0031 15.6 282 2.29

10.1a 0.0019a 0.0420 0.0420 210.0 420 3.81

10a 0.001a 9.5052 0.0038 18.8 486 3.37

a Arbitrary values: non-linear regression failed to converge. ^…i† , ^…i† : Starting values for the minimization process. ^np : Estimated number of cluster. BA: Plot basal area.

process. The `stand initiation' phase (Oliver and Larson, 1990) is probably ®nished and the inter-tree competition is generating density-dependent mortality and a decrease is stand density. Over time, plot 5 is the only plot that shows an increase in aggregation in the K-function pattern (Fig. 4). The Neyman±Scott parameter estimates con®rm this trend with a noticeable decrease in the number of clusters and more subtle increase in cluster size (^  and ^ Table 2). Nine years after implementing the clear cut, the plots are in different stage of the `stand initiation' phase. This is indicated by the variation in the tree pattern among plots and the different directions of the temporal change in this pattern. In plot 5, the recolonization is still increasing in speed and thereby increasing both stand density and tree aggregation, while in plot 4, this phase ended and a phase of more intense plant competition has probably started. Plots 1±3 seem to be in an intermediate stage when stand density and basal area still increase but tree clustering declines from occupation of more open areas by the newly established trees. The Neyman±Scott model is able to capture these trends by differences in the parameter estimates. The large variation in number of trees/ha and basal area

among plots with similar spatial pattern (plots 1, 2 and 3) indicates that plot level information that is not spatially explicit might not able to characterize this process. 5.2. Spatial pattern of ingrowth trees The IPP model suggested more random ingrowth pattern for control and BA reduction plots as a result of their more homogeneous canopy structures. In contrast, the irregular canopies of the clear cutting and lower±upper cut resulted in a more aggregated pattern in these silvicultural interventions (Figs. 5±8). In the control plots, the predicted pattern was more clustered than the observed pattern for plots 2 and 4, but was more random for plots 1 and 5 (Fig. 5). Similar variations were observed in the other treatments. The predicted patterns tended to be closer to random in BA reduction than in control, which suggests higher canopy uniformity with light silvicultural intervention than in unmanaged plots (Fig. 6). Con®dence envelopes generally predicted aggregated pattern for lower±upper cut and clear cut (Figs. 7 and 8, respectively). Although this predicted pattern seems close to the observed pattern under

304

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

Fig. 4. Established trees: confidence envelopes (95%) for the linear transformation of the K-function (L(r)- r) for the clear cut plots assuming a Neyman±Scott process (dotted lines). Solid thick line is the empirical K-function for 6 years after treatment implementation and dashed line is for 9 years. The confidence envelope is based on 19 simulations of the model using parameter estimates obtained for the pattern 6 years after treatment.

lower±upper cut plots, ingrowth in the clear plots showed an aggregation pattern that departed from the prediction. In this case, the IPP model predicted ingrowth in relatively small, randomly distributed in the plots but the observed pattern was either random

(plots 1 and 4) or had aggregation at larger scales (plots 2, 3 and 5). This result suggests that canopy openness is not the only factor controlling ingrowth distribution; other factors related to seed dispersion and environmental heterogeneity both above- and

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

305

Fig. 5. Ingrowth: confidence envelopes (95%) for the linear transformation of the K-function in control plots. Confidence envelopes are based on 19 simulations of the inhomogeneous Poisson process, with intensity function inversely proportional to the basal area point estimate (PBA) with basal area factor (F) 3.

306

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

Fig. 6. Ingrowth: confidence envelopes (95%) for the linear transformation of the K-function in BA reduction plots. Confidence envelopes are based on 19 simulations of the inhomogeneous Poisson process, with intensity function inversely proportional to the basal area point estimate (PBA) with basal area factor (F) 3.

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

307

Fig. 7. Ingrowth: confidence envelopes (95%) for the linear transformation of the K-function in lower-upper cut plots. Confidence envelopes are based on 19 simulations of the inhomogeneous Poisson process, with intensity function inversely proportional to the basal area point estimate (PBA) with basal area factor (F) 3.

308

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

Fig. 8. Confidence envelopes (95%) for the linear transformation of the K-function (L(r) - r) in clear cut plots. Confidence envelopes are based on 19 simulations of the inhomogeneous Poisson process, with intensity function inversely proportional to the basal area point estimate (PBA) of basal are factor (F) 3.

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

below-ground are also probably important. Among the plots that received the lower±upper cut treatment, the only plot that the IPP model fails to predict well is plot 1. Despite imperfect performance in the other plots, the model shows that the canopy structure does have a strong in¯uence on ingrowth pattern. This might result from the elimination of advanced regeneration in this treatment (trees with DBH<10 cm), so that the pattern of ingrowth trees was dictated by the pattern of younger regeneration, which would be more responsive to current canopy structure.

309

5.3. Tree mortality The c2 test of the model against the random mortality model shows that mortality was positively related to tree size (DBH) or absolute basal area growth (
A BA), but the strong correlation between these two variables suggests that
A BA also represents the tree size (Table 3). This relation is consistent across different types of silvicultural interventions and indicates that trees close to or in the main canopy (larger trees) are probably under high levels of com-

Table 3 Variable included in the models for prediction of individual tree mortality for the different treatment and plot (G indicates a general model for all plots within treatment) Treat./plot

A

B

C

D

E

a

1 2 3 4 5 G 1 2 3 4 5b Gb 1 2 3 4a 5a G 1 2 3 4 5 G 1 2 3 4 5 G

n

251 303 300 194 251 1299 246 253 260 325 282 1366 116 204 214 200 214 948 227 126 83 217 167 820 202 202 265 188 202 1059

m

25 60 41 34 34 194 18 26 23 28 23 118 28 16 15 24 25 108 24 9 10 25 11 79 22 22 30 20 34 128

Fitted model Tree attributes

Tree neighborhood

‡
ABA ‡DBH ÿ
RBA ‡
ABAÿ
RBA Ð ‡DBH Ð ÿ
RBA ÿ
RBA ÿ
RBA ‡DBHÿ
RBA ÿ
RBA ÿ
ABA ÿ
RBA ÿ
RBA ‡DBH‡
ABA ‡DBH ‡DBH‡
ABAÿ
RBA Ð Ð Ð ‡DBH Ð ÿDBH Ð ‡
ABA ÿ
RBA ‡DBH ‡DBHÿ
RBA ‡DBHÿ
RBA

Ð ‡
APBAF1 Ð ÿPBASF1 ‡
APBAF1 ÿPBASF1 Ð Ð Ð ÿ[
APBAF1/PBAF2] ‡[ÿPBALF1 /PBAF1] Ð ‡PBAF3 Ð ‡PBALF1 ‡PBAF1 ÿ PBALF1 ÿ PBASF1 ‡PBAF2ÿPBAF3 ÿPBASF1 ÿPBASF1 ÿPBASF1 Ð Ð ‡PBAF1 Ð Ð Ð ‡
APBAF1 Ð Ð ‡PBAF2 ÿPBALF1 ÿ PBASF1 ‡
A PBAF1

Model with multicolinearity problems. Unstable model: fitted probabilities are close to 0 or 1. n: Number of observations used in fitting the model. m: Number of dead trees observed. The last column is the p-value of the c2 test of the `thinning process' model against the random mortality model. b

c2 Pr(H0) 0.0279 0.0228 0.0056 0.0007 0.0330 0.0000 Ð 0.0225 0.0568 0.0212 0.0000 0.0000 0.0003 0.0198 0.0000 0.0009 0.0024 0.0000 0.0025 Ð Ð 0.0573 Ð 0.0300 Ð 0.0006 0.0032 0.0541 0.0026 0.0000

310

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

Fig. 9. Confidence envelopes (95%) for the linear transformation of the K-function describing the spatial pattern of mortality trees. Confidence envelopes are based on the thinning of the initial spatial pattern of the plot using the respective mortality model.

petitive/environmental stress than the smaller tree in the sub-canopy. The differences in species composition between the understory and overstory might also in¯uence this relationship; while understory trees are mainly shade tolerant, shade intolerant species can reach the forest canopy by regeneration in large gaps. Shade intolerant species tend to have a short life-span, increasing the mortality probability of larger trees.

The negative relation between mortality and relative growth rate …
R BA† probably re¯ects differences in vigor and in the stress level that individual trees are subjected to. High levels of stress are usually translated into smaller relative growth rates and hence higher probability of mortality. Variables of neighborhood description also showed signi®cant relationships with mortality in several

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

311

Fig. 9. (continued)

plots. High levels of density in the subject tree neighborhood (PBAF1 , PBAF2 PBAF3 ), the density of trees larger than the subject trees (PBALF1 and PBALF1 =PBAF2 ) and the total basal area growth in the neighborhood(
A PBAF1 ) were generally associated with higher probabilities of mortality. This result indicates that light limitation drives competition-related stress, and that inter-plant competitions is one-sided. Models ®tted at treatment level (all plots pooled within the treatment) were signi®cant when compared against the random mortality models

(Table 3) and the most consistent predictor variable of tree mortality was subject tree DBH. However, tree neighborhood variables were not always present in these models and their prediction of the spatial pattern in mortality was similar to the random mortality model, resulting in a performance worse than the plot level models. The random mortality hypothesis seems appropriate for clear cut plots, and suggests that these plots are still in a re-colonization phase with little densitydependent mortality; however, plot 4 is probably an

312

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

exception. The importance of spatial variables in predicting tree mortality is somewhat greater in control plots than in plots with partial cuts, but high variability in the models among plots within the treatments precludes generalization about the effect of reduction in stand density and tree mortality pattern. The importance of neighborhood description in the prediction of individual tree mortality was strongly dependent on the speci®c stand structure of individual plots. Fig. 9 shows 95% approximate con®dence envelopes for the K-function generated by thinned processes ®tted to mortality pattern in two plots of each treatment. In the ®rst of each pair, the random mortality hypothesis was not rejected (spatial or neighborhood variables are not present in the mortality function), and in the other plot the mortality function included spatial variables. The signi®cance of the c2 test of goodness-®t does not translate directly in maintenance of the K-function within the ®tted models' con®dence envelopes, but the model tends to produce nonrandom patterns for the con®dence envelopes, even if the pattern of the K-function is not exactly represented. The random mortality hypothesis or random thinning is applied on the original pattern of the live trees, so that it does not necessarily imply that the pattern of mortality trees will or will not follow CSR. It seems that tree mortality can be reasonably modeled as a thinning process over the pattern of live trees where the thinning function is dependent on subject tree attributes and subject tree neighborhood. Although the models failed to represent the pattern of mortality trees in some plots (for exemple, BA reduction plot 4, clear cut plot 1, and selective cut plot 3), the models produced surprisingly accurate results considering that the only information used was the tree map and DBH. Information on tree height, tree crown and forest canopy could certainly improve the performance of this model. 6. Conclusions The Neyman±Scott process is an appropriate model for simulating the pattern of live trees following clear cutting of tropical forests. These stands are still undergoing the process of re-colonization and, hence, still

show different degrees of clustering. This process represents a re®nement of the analysis of the spatial pattern because it allows the quantitative description of aggregated patterns in terms of number and size of clusters. This information allows the distinction of plots that have similar K-function graphics but that can be undergoing different ecological processes. It is also ¯exible enough to ®t the CSR patterns, and the parameter estimates clearly signal the prevalence of this spatial pattern. Other biological justi®cation for the Neyman±Scott process in modeling the spatial pattern of re-colonizing trees include  the `seed rain' tends to be clustered in these areas because many early secondary successional species and shade intolerant species are dispersed by birds that visit open areas and use established trees as perching sites;  spatial variation in density and seed viability in the soil seed bank would result in a clustered pattern of regeneration;  environmental heterogeneity that occurs naturally or that is produced by wood debris from logging operations can create microsite variation relative to seed germination and seedling establishment; These effects can result in a delay of tree establishment in certain areas of the plot, but over time, complete occupation tends to reduce or eliminate initially aggregated patterns. The inhomogeneous Poisson process seems to be a useful model to study the in¯uence of canopy structure on understory plants, and better intensity functions can be built if direct information on canopy structure is available. Geostatistics have been used to modeling canopy±understory relationships by constructing a correlation model between canopy structure and understory plant density (Moeur, 1993), or between canopy and soil variables and regeneration (Kuuluvainen et al., 1993). If canopy openness and soil properties can be considered continuous spatial variables, the representation of plant spatial pattern by spatially continuous density functions or discrete density variables might be an over-simpli®cation. IPP provides a good framework for modeling the spatial pattern of discrete events (plants) as a function of a spatially continuous variation in the environment (canopy openness, soil variables).

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

It seems that tree mortality can be reasonably modeled as a thinning process conditional on the pattern of live trees. The thinning function under such scenarios is dependent on subject tree attributes and subject tree neighborhood. Although the models failed to represent the pattern of mortality trees in some plots, the models produced biologically reasonable results that could be improved upon if information on tree height, tree crown and forest canopy were available for describing stand structure more accurately. Acknowledgements The authors thank Mr. Renato M. Jesus from the Cia. Vale do Rio Doce (CVRD) and the Linhares Reserve for making his data available and the Brazilian National Council for Scienti®c Development (CNPq) for a scholarship for the ®rst author during the development of this work. References Augspurger, C., 1983. Offspring recruitment around tropical trees: changes in cohort distance with time. Oikos 40(2), 189±260. Awasthi, A., 1990. Spatial distribution pattern of Diospyros melanoxylon. J. Trop. Ecol. 31(1), 69±74. Bariteau, M., 1992. ReÂgeÂneÂration naturelle de la foreÃt tropicale humide de Guyane: eÂtude della reÂpartition spatial de Qualea rosea Aublet, Eperua falcata Aublet et Symphonia globulifera Linnaeus f. Ann. Sci. ForestieÁre 49, 359±382. Batista, J., 1984. Spatial dynamics of trees in a Brazilian Atlantic tropical forest under natural and managed conditions. Ph.D. Thesis, University of Washington, Seattle. Busing, R., 1991. A spatial model of forest dynamics. Vegetatio 92, 167±179. Clark, D.A., Clark, D.B., 1994. Spacing dynamics of a tropical rain forest trees: evaluation of the Janzen±Connell model. Am. Naturalist 124(6), 769±788. Coley, P.D., 1983. Herbivory and defensive characteristics of tree species in a lowland tropical forest. Ecol. Monogr. 53(2), 209± 233. Collins, S.L., Klahr, S.C., 1991. Tree dispersion in oak-dominated forests along an environmental gradient. Oecologia 86(4), 471± 477. Condit, R., Hubbell, S.P., Foster, R.B., 1992. Recruitment near conspecific adults and the maintenance of tree and shrub diversity in a neotropical forest. Am. Naturalist 140(2), 261± 286. Connell, J.H., 1970. On the role of natural enemies in preventing competitive exclusion in some marine animals and in rain forest

313

trees. In: den Boer, P.J., Gradwell, G.R. (Eds.). Dynamics of Populations: Proc. Advanced Study Institute on Dynamics of Numbers in Populations. Center for Agriculture Publishing and Documentation, Wageningen, Netherlands, pp. 298±310. Cressie, N., 1991. Statistics for Spatial Data. Wiley, New York. Diggle, P., 1983. Statistical Analysis of Spatial Point Patterns. Academic Press, London. Franklin, J., Michaelsen, J., Strahler, A.H., 1985. Spatial analysis of density dependent pattern in coniferous forest stands. Vegetation 64, 29±36. Gibson, N., Brown, M.J., 1991. The ecology of Lagarostrobos franklinii (Hook. f.) Quinn (Podocarpacease) in Tasmania. 2. Population structure and spatial pattern. Austr. J. Ecol. 16(2), 223±229. Guttorp, P., 1995. Stochastic Modelling of Scientific Data. Chapman & Hall, London. Haase, P., 1995. Spatial pattern analysis in ecology based on Ripley's K-function: Introduction and methods of edge correction. J. Veg. Sci. 6, 575±582. Hubbell, S.P., 1979. Tree dispersion, abundance, and diversity in a tropical dry forest. Science 203(4387), 1299±1309. Hubbell, S.P., 1980. Seed predation and the coexistence of tree species in tropical forests. Oikos 35(2), 214±229. Husch, B., Miller, C.I., Beers, T.W., 1982. Forest Mensuration, 3rd ed. Wiley, New York. Janzen, D.H., 1970. Herbivores and the number of tree species in tropical forests. Am. Naturalist 104(940), 501±528. Jesus, R.M., 1987. Mata AtlaÃntica de Linhares: apectos florestais. In: Anais do SeminaÂrio: Desenvolvimento EconoÃmico e  reas do TroÂpico U  mido Brasileiro ± Impacto Ambinetal em A A ExcperieÃncia da CVRD. SEMA/IWRB/Cia. Vale do Rio Doce, Rio De Janeiro, pp. 35±53. Jesus, R.M., Souza, A.L., Garcia, A., 1992. Produc,aÄo SustentaÂvel de Floresta AtlaÃntica. Sociedade de Investigac,oÄes CientõÂficas, Universidade de Vic,osam, Vic,osa. Johnson, L., Riess, R., 1982. Numerical Analysis. AddisionWesley, New York. Kenkel, N.C., 1988. Pattern of self-thinning in jack pine: testing the random mortality hypothesis. Ecology 69(4), 1017±1024. Kenel, N.C., Hoskins, J.A., Hoskins, W.D., 1989. Local competition in naturally established jack pine stand. Can. J. Bot. 67, 2630±2635. Kuuluvainen, T., Hokkanen, T.J., JaÈvinen, E., Pukkala, T., 1993. Factors related to seedling growth in a boreal scots pine stand: a spatial analysis of vegetation-soil system. Canad. J. For. Res. 23, 2101±2109. Maguire, D., Batista, J.L.F., McKenzie, D., 1993. Horizontal structure of uneven-aged mixed species forests modeled as an Inhomogeneous Poisson process. In: Reynolds, K. (Ed.). Proc. IUFRO S4.11 Conf. on Stochastic Spatial Models in Forestry. University of Greenwich, Thessaloniki, Greece, pp. 163±170. McCullagh, P., Nelder, J.A., 1989. Generalized Linear Models. Chapman & Hall, Cambridge. Moeur, M., 1993. Characterizing spatial patterns of trees using stem-mapped data. For. Sci. 39(4), 756±775. Oliver, C.D., Larson, B.C., 1990. Forest stand dynamics. McGrawHill, New York.

314

J.L.F. Batista, D.A. Maguire / Forest Ecology and Management 110 (1998) 293±314

Osher, J., 1983. On estimators for the reduced second moment measure of pointo process. Statistics 14(1), 63±71. Osher, J., Stoyan, D., 1981. On the second-order and orientation analysis of planar stationary point process. Biometr. J. 23(6), 253±333. Penttinen, A., Stoyan, D., Henttonen, H.M., 1992. Marked point process in forest statistics. For. Sci. 38(4), 806±824. Ripley, B., 1977. Modelling spatial pattern. J. Roy. Statist. Soc. 39, 172±242. Rizzini, C.T., 1979. Tratado de Fitogeografia do Brasil. Higitec/ Universidade de SaÄo Paulo, SaÄo Paulo. Skarpe, C., 1991. Spatial patterns and dynamics of woody vegetation in an arid savanna. J. Veg. Sci. 2, 565±572. Sterner, R.W., Ribic, C.A., Schatz, G.E., 1986. Testing for life historical changes in spatial patterns of four tropical tree

species. J. Ecol. 74, 621±633. Szwagrzyk, J., Czerwczak, M., 1993. Spatial pattern of trees in natural forests of East-Central Europe. J. Veg. Sci. 4, 469±476. Veloso, H., GoÂes-Filho, L., 1986. Fitogeografia Brasileira: Classificac,aÄo FisionoÃmica-EcoloÂgica de Vegetac,aÄo Neotropical. RADAM Brasil, Salvador. Welden, C., Slauson, W., Ward, R., 1990. Spatial pattern and interference in PinÄon-Juniper woodlands of northwest Colorado. Great Basin Naturalist 50(4), 313±319. West, P.W., 1984. Inter-tree Competition and small-scale pattern in monoculture of Eucalyptus obilqua L'Herit. Austr. J. Ecol. 9, 405±411. Wong, M.S.J., Wright, S.P.H., Foster, R.B., 1990. The spatial pattern and reproductive consequence of outbreak defoliation in Quararibea asterolepis, a tropical tree. J. Ecol. 78, 579±588.