The attached dataset shows approximately 6000 saplings in approximately 50 variable sized forest gaps. I am interested in learning how these saplings are growing within their respective gaps (i.e. clustered, random, dispersed). As you know, a traditional approach would be to run Global Moran's I. However, aggregations of trees within aggregations of gaps seems to be an inappropriate use of Moran's I. I ran some test statistics with Moran's I using a threshold distance of 50 meters, which produced nonsensical results (i.e. p-value = 0.0000000...). The interaction among the gap aggregations are likely producing these results. I have considered creating a script to loop through individual canopy gaps and determine the clustering within each gap, although displaying these results to the public would be problematic.

What is the best approach to quantify clustering within clusters?

enter image description here

  • 1
    Aaron, you say that you tried running Moran's I, are you interested in measuring how a sapling's attribute compares with neighboring saplings' attributes (i.e. are you dealing with a marked point pattern)? The title seems to imply that you are only interested in the location of the saplings relative to one another and not their attributes.
    – MannyG
    Commented Nov 28, 2012 at 13:39
  • @MannyG Yes, I am only interested in determining if the saplings are clustered relative to the locations of other saplings within any given forest gap. There is only one species of interest and the size of the saplings are not of interest.
    – Aaron
    Commented Nov 28, 2012 at 14:01

3 Answers 3


You do not have a uniform random field, so attempting to analyze all of your data at once will violate the assumptions of any statistic you choose to throw at the problem. It is unclear from your post if your data is a marked point process (i.e, diameter or height associated with each tree location). If this data is not representing a marked point process I have no idea how you applied a Moran's-I. If the data only represents spatial locations, I would recommend using a Ripley's-K with the Besag-L transformation to standardize the null expectation at zero. This allows for a multiscale assessment of clustering. If your data has an associated value, then your best option is a local Moran's-I (LISA). I would actually look at it with both statistics. Regardless of your choice, you will still need to loop through each individual site to produce valid results. Here is some example R code for a Monte Carlo simulation of Ripley's-K/Besag's-L using the built-in redwood sapling dataset. It should be fairly straightforward to modify this to loop through your sites and produce a graph for each one.


spp <- SpatialPoints(coords(redwood))


# COERCE TO spatstat ppp OBJECT
spp.ppp=as.ppp(coordinates(spp), W)                     

area <- area.owin(W)
lambda <- spp.ppp$n/area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
bw <- seq(0, rmax, by=rmax/10)  

Lenv <- envelope(spp.ppp, fun="Kest", r=bw, i="1", j="2", nsim=99, nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE)            
plot(Lenv, main="Besag's-L", xlab="Distance", ylab="L(r)", legend=F, col=c("white","black","grey","grey"), 
    lty=c(1,2,2,2), lwd=c(2,1,1,1) )
     polygon( c(Lenv$r, rev(Lenv$r)), c(Lenv$lo, rev(Lenv$hi)), col="lightgrey", border="grey")
       lines(supsmu(bw, Lenv$obs), lwd=2)
       lines(bw, Lenv$theo, lwd=1, lty=2)
         legend("topleft", c(expression(hat(L)(r)), "Simulation Envelope", "theo"), pch=c(-32,22),
                col=c("black","grey"), lty=c(1,0,2), lwd=c(2,0,2), pt.bg=c("white","grey"))
  • 1
    But you can't just go using the convex hull as the window for your point pattern! Remember, the window is the area in which the pattern that produces the points operates. You know a-priori that the trees only grow in these set regions, and you have to set your window to reflect that. You might mitigate this by setting the range of K(r) to something very small, of the order of 0.3x size of your clearings, but you'll get biased results because of the lack of edge effect corrections. Jeffrey is using the size of the whole study area to define his rmax.
    – Spacedman
    Commented Nov 29, 2012 at 8:32
  • 1
    In my example, Yes I am using the entire region. By this is exactly why I recommended looping through each sample site (gap). Each time you subset to a specific sample area you would rerun the analysis. You cannot treat the entire study area as your random field because you do not have continuous sampling. Having only sampled gaps you, in effect, have independent plots. The Kest function I am calling, by default, uses a "border" edge correction. There are other edge correction options available. I would argue that your experimental unit is the canopy gap and should be analyzed as such. Commented Nov 29, 2012 at 15:51
  • 1
    In thinking about this a bit more. You should really be using polygons that represent each gap as your window. If you subset your problem to reflect the experimental unit then the CSR and K will be biased because the area does not reflect the actually canopy gap size. This is an issue in both my and @Spacedman's recommendations. Commented Nov 29, 2012 at 16:04
  • 2
    Note my extended example only used a coarse grid because it was a fairly simple way of creating something with roughly the right structure. Your mask should look like a map of your open forest areas. Its technically wrong to try and define the mask from the data!
    – Spacedman
    Commented Nov 29, 2012 at 16:12
  • 1
    @Spacedman. I like your approach and it is certainly efficient. My specific concern is that the canopy gaps are the experimental unit. In your approach, if two gaps are proximal the bandwidth could plausibly include observations from a different sampling units. Additionally, the resulting statistic should not reflect the "pool" of experimental units but should be representative of each unit and inference on spatial process drawn from common patterns across the experimental units. If treated globally, it represents an nonstationary intensity process which, violates statistical assumptions. Commented Nov 29, 2012 at 18:10

What you have is a point pattern with a window that is a number of small disconnected polygonal regions.

You should be able to use any of the tests in package:spatstat for CSR as long as you feed it with a correct window. This can be either a number of sets of (x,y) pairs defining each clearing or a binary matrix of (0,1) values over the space.

First lets define something that looks a bit like your data:

nclust <- function(x0, y0, radius, n) {
               return(runifdisc(n, radius, centre=c(x0, y0)))
c = rPoissonCluster(15, 0.04, nclust, radius=0.02, n=5)

and lets pretend our clearings are square cells that just happen to be this:

m = matrix(0,20,20)
imask = owin(c(0,1),c(0,1),mask = t(m)==1 )
pp1 = ppp(x=c$x,y=c$y,window=imask)

So we can plot the K-function of those points in that window. We expect this to be non-CSR because the points seem clustered within the cells. Notice I have to change the range of distances to be small - of the order of the cell size - otherwise the K-function gets evaluated over distances the size of the whole pattern.


If we generate some CSR points in the same cells, we can compare the K-function plots. This one should be more like CSR:

ppSim = rpoispp(73/(24/400),win=imask)

two point patterns in patchy windows

You can't really see the points clustered in the cells in the first pattern, but if you plot it on its own in a graphics window its clear. The points in the second pattern are uniform within the cells (and do not exist in the black region) and the K-function is clearly different from Kpois(r), the CSR K-function for the clustered data and similar for the uniform data.


In addition to Andy's post:

What you want to calculate is a measure of spatial homogeneity (ergo the hypothesis: "Are your points clustered?") such as Ripley's L and K function.

This blog post explains the how-to in R pretty well. Based on the described code i would first label each cluster in your dataset and then calculate in a loop for each cluster the critical envelope through Ripley's K

  • I've currently deleted my answer. Some brief analysis suggested that opportunistically identifying the plots based on K-means biased the local statistics to be more clustered than would be suggested by chance. This answer still applies though +1, (just the making of the windows from the data is more problematic than my original answer would suggest).
    – Andy W
    Commented Nov 28, 2012 at 17:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.