0

My question is closely related to this previous one: Simulation-based hypothesis testing on spatial point pattern hyperframes using "envelope" function in spatstat

I have obtained an mppm object by fitting a model on several independent datasets using the mppmfunction from the R package spatstat. How can I study its envelope to compare it to my observations ?

I fitted my model as such:

data <- listof(NMJ1,NMJ2,NMJ3)
data <- hyperframe(X=1:3, Points=data)
model  <- mppm(Points ~marks*sqrt(x^2+y^2), data)

where NMJ1, NMJ2, and NMJ3 are marked ppp and are independent realizations of the same experiment.

However, the envelope function does not accept inputs of type mppm:

> envelope(model, Kcross.inhom, nsim=10)
Error in UseMethod("envelope") : 
  no applicable method for 'envelope' applied to an object of class "c('mppm', 'list')"

The answer provided to the previously mentioned question indicates how to plot global envelopes for each pattern, and to use the product rule for multiple testing. However, my fitted model implies that my 3 ppp objects are statistically equivalent, and are independent realizations of the same experiment (ie no different covariates between them). I would thus like to obtain one single plot comparing my fitted model to my 3 datasets. The following code:

gamma= 1 - 0.95^(1/3)
nsims=round(1/gamma-1)
sims <- simulate(model, nsim=2*nsims)
SIMS <- list()
for(i in 1:nrow(sims)) SIMS[[i]] <- as.solist(sims[i,,drop=TRUE])
Hplus <- cbind(data, hyperframe(Sims=SIMS))

EE1 <- with(Hplus, envelope(Points, Kcross.inhom, nsim=nsims, simulate=Sims))

pool(EE1[1],EE1[2],EE1[3])

leads to the following error:

Error in pool.envelope(`1` = list(r = c(0, 0.78125, 1.5625, 2.34375, 3.125,  : 
  Arguments 2 and 3 do not belong to the class “envelope”
Camille Gontier
  • 145
  • 5

1 Answers1

1

Wrong type of subset index. Use

pool(EE1[[1]], EE1[[2]], EE1[[3]])

or just

pool(EE1)

These would have given an error message that the envelope commands should have been called with savefuns=TRUE. So you just need to change that step as well.

However, statistically this procedure makes little sense. You have already fitted a model, which allows for rigorous statistical inference using anova.mppm and other tools. Instead of this, you are generating simulated data from the fitted model and performing a Monte Carlo test, with all the fraught issues of multiple testing and low power. There are additional problems with this approach - for example, even if the model is "the same" for each row of the hyperframe, the patterns are not statistically equivalent unless the windows of the point patterns are identical, and so on.

Adrian Baddeley
  • 1,956
  • 1
  • 8
  • 7
  • Thanks a lot. Indeed, these plots are simply meant as a sanity check (my model is validated by comparing its AIC with these of null-hypothesis models) and to have a visual representation of the success of the fitting procedure (for a scientific publication). – Camille Gontier Sep 23 '21 at 11:25