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: https://i.sstatic.net/mLC0B18D.png

Any help is appreciated! Thank you