使用simmr带先验运行稳定同位素混合模型时sigma_shape变量报错

问题原因 #

当通过prior_control自定义比例先验时，会覆盖simmr的默认参数列表，导致残差标准差的先验参数sigma_shape和sigma_rate未被传入JAGS模型，触发编译错误。

解决方案 #

1. 更新simmr到最新版本

先尝试更新包，确保使用最新稳定版：

install.packages("simmr")
library(simmr)

2. 手动指定sigma_shape和sigma_rate的默认值

若更新后问题仍存在，在prior_control中显式添加这两个参数。根据simmr的默认设置，这两个参数的默认取值为sigma_shape=2、sigma_rate=1。修改模型运行代码如下：

simmr_out_informative <- simmr_mcmc(simmr_in,
                                    prior_control =
                                      list(
                                        means = prior$mean,
                                        sd = prior$sd,
                                        sigma_shape = 2,
                                        sigma_rate = 1
                                      ))

验证修改后的完整代码 #

将上述修改应用到示例代码中，完整可运行代码如下：

library(simmr)

mix <- matrix(c(
  -10.13, -10.72, -11.39, -11.18, -10.81, -10.7, -10.54,
  -10.48, -9.93, -9.37, 11.59, 11.01, 10.59, 10.97, 11.52, 11.89,
  11.73, 10.89, 11.05, 12.3
), ncol = 2, nrow = 10)
colnames(mix) <- c("d13C", "d15N")
s_names <- c("Zostera", "Grass", "U.lactuca", "Enteromorpha")
s_means <- matrix(c(-14, -15.1, -11.03, -14.44, 3.06, 7.05, 13.72, 5.96), ncol = 2, nrow = 4)
s_sds <- matrix(c(0.48, 0.38, 0.48, 0.43, 0.46, 0.39, 0.42, 0.48), ncol = 2, nrow = 4)
c_means <- matrix(c(2.63, 1.59, 3.41, 3.04, 3.28, 2.34, 2.14, 2.36), ncol = 2, nrow = 4)
c_sds <- matrix(c(0.41, 0.44, 0.34, 0.46, 0.46, 0.48, 0.46, 0.66), ncol = 2, nrow = 4)

simmr_in <- simmr_load(
  mixtures = mix,
  source_names = s_names,
  source_means = s_means,
  source_sds = s_sds,
  correction_means = c_means,
  correction_sds = c_sds,
)

proportion_means <- c(0.4, 0.3, 0.2, 0.1)
proportion_sds <- c(0.08, 0.02, 0.01, 0.02)
prior <- simmr_elicit(
  4, proportion_means,
  proportion_sds
)

# 修改后的simmr_mcmc调用，添加sigma参数
simmr_out_informative <- simmr_mcmc(simmr_in,
                                    prior_control =
                                      list(
                                        means = prior$mean,
                                        sd = prior$sd,
                                        sigma_shape = 2,
                                        sigma_rate = 1
                                      ))

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