# Maximum Pseudo-Likelihood estimator for soft-core point process

I am trying to fit a soft-core point process model on a set of point pattern using maximum pseudo-likelihood. I followed the instructions given in this paper by Baddeley and Turner

And here is the R-code I came up with

library(deldir)
library(tidyverse)
library(fields)

#MPLE

# irregular parameter k

k <- 0.4

## Generate dummy points 50X50. "RA" and "DE" are x and y coordinates

dum.x <- seq(ramin, ramax, length = 50)
dum.y <- seq(demin, demax, length = 50)
dum <- expand.grid(dum.x, dum.y)
colnames(dum) <- c("RA", "DE")

## Combine with data and specify which is data point and which is dummy, X is the point pattern to be fitted

bind.x <- bind_rows(X, dum) %>%
mutate(Ind = c(rep(1, nrow(X)), rep(0, nrow(dum))))

## Calculate Quadrature weights using Voronoi cell area

w <- deldir(bind.x$$RA, bind.x$$DE)$$summary$$dir.area

## Response
y <- bind.x$Ind/w # the sum of distances between all pairs of points (the sufficient statistics) tmp <- cbind(bind.x$$RA, bind.x$$DE) t1 <- rdist(tmp)^(-2/k) t1[t1 == Inf] <- 0 t1 <- rowSums(t1) t <- -t1 # fit the model using quasipoisson regression fit <- glm(y ~ t, family = quasipoisson, weights = w)   However, the fitted parameter for t is negative which is obviously not a correct value for a softcore point process. Also, my point pattern is actually simulated from a softcore process so it does not make sense that the fitted parameter is negative. I tried my best to find any bugs in the code but I can't seem to find it. The only potential issue I see is that my sufficient statistics is extremely large (on the order of 10^14) which I fear may cause numerical issues. But the statistics are large because my observation window spans a very small unit and the average distance between a pair of points is around 0.006. So sufficient statistics based on this will certainly be very large and my intuition tells me that it should not cause a numerical problem and make the fitted parameter to be negative. Can anybody help and check if my code is correct? Thanks very much! • Did you try to reproduce your results directly with the spatstat package? – Ege Rubak Feb 24 '19 at 20:09 • @EgeRubak I tried using ppm with softcore interaction. The parameter for simulation is$\log(\beta) = \log(6000) = 8.699$and$\sigma = 0.0065$. The fitted parameter from the simulated data is$\log(\beta) = \log(4307) = 8.37$and$\sigma = 0.0067\$. So my data is not wrong. But I really can't seem to find any mistakes in my code. – davidolohowski Feb 24 '19 at 20:42
• Not by my computer so can't offer detailed help. You both divide by weights in y and use weights in glm. Is that correct? In general my best advise is to reproduce the spatstat results as closely as possible. Make a quadrature scheme in spatstat and extract the coordinates and weights and use them in your call to glm. You can find a bunch of glm info deep inside the list returned by ppm maybe that could help with debugging also. Let me know how it goes... Also think about edge correction when trying to reproduce results... – Ege Rubak Feb 25 '19 at 4:34
• @EgeRubak Dividing y by weights and using weights in glm is just me following the procedure from the paper as I understood it and I tried removing division or removing weights, it still gives me negative values. For the edge correction, I think it might not affect the results that much as to make the fitted parameter negative since I used border correction of 0.0001 (which is very small compared to the whole observation window) in ppm and the results from the last comment was what I got. I will try to look into the results returned by ppm and see how it goes. Thanks! – davidolohowski Feb 25 '19 at 4:56

The procedure described in the paper is implemented in the spatstat package, together with additional tools to ensure numerical stability, and to provide diagnostic information.
In spatstat the soft core sufficient statistic -\sum_j d_{ij}^{-2/kappa} is rescaled to -\sum_j (sigma_0/d_{ij})^{2/kappa} where sigma_0 is a trial value estimated from the nearest-neighbour distances. This generally keeps the numerical oscillations manageable. See the help for Softcore in the spatstat package.