When fitting a Poisson process model on multiple point patters using the spatstat R package (mppm function), is there a computationally efficient way to extract the fitted expected number of points from the mppm fit? I expand with a toy example below:
Let's say that my data are point patterns over T time periods and I have a single spatio-temporal predictor in the form of an image (spatstat im class). I generate such toy data here:
set.seed(1234)
library(spatstat)
# Number of time periods. For each time period I observe a point pattern.
time_points <- 200
# spatial window is the [0, 1] square
win <- owin(xrange = c(0, 1), yrange = c(0, 1))
# A spatial variable
surf <- as.im(X = function(x, y) {exp(- (x ^ 2 + y ^ 2) * 2)}, W = win)
plot(surf)
# Scaling the spatial variable by a smooth function of time to make it
# vary across time periods.
fn_time <- log(sin(1 : time_points / 20) + 2)
log_lambda <- lapply(1 : time_points, function(t) fn_time[t] * surf)
# Generating the point patterns:
pps <- NULL
for (tt in 1 : time_points) {
pps[[tt]] <- spatstat::rpoispp(lambda = exp(log_lambda[[tt]]))
}
# Taking a look at the number of points in each pattern:
hist(sapply(pps, function(x) x$n))
Then, we can use mppm to fit a poisson process model for the point pattern data on the predictor. We can do so with the following code:
mod_dta <- hyperframe(Points = pps, Pred = log_lambda)
pp_mod <- mppm(Points ~ Pred, data = mod_dta)
In order to extract the estimated conditional intensity function we can use the predict.mppm
function:
mod_dta$fitted_cif <- predict.mppm(pp_mod, type = 'cif', ngrid = 100)$cif
which we can use to extract the expected number of points for each time period as the integral of the conditional intensity function:
plot(with(mod_dta, integral(fitted_cif)))
However, this process can be slow if the number of time periods time_points
is large, and the value for ngrid
in predict.mppm
is set to large (to improve accuracy).
I wonder if the fitted expected number of points can be extracted directly from the mppm fit (pp_mod
in the example) without having to define the conditional intensity as an image and integrate over it.
Thank you for your help.