Step One:

The function nodeFilling takes a polygon representing, in this case, a lake with a causeway and no islands (“holes”") and fills the region with a grid of equally-spaced nodes. You should select node_spacing to be small, bound only by the speed of processing the data. With small node_spacing, this analysis closely approximates a diffusion.

plot.new()
data(polygon1)
head(polygon1)
#>           [,1]      [,2]
#> [1,] 0.6421053 0.8132050
#> [2,] 0.6845247 0.4814305
#> [3,] 0.7057345 0.2858322
#> [4,] 0.6696779 0.2066025
#> [5,] 0.5190888 0.1892710
#> [6,] 0.5445405 0.4145805
nodeFillingOutput <- nodeFilling(poly=polygon1, node_spacing=0.02)
plot(nodeFillingOutput)

Step Two:

Now for the fun part. The function formLattice will ‘connect the dots’ to create a neighbor structure on the nodes. Neighbors are defined to be locations to the N, NE, E, SE, S, SW, W, and NW of each node (as long as there IS a node in that direction). You can examine a plot of the resulting object to see if all regions that should be connected are and that nodes, for instance, across a causeway are not connected.

formLatticeOutput <- formLattice(nodeFillingOutput)
plot(formLatticeOutput)

Optional Step:

If you see links that you would like to remove (for instance, a link crosses a boundary), or see links that you would like to add, use the function editLattice.

Step Three:

Next, we need to add locations for the point process (e.g. animal relocations) to the lattice object. After the data has been added, we can select a value that represents the number of steps the random walk takes around each observation. These random walks essentially produce a kernel around the observation. There is another argument that we won’t use at this time, . This is the maximum probability that the random walk moves in any single step and is set, by default, to 0.5. For instance, if a node has eight neighbors (the largest possible) and M = 0.5, the flow rate to each neighbor is 0.5/8 = 0.0625. Decreasing this value will decrease the width of a kernel for a fixed step size, so generally it should be decreased when the optimal step size turns out to be so small that large parts of the region are outside any kernel.

Internally, this package uses the function addObservations, but this is called by the functions that we will use, which are createDensity which produces a predicted density surface for a given number of steps , and crossvalDensity which uses crossvalidation over a range of value of to select the optimal step size.

Just for this example I’ll use the function csr from the package splancs to simulate a completely random point process with all observations to the east of the causeway removed. Then I’ll feed that into crossvalDensity:

Pointdata <- splancs::csr(polygon1,150)
colnames(Pointdata) <- c("x","y")
Pointdata <- Pointdata[Pointdata[,1]<0.5,]
full_polygon <- rbind(polygon1,polygon1[1,])
colnames(full_polygon) <- c("x","y")
plot(full_polygon,type="l",plt=c(0,1,0,1))
title("Simulated point process")
points(Pointdata,pch=19)

out <- crossvalDensity(formLatticeOutput,PointPattern=Pointdata, 
  M=0.5,max_steps = 150)
plot(1:150,out$ucv,type="l")

out$k
#> [1] 150

The minimum uvc (Unbiased CrossValidation criterion of Sain, Baggerly and Scott (1994)) was at 150 steps. I’ll use that many steps to get a final, smooth density estimate using createDensity.

densityOut <- createDensity(formLatticeOutput,
  PointPattern=Pointdata, k=out$k,intensity=FALSE, sparse = TRUE)
plot(densityOut)

For the ecologically-inclined, the function homerange computes the area of the 95% confidence region.

homerange(densityOut, percent=0.95)

#> Proportion of region in homerange =  0.6678922 
#> Area of homerange =  0.2166476

x_poly = c(0, 0, 0.2, 0.2, 1, 1, 1.2, 1.2, 2, 2, 0)
y_poly = c(0, 1, 1, 0.2, 0.2, 1, 1, 0.2, 0.2, 0, 0)
polyg <- cbind(x_poly, y_poly)
nodeFillingOutput <- nodeFilling(polyg, node_spacing=0.15)
plot(nodeFillingOutput)

plot(nodeFillingOutput)
text(nodeFillingOutput$nodes+rep(c(0.1,0),each=23),col=2,labels=letters[1:23])

flo <- formLattice(nodeFillingOutput)
plot(flo)

PD = rbind(c(1.5,0.1),c(1.6,0.1),c(1.1,0.5))
out <- createDensity(flo,PointPattern=PD, 
  M=0.5, k=1)
plot(out)
text(out$nodes, labels=round(out$probs,4),cex=0.6)

out <- createDensity(flo,PointPattern=PD, 
  M=0.5, k=600)
plot(out)
text(out$nodes, labels=round(out$probs,3),cex=0.6)

Error map

The output of crossvalNparReg also includes an estiamate of the underlying error \(\hat{\sigma}^2 = \sum_{i=1}^n{(y_i-\hat{y}_{i(i)})^2/n}\), which is the mean of the squared deleted residuals at all observations. In this case it equals 0.4466022.

You can also plot a countour map of spatial standard errors based on the Nadaraya-Watson kernel regression variance estimator. To do this, use the functions varianceMap and plot.varianceMapOut. Note that the standard error is not necessarily a good measure of mean squared prediction error, as the kernel estimator (especially when extrapolating in a part of the region far from the observations) can be biased.

varianceMapOut <- varianceMap(formLatticeOutput,Z,PointPattern=grid2,k=20)

plot(varianceMapOut)

Conversion to UTM

Shapefiles (and boundaries in general) are often in latitude and longitude. However, because the distance corresponding to one degree of latitude is usually not the same as the distance corresponding to one degree of longitude, the diffusion modeled by latticeDensity isn’t isotropic. To solve this problem, you can convert to UTM (Universal Transverse Mercator) coordinates. One R package that can do this conversion is rgdal:

library(rgdal)
xy <- cbind(c(118, 118.3, 118.7, 119), c(10, 25, 48, 49))
now_in_UTM <- project(xy, proj = "+proj=utm +zone=51 ellps=WGS84")
now_in_UTM
#>           [,1]    [,2]
#> [1,] -48636.65 1109577
#> [2,]  25379.91 2773187
#> [3,] 179270.67 5325252
#> [4,] 207462.87 5435168

Holes

Islands inside of lakes and other ‘holes’ are easily accomodated in latticeDensity. Each polygonal hold is stored in a different object, then the corresponding nodes are cut out of the irregular region in the nodeFilling function.

boundary <- cbind(c(0,0,0.2,1,1),c(0,0.8,1,1,0))
hole1 <- cbind(c(0.1,0.3,0.4),c(0.2,0.3,0.8))
hole2 <- cbind(c(0.8,0.9,0.9),c(0.6,0.7,0.95))
region <- nodeFilling(poly = boundary, node_spacing = 0.03, hole_list = list(hole1, hole2))
#> [1] "This function does not check to see if the holes"
#> [1] "are nonintersecting, or whether they are contained"
#> [1] "inside the boundary"
plot(region)

lattice_output = formLattice(region)
plot(lattice_output)

There are some links that are crossing the islands, but the function editLattice can be used to remove these.

NOTE: Warnings

Occasionally the message

validspamobject() is deprecated. Use validate_spam() directly

will occur. This is generated by the spam package or packages that call spam and does not seem to have an effect on the output of latticeDensity.

Also, expect that some functions (crossvalidation or deletedResid, for instance) will not run properly if you only have one or two observations.

Finally, whenever a location is outside the support of all kernels, it is given the estimated value equal to the mean of the data. Usually it is better to make the kernels cover most of the map. Usually k is reasonably large compared to the diameter of the grid, and you are fine, otherwise you might consider decreasing the movement rate by decreasing M somewhat before crossvalidating. This will tend to make the kernels have larger support.

In a similar vein, crossvalidation often seems to undersmooth data when the observations are in clusters, a common problem with spatial crossvalidation in any context. If you think the optimum number of steps (k) under crossvalidation is too low, just use a larger value of k and examine the deleted sums of squared errors to see that they don’t increase too much.