Lab 4: Determinantal and Gibbs models

library(spatstat)

Exercise 1

In this question we fit a Strauss point process model to the swedishpines data.

We need a guess at the interaction distance $R$. Compute and plot the $K$-function of the dataset and choose the value $r$ which maximises the discrepancy $|K(r)-\pi r^2|$.

We plot the above function which we want to maximize.
```
plot(Kest(swedishpines), abs(iso - theo) ~ r, main = "")
```
As seen from the plot, the maximum lies around $r = 9$ by eye. We find the optimum explicitly like follows:
```
discrep <- function(r) {
  return(abs(as.function(Kest(swedishpines))(r) - pi*r^2))
}
res <- optimise(discrep, interval = c(0.1, 20), maximum = TRUE)
print(res)
```
```
## $maximum
## [1] 9.84372
## 
## $objective
## [1] 150.6897
```
```
R <- res$maximum
```
This corresponds nicely with the plot.
Fit the stationary Strauss model with the chosen interaction distance using
```
ppm(swedishpines ~ 1, Strauss(R))
```
where R is your chosen value.

As we have assigned R, we simply write:
```
fit <- ppm(swedishpines ~ 1, Strauss(R))
```
Interpret the printout: how strong is the interaction?
```
print(fit)
```
```
## Stationary Strauss process
## 
## First order term:  beta = 0.08310951
## 
## Interaction distance:    9.84372
## Fitted interaction parameter gamma:  0.2407279
## 
## Relevant coefficients:
## Interaction 
##   -1.424088 
## 
## For standard errors, type coef(summary(x))
```
As seen, the $\gamma = `round(exp(coef(fit)[2]), 2)$ parameter is quite small. Thus there seems to be a strong negative association between points within distance R of each other. A $\gamma$ of $0$ implies the hard core process whereas $\gamma = 1$ implies the Poisson process and thus CSR.
Plot the fitted pairwise interaction function using plot(fitin(fit)).

The pairwise interaction function become:
```
plot(fitin(fit))
```

Exercise 2

In Question 1 we guesstimated the Strauss interaction distance parameter. Alternatively this parameter could be estimated by profile pseudolikelihood.

Look again at the plot of the $K$-function of swedishpines and determine a plausible range of possible values for the interaction distance.
```
plot(Kest(swedishpines), main = "")
```
A conservative range of plausible interaction distances seems to be 5 to 12 meters.
Generate a sequence of values equally spaced across this range, for example, if your range of plausible values was $[0.05, 0.3]$, then type
```
rvals <- seq(0.05, 0.3, by=0.01)
```
We generate the numbers between 5 and 12.
```
rvals <- seq(5, 12, by = 0.1)
```
Construct a data frame, with one column named r (matching the argument name of Strauss), containing these values. For example
```
D <- data.frame(r = rvals)
```
OK,
```
D <- data.frame(r = rvals)
```

Execute

fitp <- profilepl(D, Strauss, swedishpines ~ 1)

to find the maximum profile pseudolikelihood fit.

OK, let’s execute it:

fitp <- profilepl(D, Strauss, swedishpines ~ 1)

## (computing rbord)

## comparing 71 models...

## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,  71.

## fitting optimal model...

## done.

Print and plot the object fitp.

print(fitp)

## profile log pseudolikelihood
## for model:  ppm(swedishpines ~ 1,  interaction = Strauss)
## fitted with rbord = 12
## interaction: Strauss process
## irregular parameter: r in [5, 12]
## optimum value of irregular parameter:  r = 9.8

plot(fitp)

Compare the computed estimate of interaction distance $r$ with your guesstimate. Compare the corresponding estimates of the Strauss interaction parameter $\gamma$.
```
(Ropt <- reach(as.ppm(fitp)))
```
```
## [1] 9.8
```
The $r = 9.8$ is consistent with the previous guesstimate.
Extract the fitted Gibbs point process model from the object fitp as
```
bestfit <- as.ppm(fitp)
```
OK, let’s do that:
```
bestfit <- as.ppm(fitp)
```

Exercise 3

Modify Question 1 by using the Huang-Ogata approximate maximum likelihood algorithm (method="ho") instead of maximum pseudolikelihood (the default, method="mpl").

fit.mpl <- ppm(swedishpines ~ 1, Strauss(R), method = "mpl")
fit.ho  <- ppm(swedishpines ~ 1, Strauss(R), method = "ho")

## Simulating... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,  100.
## Done.

print(fit.ho)

## Stationary Strauss process
## 
## First order term:  beta = 0.1126732
## 
## Interaction distance:    9.84372
## Fitted interaction parameter gamma:  0.2257317
## 
## Relevant coefficients:
## Interaction 
##   -1.488408 
## 
## For standard errors, type coef(summary(x))

print(fit.mpl)

## Stationary Strauss process
## 
## First order term:  beta = 0.08310951
## 
## Interaction distance:    9.84372
## Fitted interaction parameter gamma:  0.2407279
## 
## Relevant coefficients:
## Interaction 
##   -1.424088 
## 
## For standard errors, type coef(summary(x))

The fits are very similar.

Exercise 4

Repeat Question 2 for the inhomogeneous Strauss process with log-quadratic trend. The corresponding call to profilepl is

fitp <- profilepl(D, Strauss, swedishpines ~ polynom(x,y,2))

## (computing rbord)

## comparing 71 models...

## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,  71.

## fitting optimal model...

## done.

fitp2 <- profilepl(D, Strauss, swedishpines ~ polynom(x,y,2))

## (computing rbord)

## comparing 71 models...

## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,  71.

## fitting optimal model...

## done.

print(fitp2)

## profile log pseudolikelihood
## for model:  ppm(swedishpines ~ polynom(x,  y,  2),  interaction = Strauss)
## fitted with rbord = 12
## interaction: Strauss process
## irregular parameter: r in [5, 12]
## optimum value of irregular parameter:  r = 9.8

(bestfit <- as.ppm(fitp2))

## Nonstationary Strauss process
## 
## Log trend:  ~x + y + I(x^2) + I(x * y) + I(y^2)
## 
## Fitted trend coefficients:
##   (Intercept)             x             y        I(x^2)      I(x * y) 
## -5.156605e+00  4.269469e-02  8.582459e-02 -5.324026e-06 -7.955428e-04 
##        I(y^2) 
## -5.506409e-04 
## 
## Interaction distance:    9.8
## Fitted interaction parameter gamma:  0.2390264
## 
## Relevant coefficients:
## Interaction 
##   -1.431181 
## 
## For standard errors, type coef(summary(x))

reach(bestfit)

## [1] 9.8

Exercise 5

Fit a stationary Gaussian DPP model to the swedishpines data using the dppm() command.

fit <- dppm(swedishpines, dppGauss())

Generate 39 simulations of the fitted model (simulations may take some time to generate).

sim <- simulate(fit, nsim = 39)

Make pointwise envelopes for a summary statistic of your choice (e.g. Kest()) using the generated 39 simulations (see the simulate argument in envelope.ppp())

env <- envelope(swedishpines, Lest, simulate = sim, nsim = 39)

## Extracting 39 point patterns from list  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,  39.
## 
## Done.

Plot the envelopes and consider whether you find the fitted model to be reasonable.

plot(env)

Lab 4: Determinantal and Gibbs models – solution

Exercise 1

Exercise 2

Exercise 3

Exercise 4

Exercise 5