LGCPs - Plot sampling
David Borchers and Finn Lindgren
Generated on 2023-02-02Source:
This practical demonstrates use of the
lgcp, which you need to use when you have observed
points from only a sample of plots in the survey region.
data(gorillas, package = "inlabru")
This dataset is a list (see
help(gorillas) for details.
Extract the the objects you need from the list, for convenience:
nests <- gorillas$nests mesh <- gorillas$mesh boundary <- gorillas$boundary gcov <- gorillas$gcov
gorillas data also contains a plot sample subset
which covers 60% of the survey region.
sample <- gorillas$plotsample
On this plot survey, only points within the rectangles are detected, but it is also informative to plot all the points here (which if it was a real plot survey you could not do, because you would not have seen them all).
plotwithall <- ggplot() + gg(boundary) + gg(sample$plots) + gg(nests, pch = "+", cex = 4, color = "blue") + geom_text(aes(label = sample$counts$count, x = sample$counts$x, y = sample$counts$y)) + gg(sample$nests, pch = "+", cex = 4, color = "red") + coord_fixed() + labs(x = "Easting", y = "Northing") plot(plotwithall)
The observed nest locations are in the SpatialPointsDataFrame
sample$nests, and the plots are in the
sample$plots. Again, we are using
the following SPDE setup:
Fit an LGCP model with SPDE only to these data by using the
samplers= argument of the function
Plot the density surface from your fitted model
lambda.sample.plot <- ggplot() + gg(lambda.sample) + gg(sample$plots) + gg(boundary, col = "yellow") + coord_fixed() lambda.sample.plot
Estimate the integrated intensity lambda. We compute both the overall integrated intensity, representative of an imagined new realisation of the point process, and the conditional expectation that takes the actually observed nests into account, by recognising that we have complete information in the surveyed plots.
Lambda <- predict(fit, ipoints(boundary, mesh), ~ sum(weight * exp(my.spde + Intercept))) Lambda.empirical <- predict( fit, rbind( cbind(ipoints(boundary, mesh), data.frame(all = TRUE)), cbind(ipoints(sample$plots, mesh), data.frame(all = FALSE)) ), ~ (sum(weight * exp(my.spde + Intercept) * all) - sum(weight * exp(my.spde + Intercept) * !all) + nrow(sample$nests)) ) rbind( Lambda, Lambda.empirical )
Fit the same model to the full dataset (the points in
gorillas$nests), or get your previous fit, if you kept it.
Plot the intensity surface and estimate the integrated intensity
Your plot should look like this:
Lambda.all should be close to each other if the plot
samples gave sufficient information for the overall prediction:
rbind( Lambda, Lambda.empirical, Lambda.all, Lambda.all.empirical = c(nrow(gorillas$nests), 0, rep(nrow(gorillas$nests), 3), rep(NA, 4)) ) #> Warning in rbind(deparse.level, ...): number of columns of result, 8, is not a #> multiple of vector length 9 of arg 4 #> mean sd q0.025 q0.5 q0.975 median mean.mc_std_err #> 1 656.8203 52.07817 571.6200 654.5181 777.4664 654.5181 5.207817 #> 2 651.1208 36.27578 587.0249 649.5293 720.0303 649.5293 3.627578 #> 3 673.2638 25.86937 624.9822 673.6706 725.9925 673.6706 2.586937 #> 4 647.0000 0.00000 647.0000 647.0000 647.0000 NA NA #> sd.mc_std_err #> 1 4.161586 #> 2 2.973416 #> 3 1.704062 #> 4 NA
Now, let’s compare the results
library(patchwork) lambda.sample.plot + lambda.all.plot + plot_layout(guides = "collect") & theme(legend.position = "left") & scale_fill_continuous(limits = range(c(0, 340)))
Do you understand the reason for the differences in the posteriors of the abundance estimates?