spatial statistics - Determining if trees within forest gaps are clustered using R?

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?

Answer

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.

# ADD REQUIRED PACKAGES
require(sp)
require(spatstat)
options(scipen=5)

# USE REDWOOD SAPLING DATASET
spp <- SpatialPoints(coords(redwood))


###################################################
###### START BESAG'S-L MONTE CARLO  ANALYSUS ######
###################################################
# CREATE CONVEX HULL FOR ANALYSIS WINDOW                       
W=ripras(coordinates(spp)) 

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


# ESTIMATE BANDWIDTH
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)  

# CALCULATE PERMUTED CROSS-K AND PLOT RESULTS       

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"))