关于R spatstat包中profilepl()构建复杂混合Gibbs模型的问询

问题 #

在R语言spatstat包中使用profilepl()构建包含硬核（Hardcore）和多个Geyer交互的混合Gibbs模型时遇到困难，尝试了如下代码：

RR <- c(expand.grid(r=seq(4,8, by=1), sat = 1:2), expand.grid(r1=seq(8,12,  by=1), sat1 = 1:2))

MS <- function(r, sat, r1, sat1) { Hybrid(A=Hardcore(NA), B=Geyer(r=r, sat = sat),  C=Geyer(r=r1, sat = sat1)) }

fit <- profilepl(RR, MS, swedishpines ~ polynom(x,y,2), correction = "isotropic", aic = FALSE)

存在以下疑问：

1. 当前写法是否正确？交互函数（RR和MS）是否有更优写法？
2. 能否仅提供一个数值区间（如4-12），让profilepl()自动选择两个最优Geyer交互半径，而非手动拆分两个范围？
3. 绘制轮廓似然图得到异常输出，有没有更好的可视化方法，尤其是能同时标识两个最优Geyer交互参数的方法？

解答 #

一、代码写法的正确性与优化 #

1. 当前写法的核心问题：

   • RR的构造错误：c(expand.grid(...), expand.grid(...))会把两个数据框按列拼接成长向量，完全不符合profilepl()对参数网格的要求——该函数需要数据框格式，每一行对应一组待拟合的参数组合，所有参数需在同一数据框的列中。
   • Hardcore(NA)写法不合理：Hardcore模型必须指定硬核距离，若想让硬核距离参与轮廓似然估计，需将其加入参数网格；若固定硬核距离，直接填入具体数值（如Hardcore(3)）即可。

2. 优化后的代码示例：
   正确构造包含所有参数的网格，统一管理待优化参数：

   # 构造参数网格：包含硬核距离h、两个Geyer的半径r/r1和饱和值sat/sat1
   RR <- expand.grid(h = seq(2, 4, by=1), 
                     r = seq(4,8, by=1), sat = 1:2,
                     r1 = seq(8,12, by=1), sat1 = 1:2)
   
   # 混合模型构造函数
   MS <- function(h, r, sat, r1, sat1) {
     Hybrid(Hardcore(h), Geyer(r=r, sat=sat), Geyer(r=r1, sat=sat1))
   }
   
   # 拟合模型
   fit <- profilepl(RR, MS, swedishpines ~ polynom(x,y,2), 
                    correction = "isotropic", aic = FALSE)
   

   若无需优化硬核距离，直接固定数值即可：

   RR <- expand.grid(r = seq(4,8, by=1), sat = 1:2,
                     r1 = seq(8,12, by=1), sat1 = 1:2)
   
   MS <- function(r, sat, r1, sat1) {
     Hybrid(Hardcore(3), Geyer(r=r, sat=sat), Geyer(r=r1, sat=sat1))
   }
   

二、自动选择两个最优Geyer半径的方法 #

profilepl()本身不支持直接从单一区间筛选两个不同的最优半径，因为它需要预先定义所有待测试的参数组合，但可以通过以下方式实现类似效果：

1. 先在4-12区间内生成足够多的半径候选值，再构造所有r < r1的两两组合（避免冗余的重复组合）：

   # 生成半径候选值
   r_candidates <- seq(4,12, by=1)
   # 构造所有r < r1的有效组合
   r_pairs <- expand.grid(r = r_candidates, r1 = r_candidates)
   r_pairs <- r_pairs[r_pairs$r < r_pairs$r1, ]
   
   # 结合sat参数构造完整参数网格
   RR <- expand.grid(r = r_pairs$r, r1 = r_pairs$r1,
                     sat = 1:2, sat1 = 1:2)
   

2. 拟合完成后，通过fit$maxprofile提取似然值最高的参数组合，其中就包含两个最优半径。

三、轮廓似然的可视化优化 #

默认的plot(fit)在多参数场景下易出现混乱，推荐以下几种可视化方法：

1. 提取似然值绘制热力图：
   用ggplot2提取参数与对应似然值，绘制热力图并标记最优参数点：

   library(ggplot2)
   
   # 提取参数网格与对应对数似然值
   prof_data <- data.frame(fit$pargrid, loglik = fit$profile)
   
   # 绘制两个Geyer半径的似然热力图
   ggplot(prof_data, aes(x=r, y=r1, fill=loglik)) +
     geom_tile() +
     scale_fill_viridis_c() +
     # 标记似然值最高的参数组合
     geom_point(data = prof_data[which.max(prof_data$loglik), ],
                aes(x=r, y=r1), color="red", size=3, shape=19) +
     labs(x="Geyer 1 半径r", y="Geyer 2 半径r1", fill="对数似然值") +
     theme_bw()
   

2. 按sat参数分面可视化：
   若需要关注饱和值sat的影响，可以分面展示：

   ggplot(prof_data, aes(x=r, y=r1, fill=loglik)) +
     geom_tile() +
     scale_fill_viridis_c() +
     geom_point(data = prof_data[which.max(prof_data$loglik), ],
                aes(x=r, y=r1), color="red", size=3, shape=19) +
     facet_grid(sat ~ sat1) +
     labs(x="Geyer 1 半径r", y="Geyer 2 半径r1", fill="对数似然值") +
     theme_bw()
   

3. 使用spatstat内置sliceplot：
   固定其他参数，仅查看两个半径的似然切片：

   sliceplot(fit, ~r + r1, cond = list(sat=1, sat1=1))
   

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