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:


# 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)

# 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.


Unfortunately I don't think you can avoid this in the general case. If your predictor in general is different from point pattern to point pattern without any relation, you will simply have to calculate the intensity in each pixel and sum up. This is what happens with the calls to predict and integral. In your specific example there is a relationship between the predictor for one pattern and the predictor for another pattern, and in some cases you may be able to exploit such a relationship to calculate the predicted intensity for one case based on the other. It depends entirely on the relation, and I won't look at any details in the given case since it is just a toy example.

  • $\begingroup$ Thank you for your response. I have found practically that getting an approximate (but good enough) value for the integral of the intensity does not require as large of an ngrid value in predict. So I can move forward with that. $\endgroup$ – GPa Dec 7 '20 at 14:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.