profilepl() for complex hybrid Gibbs models

I am having some trouble figuring out how to use profilepl() in the R spatstat package for more complicated hybrid Gibbs models. For example, I am trying to build a hybrid model with a hardcore and multiple Geyer interactions.

I have tried doing something like this:

RR <- c(expand.grid(r=seq(4,8, by=1), sat = 1:2), expand.grid(r1=seq(8,12,  by=1), sat1 = 1:2))

MS <- function(r, sat, r1, sat1) { Hybrid(A=Hardcore(NA), B=Geyer(r=r, sat = sat),  C=Geyer(r=r1, sat = sat1)) }

fit <- profilepl(RR, MS, swedishpines ~ polynom(x,y,2), correction = "isotropic", aic = FALSE)

I am not sure if I am doing this right, and it seems that there may be a better way to write the interaction function (RR and MS) that what I have done. For the Geyer interactions, I currently have to give it two "ranges" of interactions to estimate interaction radii from, and I would like to know if I could just provide one "range" of numbers (e.g. 4 - 12) and have profilepl select the two best interaction radii within this interval.

Plotting the profile likelihood also gives a weird output, not sure if there is a way to better visualize this, especially if it can indicate both optimal Geyer interaction parameters: https://i.sstatic.net/mLC0B18D.png

Any help is appreciated! Thank you

Solution

Your code is mostly correct and is producing a valid fitted object.

However, the code is not searching over all possible combinations of the parameter values, but only a few selected combinations. This is causing the strange-looking plot.

The first argument of profile should be a data frame, in which each row provides one set of parameter values to be considered. The rows of the data frame should enumerate the parameter space (all parameter values) over which you want to optimize.

In your code, the object RR contains 10 sets of values for the parameters r, sat, r1, sat1. The parameter values can also be extracted from the fitted object:

fit$param
   r sat  r1  sat1
1  4   1  8    1
2  5   1  9    1
3  6   1 10    1
4  7   1 11    1
5  8   1 12    1
6  4   2  8    2
7  5   2  9    2
8  6   2 10    2
9  7   2 11    2
10 8   2 12    2

In this data frame, the values of sat and sat1 are always equal to one another, and r1 is always equal to r + 4. If that's what you intended, then there are really only two "free variables", namely r and sat. This explains why the plot produced by plot(fit) looks weird.

It would be better to code the interaction as a function of the two free variables r and sat:

RR <- expand.grid(r=seq(4,8, by=1), sat = 1:2)
MS <- function(r, sat) { 
   Hybrid(A=Hardcore(NA), 
          B=Geyer(r=r, sat = sat),  
          C=Geyer(r=r+4, sat = sat)) }

Then refit using the same call to profilepl. This should produce a sensible looking plot.