Fitting Ecological Process Models to Spatial Patterns Using Scalewise Variances and Moment Equations
Abstract
Ecological spatial patterns are structured by a multiplicity of processes acting over a wide range of scales. We propose a new method, based on the scalewise variance—that is, the variance as a function of spatial scale, calculated here with wavelet kernel functions—to disentangle the signature of processes that act at different and similar scales on observed spatial patterns. We derive exact and approximate analytical solutions for the expected scalewise variance under different individual-based, spatially explicit models for sessile organisms (e.g., plants), using moment equations. We further determine the probability distribution of independently observed scalewise variances for a given expectation, including complete spatial randomness. Thus, we provide a new analytical test of the null model of spatial randomness to understand at which scales, if any, the variance departs significantly from randomness. We also derive the likelihood function that is needed to estimate parameters of spatial models and their uncertainties from observed patterns. The methods are demonstrated through numerical examples and case studies of four tropical tree species on Barro Colorado Island, Panama. The methods developed here constitute powerful new tools for investigating effects of ecological processes on spatial point patterns and for statistical inference of process models from spatial patterns.
Introduction
Nonrandom spatial patterns are ubiquitous in ecology and provide important information on processes structuring communities, including dispersal, competition, density dependence, and habitat associations (Law et al. 2009; McIntire and Fajardo 2009). However, similar spatial distributions can emerge from different mechanisms or different combinations of mechanisms. A first step in disentangling different processes of pattern formation is to define the scales at which the pattern occurs and thereby link the pattern to some mechanism operating at these scales (Kershaw 1963; Levin 2000). Large-scale variability (∼0.1–10 km), for example, may be driven by topographic features or soil heterogeneity operating at the same scales. Small-scale structure (<∼100 m) may be related to local dispersal and interindividual interaction.
Most of the widely used tools for characterizing spatial patterns have limited ability to separate scales. For example, Ripley’s K, pair correlation densities, and functions related to second-order product densities more generally are designed to estimate average local densities around the focal individual as a function of distance. By their cumulative nature, values at any distance reflect the combined effects of mechanisms acting at different scales (Loosmore and Ford 2006). For example, aggregation at short distances is found in a habitat specialist plant with long dispersal distance (fig. 1A) as well as a generalist with short dispersal (fig. 1B) and, of course, a specialist with short dispersal (fig. 1B). This limits our ability to investigate basic ecological mechanisms contributing to pattern formation. Contrasting spatial patterns, together with their estimated pair correlation densities and wavelet variances. A, Poisson process, a model of global dispersal in a homogenous (left) and inhomogeneous (right) habitat. B, Simulated individual-based model with short dispersal in a homogenous and inhomogeneous habitat. Note that the wavelet variance, unlike the pair correlation density, separates processes by scale. For example, habitat association at larger scales will elevate pair correlation statistics at all smaller distances, while the wavelet variance is substantially unaltered for scales at which habitat association is not operating.
To overcome some of these limitations, we need to move from a framework where variability is analyzed as a function of space to a new approach where the variance is decomposed scale by scale. Spectral analysis is such an approach, and it allows for a clear view of how different frequencies contribute to the pattern of interest (Bartlett 1964; Platt and Denman 1975). Here, we employ spectral analysis and focus in particular on the scalewise variance using wavelets, a second-order, consistent estimator of the spectral density (Percival 1995). Wavelet transforms have been employed in several studies of ecological time series (Cazelles et al. 2008; Detto et al. 2012) but less in spatial ecology (e.g., Dale and Mah 1998; Keitt and Urban 2005; Mi et al. 2005), and we are not aware of any applications in point process analysis. The mathematical and statistical properties of wavelet transforms are well understood from applications in other fields (e.g., Kumar and Foufoula-Georgiou 1997). The advantage of this technique is that it makes it possible to isolate the scales of the process of interest and analyze them separately, even in the presence of nonstationarity (i.e., habitat heterogeneity, a case that is not addressed in this study).
Individual-based, spatially explicit models are another key tool in exploring the determination of spatial patterns (Pacala 1986; Durrett and Levin 1994). Such models capture spatiotemporal dynamics with a plant’s-eye view or, more generally, an individual-based view. Dispersal and spatial interactions among individuals located a specified distance apart are represented by kernel functions describing the mechanisms that govern the interaction or the dispersal/movement event. Moment closure can be applied to describe the dynamics of mean densities and pair correlation densities and thereby obtain analytical results (Bolker and Pacala 1997; Law and Dieckmann 2000). These methods have improved our understanding of the structure of plant communities, showing how endogenous aggregation patterns affect species interactions and coexistence (Murrell and Law 2002; Bolker et al. 2003). Moment equations can also be developed—and, indeed, are more easily handled—in the Fourier transformed space (Bolker et al. 2000; Ovaskainen and Cornell 2006b). However, little has been done to apply moment methods to real data in a statistical inference framework. Even in cases when the mechanisms of pattern formation are well understood, the quantification of their effects remains a formidable challenge (McIntire and Fajardo 2009; Thompson and Katul 2011).
One of the fundamental problems in pattern analysis concerns understanding the limits to the information that can be extracted from a single random sample, including how many processes can be unambiguously discerned. For example, the same static patterns of the distributions of tropical tree species in the large forest plot (50 ha) of Barro Colorado Island have alternatively been used to investigate habitat associations (Plotkin et al. 2000; Harms et al. 2001), effective dispersal distance (Anand and Langille 2010), and negative density dependence (Bagchi et al. 2011). This raises the question of how all these processes can be analyzed simultaneously without interference of one process with another.
Here, we show how scalewise variances and, specifically, wavelet variances can be derived analytically from spatially explicit, individual-based models. We derive new analytical results on wavelet variances expected under models with local dispersal and neighborhood competition, using moment methods in homogeneous environments. We further derive the asymptotic statistical properties of these wavelet variances and show how these can be used to estimate process parameters from observed spatial patterns. A key result is the derivation of an approximate likelihood function that is consistent with the process underlying the model and makes it possible to estimate the most likely parameters and, importantly, their uncertainty. We demonstrate these methods with numerical examples and a case study of four tropical tree species. The methods presented here constitute a powerful new tool for ecologists to investigate spatial processes and patterns.
Theory
The Scalewise Variance: A Tool for Scale-By-Scale Decomposition of Spatial Patterns
Consider a spatial point process that generates a random pattern of points that are the locations of individuals in a two-dimensional space of finite area A. A realization of this process can be represented as the sum of Dirac δ functions centered at the n point locations 
and 0 everywhere else: 
.
The first spatial moment, the mean density, is simply
Application of a Fourier transform, 
, reduces the pair correlation density to a simple product minus a constant (the Fourier transform of a δ function is unity):
is the angular frequency, 
is the Fourier transform of p(x) (i.e., 
, with j the imaginary unit), and 
is its complex conjugate. The first term on the right-hand side of equation (3) is the two-dimensional spectral density function of p(x), or simply the power spectrum, 
(Bartlett 1964). Knowledge of the function 
makes it possible to assess the scalewise variance Vψ, that is, the expected variance in population density sampled with a given kernel ψ(x) (Daley and Vere-Jones 2003, p. 304):
, where 
, so equation (4) can be rewritten as 
. The functions ψ are essentially smoothing filters whose properties depend on the type of kernel and associated scaling parameter, s. They enable study of features of the process whose detail matches their scale parameters, that is, broad features for large s and fine features for small s. For this purpose, kernels should be functions with compact support in both the frequency and the spatial domain, essentially compromising between frequency and spatial resolution in order not to lose localization in physical space (in accordance with the Heisenberg uncertainty principle, which in this context states that it is impossible to achieve spectral and spatial resolution simultaneously). The scale parameter s is associated with a physical scale to allow ready interpretation; for example, for isotropic problems, it is convenient to express s as a function of the one-dimensional Fourier period 
, which represents the distance among peaks of an oscillation.Wavelets constitute a suitable family of kernel functions satisfying the above requirements (Kumar and Foufoula-Georgiou 1997). Among many choices of wavelets, the Morlet (Daubechies 1992) is one of the most commonly used wavelets and has proved effective in spatial ecological applications (Mi et al. 2005). It can assume a simplified isotropic form, and its Fourier transform is
, the scalewise variance Vψ(s) is basically identical to its Fourier counterpart when the spectrum is relatively featureless (e.g., it obeys a power law over a certain interval of frequencies).For this wavelet, the ratio between the equivalent Fourier period (λ), and the wavelet scales (s), Ff exists and can be derived analytically (Meyers et al. 1993); we find that it is
For more background on wavelets and wavelet variances, see Percival (1995), Kumar and Foufoula-Georgiou (1997), and Torrence and Compo (1998).
Statistical Properties of the Scalewise Variance Estimator
For complete spatial random (CSR) distributed populations, the normalized scalewise variance is equal to 1 at all scales, just as for a Poisson process, the variance is equal to the mean density. Eventually, at small scales, the normalized variance of any observed pattern tends to 1; as the probability to encounter two individuals in the same area becomes negligible, the variance approaches the mean, and locally the process approaches a Poisson process (Daley and Vere-Jones 2003).
The null hypothesis of complete spatial randomness can be tested using the properties of the Poisson distribution. For a CSR process, the number of individuals contained in a box of area 
is a Poisson variable with mean equal to 
. Given k independent realizations X1, X2, …, Xk, drawn from a Poisson distribution with mean 
, the sum of k normalized square deviations is approximately equal in distribution to a χ2 with 
degrees of freedom:
.Our work shows that the same approach can also be used to estimate the confidence intervals of the wavelet variance for a CSR process. In this case, the determination of the degrees of freedom is not so straightforward because the wavelet is not perfectly localized in physical space (in contrast with indicator functions). The degrees of freedom depend on the number of independent points averaged by the wavelet filter, which varies with the analyzed scale. We would expect that 
in the same manner as the box counting method. Large sample theory ensures that, under fairly general conditions (Serroukh et al. 2000), the estimator 
is asymptotically normally distributed with mean 
and finite large sample variance. For a Gaussian process, the large sample variance is equal to 
(Percival 1995), where
is distributed as a χ2 with η degrees of freedom given by
, is asymptotically (complex) Gaussian because it is a sum of independent, identically distributed random variables (for 
). Solving the integrals (4) and (8) for the Morlet wavelet with the conditions 
and 
, we find
, as expected from the box counting analogy. Monte Carlo simulations agree very well with our resulting predicted confidence intervals (fig. 2). According to equation (10), the large sample properties of the normalized wavelet variance for a CSR process are independent of the mean density (provided that there are at least five individuals in the sampled area), a result confirmed by simulations. This makes the wavelet variance a pivotal discriminant statistic for testing the null hypothesis of nonrandomness for populations with different numbers of individuals.
A, Wavelet variance estimates from 1,000 simulated point patterns of 500 individuals over 1 ha displaying complete spatial randomness (gray lines), compared with the expected value (black line) and predicted 95% confidence intervals (red lines). B–E, Predicted and observed quantiles for different scales. Note that the width of the confidence interval increases with increasing scale (red lines in A), because the number of degrees of freedom declines with increasing scale. The quantiles predicted using the χ2 approximation closely match the observed.
It is important to note that the wavelet variance estimates at different scales are not completely independent, because the functions 
overlap for small spectral intervals. This can create complications when the least squares regressions or likelihood functions are computed. One approach would be to choose the interval between scales sufficiently large to minimize this problem, but this has the drawback of losing resolution and important features of the wavelet variance. Alternatively, under the assumption that for two scales λi and λj, wavelet variances computed from independent realizations are asymptotically normally distributed with mean 
and 
, the joint probability of observed wavelet variances can be approximated by a multivariate normal distribution 
(fig. S1). We thus obtain the (new) result that for a Gaussian process, the covariance matrix can be defined, similar to the large sample variance above, as 
, where
is not known a priori, an estimate 
can be used instead (Percival 1995).Model Fitting
The statistical properties of the observed scalewise variance enable parameter estimation using maximum likelihood techniques. Let 
be the observed wavelet variance. Then, given a parametric model defining the expected wavelet variance 
as a function of a parameter vector β, and assuming that the estimated wavelet variance deviations 
are jointly normal distributed with dispersion matrix Σ, we find that the likelihood function L is directly proportional to
is known analytically or by numerical integration, maximum likelihood estimation can be done using a generalized nonlinear least squares algorithm (app. B). More generally, it can be done using stochastic methods, that is, by simulating point patterns for specified parameters (for a review of such stochastic methods, see Hartig et al. 2011).The dispersion matrix Σ is known for Gaussian processes from 
and equation (11); however, it is unknown for the general case. A first estimate of the parameters can be made using the Gaussian assumption. In a second step, simulations should give the correct Σ. This process can be repeated iteratively until model parameters and their uncertainties converge. We found that one iteration was usually enough. The likelihood function (eq. [12]) is a new result that constitutes the foundation of the fitting method.
Models
A General Spatial Logistic Model for Population Dynamics and Its Expected Scalewise Variance
Consider a spatial-temporal process p(x, t) that represents the dynamics of a single population of identical sessile individuals, such as plants, on a two-dimensional spatial domain. Individuals die at rate m and reproduce at rate f, and propagules (e.g., seeds) are dispersed from the parent according to dispersal kernel D(x). The establishment probability of a propagule depends on the positions of all the conspecific individuals in its neighborhood, with the influences of neighbors described by interaction kernel K(x). At each time t, the process is regarded as a realization of a point pattern process and described by the function p(x, t). The probability of finding a new individual at location x is given by the convolution 
. The density-dependent probability of establishment of new individuals due to competition with neighbors (including parents) is expressed by the convolution 
.
Assuming isotropy of the dispersal and competition kernels and a homogeneous environment, the temporal dynamics of the first and second moments (Bolker et al. 2000; Dieckmann and Law 2000) follow
is the third moment (for details, see app. A).We obtain the exact analytical solution for the expected scalewise variances of this model under the simplified case, in which establishment is purely spatial independent (model I), and the approximate solution for the full model (model II).
Model I: Dispersal Only
If establishment is purely spatial independent (e.g., 
), then the spatial pattern is due entirely to dispersal. The equation for the first moment reduces to the classic logistic model
, is
).Equation (17) can be used to estimate the dispersal parameter from the observed scalewise variance for any given dispersal kernel. This is completely straightforward if the Fourier transform of the dispersal kernel (i.e., its characteristic function) is known analytically (Table 1) and can be done numerically in other cases. The wavelet variance provides a clear signature of the effects of different dispersal parameters (fig. 3A) and also differs among different dispersal kernels (fig. 3B).
| Name | Dispersal kernel  | Fourier transform of kernel  | Gaussian factor Gf |
|---|---|---|---|
| Gaussian |  |  | 1 |
| Exponential |  |  |  |
| Laplace |  |  |  |
| Uniform |
|  | 2 |
| Epanechnikov |
|  |  |
| Biweight |
|  |  |
Wavelet variances for model I, in which spatial patterns are driven entirely by dispersal. A, Analytical solutions for a two-dimensional Laplace dispersal kernel under different dispersal parameters. B, Analytical solutions for different dispersal kernels as a function of the normalized Fourier period, where Gf is the equivalent Gaussian factor (Table 1). C, Observed wavelet variances for 1,000 simulations on a square arena of 
ha with Gaussian dispersal kernel, with parameter values 
m and 
(gray lines), with their average (black line) and 95% confidence intervals (CIs; black dashed lines under red line), compared with the analytical model (red line). D, 95% CIs of 1,000 simulations with 
m and 
ha, relative to the theoretical 95% CIs for Gaussian processes (dashed red lines).
One important property of this dispersal-only model is that for large scales (
), the logarithmic slope of the variance is independent of the type of kernel and asymptotically approaches 2. This suggests defining an equivalent Gaussian factor for each kernel such that for large scales all the curves collapse to the solution obtained with the Gaussian kernel. The Gaussian factor for a particular kernel 
is obtained analytically from the limit for 
of the ratio
This model is simulated starting from an initial random pattern of n individuals on a torus of area A (or larger to minimize edge effects). A new individual is generated at a random distance, chosen according to the dispersal kernel distribution D(r, cD), from a randomly chosen parent. Then, a randomly chosen individual is removed from the point pattern. Simulation of births and deaths continues until the second-order spatial moment is approximately constant.
Simulations show that the analytical model correctly captures the mean spatial pattern (fig. 4C). Unlike the Poisson process, for patterns generated by this cluster process, the confidence intervals around the wavelet variance estimates—and thus around Wavelet variance for model II, in which both dispersal and density-dependent establishment affect spatial patterns. A, B, Analytical solutions under varying values of the nondimensional parameters 
—depend on the number of individuals n (fig. 3D). As 
, the assumption of normality is increasingly well satisfied, the dispersion matrix Σ for the likelihood function approaches the analytical expression derived above (eq. [11]), and the confidence intervals approach those obtained analytically (fig. 3D). Maximum likelihood fits to simulated data show that there is a small bias in parameter estimates (1%–4%), especially for small data sets (Table 2). The standard deviations of the estimates, 
, decline with increasing n. Observed 
are lower than expected values for small n and very accurate for 
(Table 2; app. A).
and 
, compared with the equivalent model I (black dashed lines). C, Observed wavelet variance for 1,000 simulated processes, with 
, 
, 
m, 
ha (gray lines), their average (black solid line), and their 95% confidence intervals (black dashed line), relative to the solution for model I with 
(red line). D, Log likelihood as a function of the model parameters, when the true value of the dispersal distance (
) is known. Dashed lines indicate the true parameter values: 
, 
.
| n |  |  |  |
|---|---|---|---|
| 25 | 3.12 | .67 | .75 |
| 50 | 3.09 | .59 | .64 |
| 200 | 3.09 | .36 | .39 |
| 500 | 3.06 | .29 | .30 |
| 1,000 | 3.03 | .24 | .25 |
Note. 
and 
are the average and standard deviation of the parameter estimates from 1,000 simulations (
ha, Gaussian kernel 
m). 
is the average estimated standard deviation using 
, where J is the Jacobian and 
, the covariance matrix of the wavelet variances (for details, see app. B).
Model II: Dispersal and Conspecific Density-Dependent Establishment
When spatial density-dependent effects are present, equations (13) do not have a closed analytical solution; an approximate solution can be obtained through moment closure (app. A; Bolker et al. 2000; Murrell et al. 2004):
It is useful to further simplify this model for the purposes of examining its qualitative behavior. If we assume long-range dispersal and competition, then 
, and the scalewise variance becomes
, which reflects the importance of negative density dependence (NDD), and a spatial parameter 
, which reflects the relative scale of competition to dispersal.The influence of the density dependence parameter Π1 on the wavelet variance clearly shows how NDD reduces aggregation relative to what would be expected from dispersal alone (eq. [21]; fig. 4A). When competition and dispersal scales are equal and density-dependent effects equal density-independent ones (
), the distribution is indistinguishable from a CSR process (fig. 4A, green line). When density-dependent mortality is greater than density-independent mortality (
), spatial distributions are disaggregated, because the probability of finding new individuals increases with distance from parents for distances smaller than cD, exactly the pattern posited by Janzen (1970) in his classic article on the influences of natural enemies (fig. 4A, red line).
The influence of parameter Π2 on the wavelet variance shows how the relative scales of dispersal and competition determine the shape of the wavelet variance function (fig. 4B). If competition occurs over shorter distances than dispersal (
), individuals can be disaggregated at small scales, while at larger scales they start to form patterns because of dispersal limitation (fig. 4B, blue line). If competition occurs at longer scales than dispersal (
), the spatial pattern for 
is insensitive to NDD and the solution resembles that for dispersal only, while at larger scales, clustering is reduced far below that expected under dispersal alone (fig. 4B, red line). The latter case results in the formation of a pattern with a characteristic scale (corresponding to the peak of the spectrum), just as occurs in reaction-diffusion systems when a slow diffusing activator (here dispersal) combines with a fast diffusing inhibitor (here competition; Turing 1952).
These qualitative insights are useful, but unfortunately, the underlying equation for the scalewise variance is just an approximation that is valid only for a limited range of parameter values. Numerical simulations have shown that the power-1 closure used to obtain the approximate solution in equations (19) and (20) may fail for moderate overdispersion (Murrell et al. 2004) and high density (Raghib et al. 2011). Moreover, the simplified model (eq. [21]) explored above is limited to long-range dispersal or competition (relatively large cD or cK) or moderate NDD (small Π1). Thus, we next consider the exact behavior and implications for parameter estimation using numerical examples.
This process can be simulated through a small variation on the previous algorithm. Offspring are still generated from a randomly chosen parent and displaced according to the dispersal kernel, but they establish successfully (and thus are retained) only with probability 
. If the offspring fails to establish, another parent is randomly chosen and another offspring randomly displaced.
Simulations confirm the qualitative insights from the analytical approximation hold also for shorter-range dispersal and competition (e.g., fig. 4C). They also demonstrate the feasibility of estimating parameters of model II from wavelet variance data using maximum likelihood, since the likelihood is maximized around the true parameter values (fig. 4D). However, there is a relatively large region of parameter space with very similar likelihoods, even when one of three parameters is assumed known. Estimation of all three parameters simultaneously leads to even larger regions of equifinality and thus wider confidence intervals.
Case Study
We analyzed spatial patterns of four tropical tree species in the 50-ha forest dynamics plot on Barro Colorado Island, Panama (Leigh et al. 1996). All freestanding woody stems >1 cm in diameter were identified to species and mapped to 0.1 m, with censuses in 1982, 1985, and every 5 years through 2010 (Hubbell and Foster 1986; Condit 1998). Barro Colorado data are deposited in the National Center for Ecological Analysis and Synthesis data repository (
Methods
For each species, a raster was created by counting the number of individuals in each 
-m grid cell. The Morlet wavelet variance (
) and the abundance (first and second moment) were calculated for each census for 39 scales log-evenly spaced between 2 and 100 m. We averaged the seven individual census normalized wavelet variances and analyzed this ensemble average as a single realization. (If the spatial patterns in different censuses were statistically independent, the ensemble average would follow a χ2 distribution with degrees of freedom equal to the sum of the degrees of freedom of individual realizations. Turnover rates average 10% per census interval, so the patterns in different censuses are far from independent.)
We first tested whether spatial patterns deviated significantly from random and, if so, at what scales, by comparing observed wavelet variances with the expectation under CSR.
Model I (eq. [17]) and the simplified model II (eq. [21]) were then fitted for each species, using a Laplace kernel for dispersal and competition (Van den Bosch 1988; Bolker and Pacala 1997; Ovaskainen and Cornell 2006a). The models, 
, were fitted using generalized nonlinear least squares (details in app. B). The fitting scheme differed slightly from the numerical examples because the dispersion matrix Σ was not known a priori. Thus, we obtained a first estimate of the parameters using 
. We then simulated 1,000 processes using these estimated parameters, used the simulations to compute the exact matrix Σ, and reestimated the model parameters and their uncertainties. The differences between the first and second estimates were small because the likelihood functions are qualitatively very similar (but the parameter uncertainties are underestimated using the first matrix Σ, especially for small N); thus, no further iterations were necessary. In fitting model II, we set the initial value of 
as equal to the best-fit dispersal parameter from model I, and 
(thus, also 
).
We compared the models in terms of their modified Akaike information criterion (AIC; Burnham and Anderson 2002):
Results
All four species showed highly significant spatial structure, with wavelet variances significantly >1 for spatial scales ≫10 m (fig. 5). Model I reproduced the basic features of the wavelet variance for most species. However, the slope of the observed wavelet variance was generally <2, the value expected under model I; thus, model II was a better fit for all four species (fig. 5). This implies that NDD reduces clustering created by dispersal. Dispersal distance estimates under model II were invariably less than those for model I, and the difference increased as Π1 increased. Observed spatial patterns (left) of four tropical tree species in the 50-ha plot on Barro Colorado Island, Panama, along with their observed (right; dots) and fitted (right; solid lines) wavelet variances. The focal species were chosen to represent the range of patterns observed among the 100 most abundant species. The observed wavelet variances are compared with 95% confidence intervals for a complete spatial randomness process (red dashed lines). Fitting statistics are shown on the right; best-fit parameter values are shown with their standard deviations; all distances are in meters. In all cases, model II provides a significantly better fit than model I (negative ΔAICII-I); the difference is greatest by far in the case of Desmopsis panamensis.
For Desmopsis panamensis, an abundant understory species, the wavelet variance was significantly higher than 1 at scales >10 m but increased less rapidly than predicted by model I. Furthermore, there was disaggregation at small scales, which was correctly captured by model II but not model I (fig. 5). This is the classic case that fairly strong (
) conspecific NDD acts at smaller scales than dispersal (
).
In contrast, Guatteria dumetorum, a less abundant canopy species, appeared less aggregated than D. panamensis. Small scales were entirely bounded by the 95% confidence intervals for a CSR. Both models estimated longer dispersal for Guatteria than for Desmopsis. Though NDD was significant (model II does better), it was much weaker than in Desmopsis (fourfold lower Π1).
Coussarea curvigemmia and Protium tenuifolium, both small trees of moderate abundance, showed the highest values of the wavelet variance, especially at larger scales. Coussarea was estimated to have only weak NDD, while Protium had strong NDD whose effect was particularly obvious at large scales.
Discussion
Ecologists have long examined spatial patterns for insights into ecological processes (e.g., Kershaw 1963). A first focus has been to evaluate whether spatial patterns are significantly nonrandom. The test developed here has two key advantages over previous approaches. First, it is analytical; it does not required Monte Carlo simulations. Second, it identifies at which scales the process departs significantly from CSR. In contrast, randomization tests based on Ripley’s K (e.g., Plotkin et al. 2000) or pair density correlations (e.g., Condit et al. 2000) cannot distinguish whether, for example, departures from CSR at small distances are due to causes acting at small scales or at larger scales (Loosmore and Ford 2006).
Once spatial patterns are established to be nonrandom, the next challenge is to extract information about the processes underlying this nonrandomness. Here, we provide new tools for calculating the likelihood of particular patterns under alternative models and alternative parameter values and, thereby, a method to fit models and choose among them on the basis of static patterns. This represents a major advance over previous approaches because it provides a clear analytical link to models formulated explicitly in terms of ecological processes. Some previous studies have related patterns to classic cluster algorithms (Waagepetersen and Guan 2009; Wiegand et al. 2009), but such linkages provide little insight into, for example, how dispersal and density dependence contribute to pattern formation. Linking to individual-based models derived using basic ecological principles (e.g., reproduction, mortality, dispersal, competition) allows easier interpretation (e.g., Anand and Langille 2010). Further, our methods enable quantification of uncertainties in resulting parameter estimates, critical information for comparing and integrating results of multiple studies. An important caveat, as mentioned above, is that the approximation used in the moment equation is not uniquely defined and the error may be dependent on the model parameter (Murrell et al. 2004), a clear drawback for inferential analysis that should be more carefully addressed in future work.
A major advantage of linking ecological patterns with individual-based, spatially explicit models is that all model assumptions are clearly stated in terms of ecological processes, and these can easily be kept in mind when results are interpreted. Taking the above case study results as an example, we noted that dispersal distance estimates under model II are invariably less than those for model I, and the difference increases as Π1 increases. This makes sense, since the dispersal estimates from model I are best thought of as estimates of effective dispersal distance—that is, the dispersal distance of successfully establishing individuals—which becomes ever larger than actual dispersal distance as conspecific NDD becomes stronger. Even for model II, it is important to note that estimated dispersal distances are not expected to correspond exactly to real dispersal distances, because the fitted model assumes all individuals are reproductive. Given that in reality only the larger individuals are reproductive, real dispersal distances must necessarily be longer than those estimated under this model. Similarly, estimates of NDD obtained from fitting model II may be biased because in the real world competition for resources is size dependent, while the fitted model lacks size structure (for an age-structured model, see Murrell 2009). Interpretation must also consider the suite of alternative models examined and unexamined. Models that incorporate NDD in survival, growth, or fecundity are likely to produce similar patterns to model II, which has NDD establishment.
The models analyzed here do not include habitat heterogeneity, which affects spatial structure in many populations. Where habitat variability operates scales much larger than the scales of dispersal and competition, the spatial patterns at these smaller scales would not be affected, and to a first-order approximation, habitat effects can be neglected in estimating dispersal and competition parameters. However, in other cases, habitat heterogeneity may act on scales comparable to dispersal and competition, in which case the spatial pattern will depend on all three processes (habitat association, dispersal, and competition) and it will be difficult to disentangle their influences. The Barro Colorado Island plot can be divided into five broad habitat types: slope, low plateau, high plateau, streamside, and swamp (Harms et al. 2001). Of our focal species, Desmopsis panamensis is positively associated with high plateau, Guatteria dumetorum is positively associated with slope and negatively with high plateau, and Protium tenuifolium is negatively associated with the swamp, according to conservative habitat randomization tests (Harms et al. 2001). Given the large spatial scale of these habitats, they are expected to affect wavelet variances at scales of ∼100 m and greater. In contrast, dispersal and competition operate mostly on smaller scales in our system. Thus, it is not surprising that we are able to obtain a good fit to small-scale spatial structure in these species with a model that ignores habitat heterogeneity.
The models here also ignore heterospecific interactions as determinants of population spatial structure, despite the fact that heterospecific neighbors also have strong influences on the growth and survival of individuals (Uriarte et al. 2004; Comita et al. 2010). A multispecies model that includes heterospecific as well as conspecific interactions has been developed by Law and Dieckmann (2000). We can use this model to simultaneously fit the wavelet variance of two individual species and their covariance. However, in species-rich communities such as Barro Colorado, with ∼300 tree species in 50 ha, a dilution effect yields approximate species independence (Wiegand et al. 2012). As a result, there is little improvement, if any, in fitting the multivariate model compared with fitting each species independently. Clearly, there are other communities in which heterospecific interactions are relatively more important in shaping spatial patterns and in which a multispecies model would thus provide a better fit and yield more information about the processes driving dynamics in the system.
Our methods have certain limitations, which are general to any analyses of spatial patterns. First, multiple processes can produce similar spatial patterns, providing limits to what any analysis of static patterns alone can yield. For example, strong NDD can counterbalance dispersal limitation effects, producing a pattern that is similar or in extreme cases indistinguishable from a random pattern (long dispersal). In general, this problem arises when the scales of different ecological processes overlap, making it hard to quantify, for these scales, how much variance is produced (or reduced) by one process or another. Second, analyses of processes with larger scales require data with adequate spatial extension. Thus, for example, estimates of dispersal distance that are long relative to the scale of the data (
) are less accurate because asymptotic trends at larger scales are not well resolved or are absent in the scalewise variances.
Recommendations and Future Directions
Here, we developed new tools for investigating the effects of ecological processes on spatial point patterns and deriving information on ecological processes from these patterns. Future theoretical work should investigate how habitat heterogeneity is expected to influence scalewise variances (Ovaskainen and Cornell 2006a) and evaluate the potential for linking species spatial patterns to habitat. Exploration of expected scalewise variances for additional dispersal and interaction kernels, anisotropic effects, different types of NDD (and potentially positive density dependence), and more complex models incorporating age and size structure (Murrell 2009) and/or multiple species (Law and Dieckmann 2000) could further shed light on the impacts of different ecological processes on spatial patterns. These methods can also be extended beyond sessile organisms by considering adding a movement kernel to the model, following Law and Dieckman (2000).
Maximum insight into underlying ecological processes will be gained by judiciously combining the spatial analyses pioneered here with other sources of information. Fitting even moderately complex models to spatial data alone results in wide confidence intervals (e.g., fig. 4) because the scales of many ecological processes often overlap (e.g., short dispersal and NDD), with the consequence that similar patterns are produced with different combinations of the parameters. The likelihood function for the normalized wavelet variances can be combined with prior information or analyses of other data sets to narrow confidence intervals. For example, a study of seed dispersal might inform prior distributions on the dispersal parameter. Data from multiple censuses may provide further information about the spatial-temporal dynamics (Ovaskainen and Cornell 2006b), enabling, for example, independent estimates of fecundity and mortality. Wide application of the methods developed here could shed light on the prevalence and strength of NDD, the relative scales of dispersal and competition, and the distribution of effective dispersal distances within and among populations and communities.
Acknowledgments
We gratefully acknowledge the financial support of the Smithsonian Institution Global Earth Observatories, and we thank R. Condit and S. Hubbell for providing the Barro Colorado Island tree data and J. Calabrese and R. Chrislom for useful comments. The Barro Colorado Island forest dynamics research project was made possible by National Science Foundation grants to S. P. Hubbell (DEB-0640386, DEB-0425651, DEB-0346488, DEB-0129874, DEB-00753102, DEB-9909347, DEB-9615226, DEB-9615226, DEB-9405933, DEB-9221033, DEB-9100058, DEB-8906869, DEB-8605042, DEB-8206992, and DEB-7922197), support from the Center for Tropical Forest Science, the Smithsonian Tropical Research Institute, the John D. and Catherine T. MacArthur Foundation, the Mellon Foundation, and the Small World Institute Fund.
Literature Cited
Abramowitz, M., and I. A. Stegun. 1965. Handbook of mathematical functions. Dover, New York. Anand, M., and A. Langille. 2010. A model-based method for estimating effective dispersal distance in tropical plant populations. Theoretical Population Biology 77:219–226. Bagchi, R., P. Henrys, P. Brown, and D. Burslem. 2011. Spatial patterns reveal negative density dependence and habitat associations in tropical trees. Ecology 92:1723–1729. Bartlett, M. 1964. The spectral analysis of two-dimensional point processes. Biometrika 51:299–311. Bolker, B., and S. Pacala. 1997. Using moment equations to understand stochastically driven spatial pattern formation in ecological systems. Theoretical Population Biology 52:179–197. Bolker, B. M., S. W. Pacala, and S. A. Levin. 2000. Moment methods for ecological processes in continuous space. Pages 388–411 in U. Dieckmann, R. Law, and J. A. J. Metz, eds. The geometry of ecological interactions. Cambridge University Press, Cambridge. Bolker, B. M., S. W. Pacala, and C. Neuhauser. 2003. Spatial dynamics in model plant communities: what do we really know? American Naturalist 162:135–148. Burnham, K. P., and D. R. Anderson. 2002. Model selection and multimodel inference: a practical information-theoretic approach. 2nd ed. Springer, New York. Cazelles, B., M. Chavez, and D. Berteaux. 2008. Wavelet analysis of ecological time series. Oecologia (Berlin) 156:287–304. Comita, L. S., H. C. Muller-Landau, S. Aguilar, and S. P. Hubbell. 2010. Asymmetric density dependence shapes species abundance in a tropical tree community. Science 329:330–332. Condit, R. 1998. Tropical forest census plots. Springer, Berlin. Condit, R., P. S. Ashton, P. Baker, S. Bunyavejchewin, S. Gunatilleke, N. Gunatilleke, S. P. Hubbell, et al. 2000. Spatial patterns in the distribution of tropical tree species. Science 288:1414–1418. Dale, M., and M. Mah. 1998. The use of wavelets for spatial pattern analysis in ecology. Journal of Vegetation Science 9:805–814. Daley, D., and D. Vere-Jones. 2003. An introduction to the theory of point processes. I. Elementary theory and methods. 2nd ed. Springer, New York. Daubechies, I. 1992. Ten lectures on wavelets. Society for Industrial and Applied Mathematics, Philadelphia. Detto, M., A. Molini, G. Katul, P. Stoy, S. Palmroth, and D. Baldocchi. 2012. Causality and persistence in ecological systems: a nonparametric spectral Granger causality approach. American Naturalist 179:524–535. Dieckmann, U., and R. Law. 2000. Relaxation projections and the method of moments. Pages 412–455 in U. Dieckmann, R. Law, and J. A. J. Metz, eds. The geometry of ecological interactions. Cambridge University Press, Cambridge. Durrett, R., and S. Levin. 1994. The importance of being discrete (and spatial). Theoretical Population Biology 46:363–394. Harms, K. E., R. Condit, S. P. Hubbell, and R. B. Foster. 2001. Habitat associations of trees and shrubs in a 50-ha Neotropical forest plot. Journal of Ecology 89:947–959. Hartig, F., J. M. Calabrese, B. Reineking, T. Wiegand, and A. Huth. 2011. Statistical inference for stochastic simulation models: theory and application. Ecology Letters 14:816–827. Hubbell, S. P., and R. B. Foster. 1986. Biology, chance, and history and the structure of tropical rain forest tree communities. Pages 314–329 in J. Diamond and T. J. Case, eds. Community ecology. Harper & Row, New York. Janzen, D. H. 1970. Herbivores and the number of tree species in tropical forests. American Naturalist 104:501–528. Keitt, T. H., and D. L. Urban. 2005. Scale-specific inference using wavelets. Ecology 86:2497–2504. Kershaw, K. A. 1963. Pattern in vegetation and its causality. Ecology 44:377–388. Kumar, P., and E. Foufoula-Georgiou. 1997. Wavelet analysis for geophysical applications. Reviews of Geophysics 35:385–412. Law, R., and U. Dieckmann. 2000. A dynamical system for neighborhoods in plant communities. Ecology 81:2137–2148. Law, R., J. Illian, D. F. R. P. Burslem, G. Gratzer, C. V. S. Gunatilleke, and I. A. U. N. Gunatilleke. 2009. Ecological information from spatial patterns of plants: insights from point process theory. Journal of Ecology 97:616–628. Leigh, E. G., Jr., A. S. Rand, and D. M. Windsor. 1996. The ecology of a tropical forest: seasonal rhythms and long-term changes. Smithsonian Institution, Washington, DC. Levin, S. A. 2000. The problem of pattern and scale in ecology. Ecology 73:1943–1967. Loosmore, N. B., and E. D. Ford. 2006. Statistical inference using the G or K point pattern spatial statistics. Ecology 87:1925–1931. McIntire, E. J. B., and A. Fajardo. 2009. Beyond description: the active and effective way to infer processes from spatial patterns. Ecology 90:46–56. Meyers, S. D., B. Kelly, and J. O’Brien. 1993. An introduction to wavelet analysis in oceanography and meteorology: with application to the dispersion of Yanai waves. Monthly Weather Review 121:2858–2866. Mi, X., H. Ren, Z. Ouyang, W. Wei, and K. Ma. 2005. The use of the Mexican Hat and the Morlet wavelets for detection of ecological patterns. Plant Ecology 179:1–19. Murrell, D. J. 2009. On the emergent spatial structure of size-structured populations: when does self-thinning lead to a reduction in clustering? Journal of Ecology 97:256–266. Murrell, D. J., U. Dieckmann, and R. Law. 2004. On moment closures for population dynamics in continuous space. Journal of Theoretical Biology 229:421–432. Murrell, D. J., and R. Law. 2002. Heteromyopia and the spatial coexistence of similar competitors. Ecology Letters 6:48–59. Ovaskainen, O., and S. J. Cornell. 2006a. Asymptotically exact analysis of stochastic metapopulation dynamics with explicit spatial structure. Theoretical Population Biology 69:13–33. ———. 2006b. Space and stochasticity in population dynamics. Proceedings of the National Academy of Sciences of the USA 103:12781–12786. Pacala, S. 1986. Neighborhood models of plant population dynamics. 2. Multi-species models of annuals. Theoretical Population Biology 29:262–292. Percival, D. P. 1995. On estimation of the wavelet variance. Biometrika 82:619–631. Platt, T., and K. L. Denman. 1975. Spectral analysis in ecology. Annual Review of Ecology and Systematics 6:189–210. Plotkin, J. B., M. D. Potts, N. Leslie, N. Manokaran, J. Lafrankie, and P. S. Ashton. 2000. Species-area curves, spatial aggregation, and habitat specialization in tropical forests. Journal of Theoretical Biology 207:81–99. Raghib, M., N. A. Hill, and U. Dieckmann. 2011. A multiscale maximum entropy moment closure for locally regulated space-time point process models of population dynamics. Journal of Mathematical Biology 62:605–653. Serroukh, A., A. T. Walden, and D. B. Percival. 2000. Statistical properties and uses of the wavelet variance estimator for the scale analysis of time series. Journal of the American Statistical Association 95:184–196. Thompson, S. E., and G. G. Katul. 2011. Inferring ecosystem parameters from observation of vegetation patterns. Geophysical Research Letters 38:2–5. Torrence, C., and G. P. Compo. 1998. A practical guide to wavelet analysis. Bulletin of the American Meteorological Society 79:61–78. Turing, A. 1952. The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society B: Biological Sciences 327:37–72. Uriarte, M., R. Condit, C. D. Canham, and S. Hubbell. 2004. A spatially explicit model of sapling growth in a tropical forest: does the identity of neighbours matter? Journal of Ecology 92:348–360. Van den Bosch, F., J. Zadoks, and J. Metz. 1988. Focus expansion in plant disease. II. Realistic parameter-sparse models. Phytopathology 78:59–64. Waagepetersen, R., and Y. Guan. 2009. Two-step estimation for inhomogeneous spatial point processes. Journal of the Royal Statistical Society B 71:685–702. Wiegand, T., A. Huth, S. Getzin, X. Wang, Z. Hao, C. V. S. Gunatilleke, and I. A. U. N. Gunatilleke. 2012. Testing the independent species’ arrangement assertion made by theories of stochastic geometry of biodiversity. Proceedings of the Royal Society B: Biological Sciences 279:3312–3320. Wiegand, T., I. Martínez, and A. Huth. 2009. Recruitment in tropical tree species: revealing complex spatial patterns. American Naturalist 174:E106–E140.
