如何在spatstat中用effectfun计算GAM模型的标准误？

问题背景 #

我用spatstat 3.0.2分析史前五个时期考古遗址的分布模式，构建了包含海拔、坡度等协变量的点过程模型，其中涵盖GAM模型，需要对比GAM与非GAM模型的结果。但发现effectfun()函数对GAM模型无法计算标准误，仅能输出拟合线，已用bei数据集复现该问题：

data(bei)
data(bei.extra)
elev <-  bei.extra$elev
fitNoGam <- ppm(bei~elev)
fitGam <- ppm(bei~s(elev),use.gam=TRUE)
par(mfrow=c(1,2))
plot(predict(fitNoGam))
plot(predict(fitGam))
plot(effectfun(fitNoGam,"elev",se.fit=T))
plot(effectfun(fitGam,"elev")) # 无法计算标准误

解决方案 #

spatstat中用use.gam=TRUE拟合的ppm模型本质是mgcv包的GAM对象，effectfun()目前不支持为这类模型计算标准误，我们可以借助mgcv的predict.gam()函数手动获取预测值和标准误，再绘制带置信区间的效果曲线：

# 生成elev的取值序列（覆盖数据的全范围）
elev_seq <- seq(min(elev), max(elev), length.out = 100)
# 构建预测用的数据框
newdata <- data.frame(elev = elev_seq)
# 从GAM模型中获取预测值与标准误（type="link"对应线性预测器，和effectfun输出逻辑一致）
pred_result <- predict(fitGam, newdata = newdata, se.fit = TRUE, type = "link")
# 计算95%置信区间
lower_ci <- pred_result$fit - 1.96 * pred_result$se.fit
upper_ci <- pred_result$fit + 1.96 * pred_result$se.fit

# 绘制拟合线与置信区间
plot(elev_seq, pred_result$fit, type = "l", 
     xlab = "海拔", ylab = "线性预测器", 
     main = "GAM模型中海拔的效应（带95%置信区间）")
lines(elev_seq, lower_ci, lty = 2, col = "gray50")
lines(elev_seq, upper_ci, lty = 2, col = "gray50")

补充说明 #

• 上述代码输出和effectfun()的拟合线逻辑一致，均基于线性预测器的结果；
• 如果需要在原始强度尺度上展示结果，只需将type="link"改为type="response"，同时对置信区间做对应逆变换（比如用exp()，因为ppm默认使用对数链接）。

内容的提问来源于stack exchange，提问作者A.M. Xenaki