CODEWITHSUNDEEP
rspatialspatstat
proximacentauri
proximacentauri

Reputation: 1879

Simulate a list of kppm objects in R spatstat

I would like to use envelope simulation on a list of kppm objects.

I have a point dataset that is in a grid pattern. I would like to shift each point such that its position in the cell is random, then run a spatstat simulation to determine if the points are randomly distributed (obviously will be at sub cell size scale). Then randonly shift shift each point inside the cell, and re-run the simulation - repeating n times and generate an envelope of all the simulations. There are 2 ways to do this 1) using a list of point sets and 2) using simulate.

Method 1) - on a list called my_kppm_list

Lmat = envelope(fitM, Lest, nsim = 19, global=FALSE, simulate=my_kppm_list)
plot(Lmat)

How to create a list of cluster process models (kppm)?

Naive way:

my_list <- list(kppm_0, kppm_1)

This fails when trying to run simulation:

Lmat = envelope(fitM, Lest, nsim = 19, global=FALSE, simulate=my_list)

Error in envelopeEngine(X = X, fun = fun, simul = simrecipe, nsim = nsim,  : 
  ‘simulate’ should be an expression, or a list of point patterns of the same kind as X

I can convert a list to .ppm

fitM_0 <- kppm(create_pts(), ~1, "MatClust")
fitM_1 <- kppm(create_pts(), ~1, "MatClust")
my_list <- list(fitM_0, fitM_1)
ppm_list <- lapply(my_list, as.ppm)

But trying to convert to kppm fails with an error

kppm_list <- lapply(my_list, as.kppm)

Method 2) apply a function in simulation such that a random shift is applied to each point then simulation run, and envelope of all simulations is used (see page 399 of Baddelley et al. Spatial Point Patterns book (2016)):

e_rand <- function(){
  j_x <- matrix(unlist(runif(dim(c_df)[1], min=-10, max=10)), ncol = 1, byrow = T)
  j_y <- matrix(unlist(runif(dim(c_df)[1], min=-10, max=10)), ncol = 1, byrow = T)
  x_j<- c_df[,1]+j_x
  y_j<- c_df[,2]+j_y
  
  c_j <- ppp(x = x_j, y = y_j, window = window)
  return(c_j)
}

Lmat = envelope(fitM, Lest, nsim = 19, global=TRUE, simulate=e_rand)

However I found that the null model (red dashed line in output plot) had kinks in it when simulate is added - kinks that do not appear without simulate.

Upvotes: 0

Views: 113

Answers (1)

Adrian Baddeley
Adrian Baddeley

Reputation: 2973

How to create a list of cluster process models (kppm)?

This is not the problem. In your example you are successfully creating a list of objects of class kppm. If fit1 and fit2 are fitted models of class kppm, then

m <- list(fit1, fit2)

is a list of objects of class kppm.

The problem is that you then pass this list as an argument to the function envelope which does not accept this format. The error message from envelope.ppp or envelope.kppm says that the argument simulate must be either a fitted model, or a list of point patterns.

An envelope is constructed by generating simulated point patterns, computing a summary function for each simulated pattern, and computing the upper and lower extremes of these summary functions. The argument simulate (if it is given) is a recipe for generating the simulated point patterns. If simulate is an expression like expression(runifpoint(42)) then this expression will be evaluated nsim times to produce nsim different point patterns. If simulate is a fitted model, then nsim simulated realisations of the model will be generated. If simulate is a list of point patterns, then they will be taken as the simulated random patterns.

It is unclear what you want to do with your list of models. Do you want to construct a single envelope, or a list of envelopes?

  • Option 1: you have a list m of models of class kppm. For each of these models m[[i]], you want to construct an envelope e[[i]], where the limits are determinedby simulation from m[[i]]`.
  • Option 2: you have a list m of models of class kppm. For each model m[[i]] you want to generate one point pattern, say X[[i]], and build an envelope e using these patterns.

For option 1, type something like

e <- anylapply(m, function(fit) { 
        envelope(Y, Lest, nsim = 19, global=FALSE, simulate=fit)})

For option 2,

X <- solapply(m, simulate, nsim=1, drop=TRUE)
e <- envelope(Y, Lest, nsim=19, global=FALSE, simulate=X)

If you wanted something else, please clarify.

Upvotes: 1

Related Questions