model selection for Mixed Poisson model for replicated point pattern analysis

Question

I have fitted a mixed Poisson Point Process Model to Several Point Patterns using the mppm function from the R package spatstat.

I am able to fit the model to a full model with several predictors which are im objects in which pixels contain either distance to or density of the predictor.

It seems possible to do some model selection if the model does not include a random part. On the contrary, the introduction on the Random component seems to complicate things and that available methods stop working.

How model selection should be performed when analysing Poisson processes on replicated point pattern? In particular how should it be performed on mixed models? Is it possible

While this post arise from technical issues in using spatstat functions, it doesn't want to be only a programming question, rather I'm interested in finding solution for do model selection in this context.

Following, a reproducible example. Data can be found here:

library(spatstat)
load('data4stack-model_sel-20200516.RData')

urc

# without Random factor
# ------------------------

m50a<-mppm(urc_ppp~ pref_dist+kelp50+bushy50+ota50+star_dens+dsh_dens + are_dens+ dol_dens, data=urc) 
m150a<-mppm(urc_ppp~ pref_dist+kelp150+bushy150+ota150+dsh_dens + are_dens+ dol_dens, data=urc)
# Warning messages:
# 2: Values of the covariates ‘dsh_dens’, ‘are_dens’, ‘dol_dens’ were NA or undefined at 0.64% (34 out of 5344) of the quadrature points. Occurred while executing: 
# 3: data contain duplicated points

Regarding the warning on duplicated points: I personally digitized the points, double checked them, and verified their coordinates. There are no duplicated points, however for some reason (rounding?) it give this warning.

Regarding the warning on NA in covariates: I believe it is due to the fact that some covariates are not defined in some areas where point of the response variable may be. I could do the analysis without those points, so if the model fitting does not include those points I have no problems. I guess that other NAs, my originate in some discrepancy introduced when discretizing into pixels to produce the im variables.

In anycase these 2 warnings do not seem to affect model selection for models with no random component. See the rest:

anova(m50a, test='LRT')
m50b<-update(m50a, .~.-ota50)

anova(m50a,m50b, test='LRT') # This works

step(m50a) # This seems to work

AIC(m50a)
extractAIC(m150a) 
AIC(m50a,m50b,m150a) # Just give AIC for the first model
# Work around for AIC-based model comparison ?
AIC(m50a); AIC(m50b); AIC(m150a)


# with Random factor
# --------------------

m50r_a<-mppm(urc_ppp~ pref_dist+kelp50+bushy50+ota50+dsh_dens+are_dens+dol_dens, random=~1|ID,data=urc) 
AIC(m50r_a) ; extractAIC(m50r_a)  # both give NA

drop1(m50r_a) # it runs but it doesn't do any selection since AIC are not provided
step(m50r_a) # it doesn't even run

# hypothesys testing:
m50r_b<-update(m50r_a, .~.-are_dens)
anova(m50r_a,m50r_b, test='Chisq') 
# Error in anova.lme(structure(list(modelStruct = structure(list(reStruct = structure(list( : 
#   objects must inherit from classes "gls", "gnls", "lm", "lmList", "lme", "nlme", "nlsList", or "nls"
# Error in anova.mppm(m50r_a, m50r_b, test = "Chisq") : anova failed

Any suggestion? Thanks

Answer 1

This question is mostly about the spatstat package.

In the function mppm, models with random effects are fitted by penalised quasi-likelihood using the function glmmPQL from the MASS package. The resulting fitted objects returned from glmmPQL have likelihood equal to NA. Attempting to apply anova to the object returned from glmmPQL yields an error message that anova is not supported for these objects. (That's a message from the MASS package.)

The issue is that a Penalised Quasilikelihood is not a likelihood, so it is unclear how to do inference (including model selection) for models fitted by PQL. This is a research question.

In previous versions of lme, objects fitted by MASS::glmmPQL were accepted by anova.lme and other methods, and the commands you typed in your question would have executed without failure. However, they are no longer accepted, which is leading to this error message. I presume that the authors of lme have done this to prevent exactly this kind of usage - because the calculations in anova.lme are strictly speaking invalid for a glmmPQL fit.

Out of respect for the wisdom of the authors of lme, I will not try to work around their code. Instead I will tweak the code in spatstat so that it gives a more useful error message.

Incidentally - your data do include duplicated points:

with(urc, anyDuplicated(urc_ppp))
#  1      5       22
#  FALSE  FALSE   TRUE

According to with(urc, summary(urc_ppp)) the coordinates were recorded to the nearest 0.1 units (in windows about 20 units across) so this is likely to be the reason for the existence of duplicated points.

However, the software is supposed to cope with this.