how to visualize the effect of a variable in a multiple point process model fitted using gam using spatstat?

I am fitting multiple point process model fitted using gam in spatstat (version 3.0-7). To be clear I am using a call of the general form:

mppm(formula, data=hyperframe, eps=0.5, rbord=0.5, use.gam=TRUE)

Thus my model is Poisson process.

What I would like to do when the model is fitted is to plot:

the effect of a given predictor variable in an mgcv::plot.gam fashion
maps of the spatial trend.

However I get these error messages When trying to use functions

effectfun() : Error: First argument 'model' should be a fitted model of class ‘ppm’, ‘kppm’, ‘lppm’, ‘dppm’, ‘rppm’ or ‘profilepl’
plot.mppm() : Error: This calculation is not supported for GAM fits

I am considering the following workaround options:

Subfit approach in spatstat

re-fit the model to each pattern in my hyperframe individually (because subfits does not work for gam): I apply the same model formula used in the mppm model and run it in a ppm for each pattern. But in doing so, I might modify the k parameters independently for each pattern because the k-value used in mppm might give an error with some variables.
use effectfun and plot on each individually fitted ppm models. Here I am afraid that I could possibly end up with splines shapes different for each pattern, thus complication the interpretation

hack mppm and use mgcv

run the mppm with use.gam=T model I want.
extract the data.frame from the mppm object using my_model$Fit$moadf
fit the model with mgcv::gam using the call that is used internally by mppm:

gam(fmla, family = quasi(link = log, variance = mu), weights = .mpl.W * caseweight, data = moadf, subset = (.mpl.SUBSET == "TRUE"), control = ctrl)

at that point, I could use the mgcv functions to visualize the plots I need (at least the effect of a given predictor variable.

My questions

Is any of my described approaches good? One of them to be preferred? Or, are there other (better) options to achieve my goal?

Solution

Short answer: extract the fitted gam object directly as

m <- my_model$Fit$FIT

Then do plot(m, ...) which invokes mgcv::plot.gam.

Long answer:

Currently, gam fits (produced when use.gam=TRUE) are not fully supported in mppm.

Currently we use the subfits mechanism to perform many tasks for mppm objects, including plotting and predicting the model. This is our way to avoid (for the moment) having to re-write all of the code that supports the ppm class, all over again to support the mppm class.

The function subfits is a hack which involves making "fake" objects of class ppm (they are fake because they were not actually obtained by fitting a model to a single point pattern). This fakery has limitations and it does not work at all when the model is fitted by gam.

Ultimately the code supporting the mppm class will be rewritten as native code and this problem will be resolved. Until then, I suggest you use the approach described in the short answer.