Here is some basic code to get started with spatstat using the Beilschmiedia pendula dataset.
install.packages("spatstat")
library(spatstat)
data(bei)
plot(bei)
plot(bei.extra$elev)
plot(bei.extra$grad)
fit.1 = ppm(bei, ~poly(elev, grad, degree = 2), covariates = bei.extra)
plot(predict(fit.1))
plot(bei, add = TRUE, cex = 0.4)
diagnose.ppm(fit.1)
AIC(fit.1)
fit.2 = ppm(bei, ~poly(elev, grad, degree = 2) + poly(x, y, degree = 2), covariates = bei.extra)
plot(predict(fit.2))
plot(bei, add = TRUE, cex = 0.4)
diagnose.ppm(fit.2)
AIC(fit.2)
Can you find a model which removes the main trend in residual plots?
To view the tree locations atop of a map of gradient, type:
ReplyDeleteplot(bei.extra$grad)
plot(bei, add = TRUE, cex = 0.4)
To view the tree locations atop of a map of elevation, type:
ReplyDeleteplot(bei.extra$elev)
plot(bei, add = TRUE, cex = 0.4)
To determine significance of model parameters:
ReplyDeletesummary(fit.1)
fit.2 still had some problems for me but it is a lot better than fit.1.
ReplyDeleteYou can compare them using
anova(fit.1,fit.2,test="Chisq")
if you believe the model is plausible. Big if though.
I reckon a key thing we are missing here is aspect - eyeballing the raw data maps this species appears to be found on southerly aspects as well as at high elevation?
The other thing is point interactions (dispersal etc...)
It looks like aspect might be a useful covariate. Loving the residual plot functionality spatstat!
ReplyDeleteTo test for dependence among point locations:
ReplyDeleteenv.1 = envelope(fit.1)
env.2 = envelope(fit.2)
I'm going to play devil's advocate and say why didn't I just generate some random "absence" points (since we have censused all the presences, and any point that isn't a presence is by definition an absence) and do a t-test against elevation???
ReplyDeleteOkay a bit slow on getting the comment up... but here's the AICs for both fits:
ReplyDeletefit.1 = fit.1 = ppm(bei, ~poly(elev, grad, degree = 2), covariates = bei.extra)
AIC = 41563.16
fit.2 = ppm(bei, ~poly(elev, grad, degree = 2) + poly(x, y, degree = 2), covariates = bei.extra)
AIC = 40739.85
A t-test might be useful for determining whether elevation is important, but how can it be used to make predictions of species intensity based on elevation? A point process model takes into account the spatial nature of the data such that maps of predicted intensities, etc., can be made in addition to testing the significance of covariates.
ReplyDeleteAdditionally, analyzing the data as a point pattern does not require a subjective subset of absences, and further allows users to access tools for assessing goodness of fit.
ReplyDelete