# 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? • 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 Nov 28 '12 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 Nov 28 '12 at 14:01

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"))
``````
• 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 Nov 29 '12 at 8:32
• 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. – Jeffrey Evans Nov 29 '12 at 15:51
• 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. – Jeffrey Evans Nov 29 '12 at 16:04
• 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 Nov 29 '12 at 16:12
• @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. – Jeffrey Evans Nov 29 '12 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:

``````set.seed(310366)
nclust <- function(x0, y0, radius, n) {
}
c = rPoissonCluster(15, 0.04, nclust, radius=0.02, n=5)
plot(c)
``````

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

``````m = matrix(0,20,20)
m[1+20*cbind(c\$x,c\$y)]=1
plot(pp1)
``````

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.

``````plot(Kest(pp1,r=seq(0,.02,len=20)))
``````

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)
plot(ppSim)
plot(Kest(ppSim,r=seq(0,.02,len=20)))
`````` 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.