如何用spatstat的kppm获取LGCP模型scale和var参数的置信区间

问题背景 #

使用spatstat包的kppm()拟合对数高斯Cox过程（LGCP）模型时，调用confint()仅能返回固定效应（如Intercept、协变量系数）的置信区间，无法获取协方差参数scale和var的区间。示例代码：

library(spatstat)
mod <- kppm(redwood ~ x, "LGCP", model = "matern", nu = 0.5)
confint(mod)

返回结果仅包含Intercept和x的系数区间，缺失scale和var的置信区间。

你尝试过参数自助法模拟获取区间，但对合理性存疑且大数据下耗时过长：

sim_pp <- simulate(mod, nsim = 200)

mod_sim <- lapply(sim_pp, function(Y){
  kppm (Y ~ x, "LGCP", model = "matern", nu = 0.5)
})

cov_pars <- lapply(mod_sim, function(X){
  var <- X[["par"]][1]
  scale <- X[["par"]][2]
  return(data.frame(Var = var, Scale = scale))
})

cov_pars <- Reduce(rbind, cov_pars)

var_CI <- round(quantile(cov_pars$Var, probs = c(0.025, 0.975)), 6)
scale_CI <- round(quantile(cov_pars$Scale, probs = c(0.025, 0.975)), 6)

可行解决方案 #

1. 渐近正态近似法（最快）

利用模型的渐近方差-协方差矩阵，基于正态分布假设计算置信区间，大样本下可靠：

# 获取包含所有参数的方差-协方差矩阵
vcov_full <- vcov(mod, type = "vcov")
# 查看所有参数的名称及顺序
print(names(coef(mod)))
# 提取协方差参数的索引（根据实际输出调整）
cov_par_idx <- which(names(coef(mod)) %in% c("var", "scale"))
# 计算95%置信区间
cov_pars_ci <- coef(mod)[cov_par_idx] + qnorm(c(0.025, 0.975)) %o% sqrt(diag(vcov_full)[cov_par_idx])
colnames(cov_pars_ci) <- c("2.5%", "97.5%")
print(cov_pars_ci)

注意：小样本下正态假设可能不成立，结果偏差较大。

2. 优化参数自助法

你的模拟方法本身是合理的（参数自助法），可通过并行计算大幅提速：

library(parallel)
# 使用除1个核心外的所有CPU核心
n_cores <- detectCores() - 1
# 并行拟合模拟数据集
sim_pp <- simulate(mod, nsim = 100) # 可根据精度需求调整模拟次数
mod_sim <- mclapply(sim_pp, function(Y){
  kppm(Y ~ x, "LGCP", model = "matern", nu = 0.5)
}, mc.cores = n_cores)
# 提取协方差参数并计算分位数区间
cov_pars <- do.call(rbind, lapply(mod_sim, function(X){
  data.frame(Var = X[["par"]]["var"], Scale = X[["par"]]["scale"])
}))
var_CI <- round(quantile(cov_pars$Var, c(0.025, 0.975)), 6)
scale_CI <- round(quantile(cov_pars$Scale, c(0.025, 0.975)), 6)

合理性说明：参数自助法基于原模型正确设定的假设，若原模型拟合良好，区间可靠性高；若模型拟合差，结果不可靠。

3. 剖面似然法（最稳健）

剖面似然法不依赖正态假设，结果更稳健，spatstat的profile.plppm()支持此功能：

# 计算var参数的剖面似然并获取置信区间
prof_var <- profile.plppm(mod, which = "var")
ci_var <- confint(prof_var)

# 计算scale参数的剖面似然并获取置信区间
prof_scale <- profile.plppm(mod, which = "scale")
ci_scale <- confint(prof_scale)

print("var的95%置信区间：")
print(ci_var)
print("scale的95%置信区间：")
print(ci_scale)

优势：无需正态假设，对小样本更友好；计算量通常介于渐近法和自助法之间。

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