Determinantal and Gibbs models

Models for inhibited data

library(spatstat)
plot(cells)

jp <- residualspaper$Fig1
plot(jp)

Determinantal point process models

Determinantal point process models are a relatively new class of models with many attractive theoretical properties.

These models can be fitted in spatstat using the function dppm.

dppm(jp ~ x+y, dppGauss)
## Inhomogeneous determinantal point process model
## Fitted to point pattern dataset 'jp'
## Fitted by minimum contrast
##  Summary statistic: inhomogeneous K-function
## 
## Log intensity:  ~x + y
## 
## Fitted trend coefficients:
## (Intercept)           x           y 
## 0.591839808 0.014329205 0.009643885 
## 
## Fitted DPP model:
## Gaussian determinantal point process model
## The parameters are: lambda = an image, alpha = 0.06402, d = 2
## The parameter lambda specifies the intensity of the process.
## The parameter d specifies the dimension of the state space.
fit <- dppm(jp ~ polynom(x,y,2),
            dppMatern,
        statistic="pcf", statargs=list(stoyan=0.2))
## Warning: The value of the empirical function 'pcf' for r= 0 was Inf. Range of r
## values was reset to [0.0048828125, 2.5]
plot(predict(fit))      

plot(simulate(fit))

Gibbs models

Gibbs models were developed in theoretical physics to describe the behaviour of molecular gases. A point pattern \(x\) represents a spatial configuration of molecules. The probability of a particular configuration \(x\) is \[ p(x) = Z \exp(- U(x)) \] where \(U(x)\) is the potential energy of the configuration, and \(Z\) is a normalising constant. In fact \(p(x)\) is a probability density relative to the completely random (homogeneous Poisson) point process.

To visualise this, imagine that we first generate an infinite “ensemble” of realisations of the homogeneous Poisson process. Then each realisation is either deleted or retained (in its entirety) depending on its potential energy; a realisation \(x\) is retained with probability \(\exp(-U(x))\). Then what remains is an ensemble of realisations of the Gibbs process.

The simplest example is the hard core process in which the points represent the centres of discs of diameter \(d\) which cannot overlap. A realisation \(x\) has potential energy \(U(x) = -\infty\) if any pair of points in \(x\) lies closer than distance \(d\); otherwise it has potential \(U(x) = 0\). Now generate an infinite ensemble of realisations of the Poisson process. Then delete any configuration which contains a pair of points closer than distance \(d\). The remaining realisations are an ensemble of realisations of the hard core process.

Gibbs models can be fitted to point pattern data by maximising Besag’s pseudolikelihood. This is performed by ppm.

ppm(cells~ 1, Hardcore())
## Stationary Hard core process
## 
## First order term:  beta = 282.7782
## 
## Hard core distance:  0.08168525
## 
## For standard errors, type coef(summary(x))
minnndist(cells)
## [1] 0.08363014
ppm(cells ~ 1, Strauss(0.1))
## Stationary Strauss process
## 
## First order term:  beta = 1138.136
## 
## Interaction distance:    0.1
## Fitted interaction parameter gamma:  0.0050219
## 
## Relevant coefficients:
## Interaction 
##    -5.29395 
## 
## For standard errors, type coef(summary(x))
fit <- ppm(cells ~ 1, Strauss(0.1))
plot(simulate(fit, nsim=4))
## Generating 4 simulated patterns ...1, 2, 3,  4.

plot(pcfinhom(jp))

minnndist(jp)
## [1] 0.1410937
ppm(jp ~ x+y, Strauss(0.2))
## Nonstationary Strauss process
## 
## Log trend:  ~x + y
## 
## Fitted trend coefficients:
## (Intercept)           x           y 
## 0.679292879 0.006806174 0.024664271 
## 
## Interaction distance:    0.2
## Fitted interaction parameter gamma:  0.6367406
## 
## Relevant coefficients:
## Interaction 
##   -0.451393 
## 
## For standard errors, type coef(summary(x))
ppm(jp ~ x+y, Strauss(0.5))
## Nonstationary Strauss process
## 
## Log trend:  ~x + y
## 
## Fitted trend coefficients:
##  (Intercept)            x            y 
##  0.364642729 -0.004808391  0.033863194 
## 
## Interaction distance:    0.5
## Fitted interaction parameter gamma:  1.1799626
## 
## Relevant coefficients:
## Interaction 
##   0.1654827 
## 
## For standard errors, type coef(summary(x))
## 
## *** Model is not valid ***
## *** Interaction parameters are outside valid range ***