Parallelization & processing of ppm() in spatstat?

Question

I am attempting to fit a point pattern with ppm() in spatstat, version 1.60-1, using a combination of first-order covariates (as im objects) and a second-order interaction term, in this case an Area-Interaction model.

For my data, comprising: ~300 points, 3 covariate im.objects of ~110 MB each, and the area interaction term at r = 15000 map units (meters), processing time is proving to be prohibitively long - over 72 hours at this point.

Using some of the data included in spatstat itself:

require(spatstat)
res <- vector()

for (i in seq(1, 9, 0.5)) {
    start.time <- Sys.time()
    ppm(swedishpines, ~1, AreaInter(r=i))
    end.time <- Sys.time()
    time.taken <- end.time - start.time
    res[i] <- time.taken
}

plot(res)

Output:

it can be seen that the processing time increases roughly linearly. I imagine that the size of the point pattern, window, inclusion of covariates, and apparently most importantly, the value of r, heavily affect this too.

My reason for believing this is that I previously ran ppm() with a larger point pattern, similar covariates, but a much smaller value of r (~7000 instead of ~15000)

It appears to be relatively simple to parallelize the envelope() function and pool the results with the corresponding function, but I am looking for advice on how to parallelize ppm() and be able to combine the output fitted models? There doesn't seem to be a direct implementation of this task for ppm objects.

Alternatively, if someone is more familiar with the implementation of the Area Interaction model in spatstat, I'd be grateful to know if I am just chasing cars with such a large value of r, and if I would be better off trying a different approach.

This is just on my home machine, running on an old i5-3570K quad core with 8 GB of RAM. — p.riris

Adrian Baddeley Adrian Baddeley · Accepted Answer · 2019-07-09T05:25:06

In the above code, you are selecting the value of the disc radius parameter r by making repeated calls to ppm in a loop. This is wasteful, because information is re-computed each time. It would be much faster to use the spatstat function profilepl which is designed for this purpose. It uses efficient shortcuts, including special tricks for the area interaction model.

The area-interaction model is not a second order interaction; it has infinite order. It is therefore quite complicated to compute.

Another strategy for selecting r would be to temporarily replace the area interaction by the Penttinen interaction, a second order interaction, which is close to the area interaction. Being a second order interaction the Penttinen model can be fitted very quickly so you can rapidly determine a good value for r.

The development of faster algorithms for AreaInter is on our list of things to do. But no-one is funding us for this..

Parallelization & processing of ppm() in spatstat?

Output:

1 Answers