R中spatstat的kppm拟合LGCP时置信区间全为NA的问题求助

云上老周

2026-7-3

解决kppm拟合LGCP后confint返回全NA的问题

核心问题分析

使用kppm()以最小对比度法（method="mincon"）拟合对数高斯Cox过程（LGCP）后，confint()返回全NA、var()/cvar()无法正常计算协方差矩阵，主要原因有三点：

最小对比度法的局限性：该方法通过匹配非齐次K函数拟合模型，不会自动计算系数的方差-协方差矩阵，因此无法直接生成置信区间。
协变量NA值干扰：拟合警告显示0.57%的积分点（quadrature points）对应协变量为NA，破坏了方差计算的数值稳定性。
查询点超出图像域：部分点超出协变量图像范围，插值引入的误差进一步影响计算。

解决方案

1. 改用似然拟合方法（优先推荐）

LGCP的似然拟合方法（method="palm"或method="ho"）会自动计算系数的协方差矩阵，支持confint()、var()等函数正常工作。修改拟合代码：

kppm_all_fes <- kppm(
  X = quad_fes,
  trend = ~log(daylight_P50sp_im) + log(daylight_Hits_im) + log(d_fes_im) + 
    SpeedLimit_im + strada_tip_im + log(P1_im) + log(E1_im),
  clusters = "LGCP", 
  method = "palm"  # 或选择 method="ho"（Hoeting-Geyer似然）
)
# 现在可正常获取置信区间
confint(kppm_all_fes)

说明：palm似然计算精度较高但对资源要求稍高；ho似然计算更快，适合大样本或复杂协变量场景。

2. 清理协变量与积分点

先处理拟合警告中的NA和域范围问题，再重新拟合：

# 1. 移除协变量为NA的积分点
Q <- quad_fes
# 提取所有积分点的协变量值
covar_df <- data.frame(
  P50 = as.numeric(Q$dummy %over% daylight_P50sp_im),
  Hits = as.numeric(Q$dummy %over% daylight_Hits_im),
  d_fes = as.numeric(Q$dummy %over% d_fes_im),
  SpeedLimit = as.numeric(Q$dummy %over% SpeedLimit_im),
  strada = as.numeric(Q$dummy %over% strada_tip_im),
  P1 = as.numeric(Q$dummy %over% P1_im),
  E1 = as.numeric(Q$dummy %over% E1_im)
)
# 筛选无NA的积分点
valid_dummy <- Q$dummy[complete.cases(covar_df)]
# 重建积分方案
quad_clean <- quadscheme(data = ppp_fes, dummy = valid_dummy)

# 2. 扩展协变量图像覆盖整个点模式窗口
daylight_P50sp_im_ext <- extend.im(daylight_P50sp_im, ppp_fes$window)
daylight_Hits_im_ext <- extend.im(daylight_Hits_im, ppp_fes$window)
# 其他协变量重复上述extend.im操作

# 用清理后的输入重新拟合
kppm_clean <- kppm(
  X = quad_clean,
  trend = ~log(daylight_P50sp_im_ext) + log(daylight_Hits_im_ext) + log(d_fes_im) + 
    SpeedLimit_im + strada_tip_im + log(P1_im) + log(E1_im),
  clusters = "LGCP",
  method = "palm"  # 仍推荐用似然方法
)

3. Bootstrap手动计算置信区间（仅当必须用mincon时）

若因业务需求必须使用最小对比度法，可通过bootstrap模拟生成置信区间：

library(boot)
# 定义bootstrap统计量函数
boot_coef_func <- function(data, idx) {
  # 重采样生成LGCP点模式
  resampled_ppp <- rLGCP(kppm_all_fes)
  # 重建积分方案
  boot_quad <- quadscheme(data = resampled_ppp, dummy = test_dummy_points)
  # 重新拟合模型
  boot_fit <- kppm(
    X = boot_quad,
    trend = ~log(daylight_P50sp_im) + log(daylight_Hits_im) + log(d_fes_im) + 
      SpeedLimit_im + strada_tip_im + log(P1_im) + log(E1_im),
    clusters = "LGCP",
    method = "mincon"
  )
  # 返回系数
  coef(boot_fit)
}

# 运行bootstrap（R为重复次数，建议至少100次）
boot_result <- boot(data = ppp_fes, statistic = boot_coef_func, R = 100)

# 对每个系数计算百分位数置信区间
for(i in 1:length(coef(kppm_all_fes))){
  cat(names(coef(kppm_all_fes))[i], "的95%置信区间：\n")
  print(boot.ci(boot_result, index = i, type = "perc"))
}

说明：bootstrap计算量较大，可根据计算资源调整重复次数。

内容的提问来源于stack exchange，提问作者Marta