在JAGS模型内从sims.list创建派生参数的可行性咨询

云原生阿岚

2026-7-5

问题描述

我使用JAGS开展整合模型/种群生存力分析，需创建派生参数用于种群预测。当前无报错模型中，依赖时间的生存参数由重捕标记模块生成，用于有数据年份的种群状态空间模型；预测年份则使用该参数的均值作为派生参数，但我认为该均值的方差过小，影响了预测种群规模。因此尝试从时间依赖生存参数的sims.list中随机抽样生成派生参数以更好地捕捉方差，但在模型内引用sims.list时出现$附近的语法错误。现咨询：能否在JAGS模型运行过程中从其sims.list中抽样？还是仅能在模型运行完成后访问该列表？

报错信息

Error in checkForRemoteErrors(val) : 3 nodes produced errors; first error: Error parsing model file: syntax error on line 103 near "$"

简化示例代码

# Load woodchat shrike data and produce data overview
library(IPMbook)
data(woodchat10)
str(woodchat10)

K <- 15    # Number of forecast years
jags.data <- list(marr.j=woodchat10$marr.j, marr.a=woodchat10$marr.a, n.occasions=ncol(woodchat10$marr.j), rel.j=rowSums(woodchat10$marr.j), rel.a=rowSums(woodchat10$marr.a), J=woodchat10$J, B=woodchat10$B, count=woodchat10$count, pNinit=dUnif(1, 50), K=K)
str(jags.data)


K <- 15    # Number of forecast years
jags.data <- list(marr.j=woodchat10$marr.j, marr.a=woodchat10$marr.a, n.occasions=ncol(woodchat10$marr.j), rel.j=rowSums(woodchat10$marr.j), rel.a=rowSums(woodchat10$marr.a), J=woodchat10$J, B=woodchat10$B, count=woodchat10$count, pNinit=dUnif(1, 50), K=K)
str(jags.data)


# Write JAGS model file
cat(file="model1.txt", "
model { 
# Priors and linear models
mean.logit.sj <- logit(mean.sj)
mean.sj ~ dunif(0, 1) 
mean.logit.sa <- logit(mean.sa)
mean.sa ~ dunif(0, 1) 
mean.p ~ dunif(0, 1)
mean.log.f <- log(mean.f)
mean.f ~ dunif(0, 10)

for (t in 1:(n.occasions-1)){
  p[t] <- mean.p
}

for (t in 1:(n.occasions-1+K)){  # Here we extend the loop to K more years
  logit.sj[t] <- mean.logit.sj + eps.sj[t]
  eps.sj[t] ~ dnorm(0, tau.sj)
  sj[t] <- ilogit(logit.sj[t])
  logit.sa[t] <- mean.logit.sa + eps.sa[t]
  eps.sa[t] ~ dnorm(0, tau.sa)
  sa[t] <- ilogit(logit.sa[t])
}

for (t in 1:(n.occasions+K)){   # Extended loop also here
  log.f[t] <- mean.log.f + eps.f[t]
  eps.f[t] ~ dnorm(0, tau.f)
  f[t] <- exp(log.f[t])
}

sigma.sj ~ dunif(0, 10)
tau.sj <- pow(sigma.sj, -2)
sigma.sa ~ dunif(0, 10)
tau.sa <- pow(sigma.sa, -2)
sigma.f ~ dunif(0, 10)
tau.f <- pow(sigma.f, -2)

sigma ~ dunif(0.5, 50)
tau <- pow(sigma, -2)

# Population count data (state-space model)
# Model for the initial population size: uniform priors
N[1,1] ~ dcat(pNinit)
N[2,1] ~ dcat(pNinit)

# Process model over time: our model of population dynamics
for (t in 1:(n.occasions-1)){  # Note extended loop
  N[1,t+1] ~ dpois(sj[t] * f[t] * (N[1,t] + N[2,t])) 
  N[2,t+1] ~ dbin(sa[t], (N[1,t] + N[2,t]))
}

for (t in n.occasions:(n.occasions-1+K)){  # Note extended loop
  N[1,t+1] ~ dpois(sj.f * f[t] * (N[1,t] + N[2,t])) 
  N[2,t+1] ~ dbin(sa[t], (N[1,t] + N[2,t]))
}

# Observation model
for (t in 1:n.occasions){
  count[t] ~ dnorm(N[1,t] + N[2,t], tau)
}

# Productivity data (Poisson regression model)
for (t in 1:n.occasions){
  J[t] ~ dpois(f[t] * B[t])
}

# Capture-recapture data (CJS model with multinomial likelihood)
# Define the multinomial likelihood
for (t in 1:(n.occasions-1)){
  marr.j[t,1:n.occasions] ~ dmulti(pr.j[t,], rel.j[t])
  marr.a[t,1:n.occasions] ~ dmulti(pr.a[t,], rel.a[t])
}
# Define the cell probabilities of the m-arrays
# Main diagonal
for (t in 1:(n.occasions-1)){
  q[t] <- 1-p[t]   # Probability of non-recapture
  pr.j[t,t] <- sj[t] * p[t]
  pr.a[t,t] <- sa[t] * p[t]
  # Above main diagonal
  for (j in (t+1):(n.occasions-1)){
    pr.j[t,j] <- sj[t] * prod(sa[(t+1):j]) * prod(q[t:(j-1)]) * p[j]
    pr.a[t,j] <- prod(sa[t:j]) * prod(q[t:(j-1)]) * p[j]
  } #j
  # Below main diagonal
  for (j in 1:(t-1)){
    pr.j[t,j] <- 0
    pr.a[t,j] <- 0
  } #j
} #t
# Last column: probability of non-recapture
for (t in 1:(n.occasions-1)){
  pr.j[t,n.occasions] <- 1-sum(pr.j[t,1:(n.occasions-1)])
  pr.a[t,n.occasions] <- 1-sum(pr.a[t,1:(n.occasions-1)])
}

# Derived parameters
# Total population size
for (t in 1:(n.occasions+K)){
  Ntot[t] <- N[1,t] + N[2,t]
}

# Random draw of survival from 1 to n.occassions-1
sj.f <- sample(out1$sims.list$sj, 1)

# Check whether the population is extinct in the future
for (t in 1:K){
  extinct[t] <- equals(Ntot[n.occasions+t], 0)  
}
} # end model
")

# Initial values
inits <- function(){list(mean.sj=runif(1, 0, 0.5), mean.sa=runif(1, 0.4, 0.6), mean.f=runif(1, 1.3, 2))}

# Parameters monitored
parameters <- c("mean.sj", "sigma.sj", "mean.sa", "sigma.sa", "mean.p", "mean.f", "sigma.f", "sigma", "sj", "sa", "f", "N", "Ntot", "extinct")

# MCMC settings
ni <- 20000; nb <- 5000; nc <- 3; nt <- 3; na <- 1000

# Call JAGS (ART 1 min), check convergence and summarize posteriors
out1 <- jags(jags.data, inits, parameters, "model1.txt", n.iter=ni, n.burnin=nb, n.chains=nc, n.thin=nt, n.adapt=na, parallel=TRUE)

解决方案

JAGS运行过程中无法访问sims.list：sims.list是JAGS运行完毕后返回的R环境对象，JAGS模型内部仅能使用模型自身定义的变量、输入数据以及JAGS内置函数，无法调用R的对象或语法（比如$访问列表元素、R版sample()）。你模型里的sj.f <- sample(out1$sims.list$sj, 1)是R代码，放到JAGS模型定义中必然触发语法错误。
要保留预测期生存参数的时间变异，有两种可行方案：
1. 沿用现有随机结构生成预测值：你模型中sj[t]已采用logit.sj[t] <- mean.logit.sj + eps.sj[t]的正态分布结构，预测年份的sj可以直接沿用该结构（你的代码中已经把循环扩展到了n.occasions-1+K），无需额外定义sj.f，JAGS会在MCMC迭代中自然生成带方差的预测值。
2. Bootstrap抽样观测年份的生存参数：如果想直接从已有观测年份的sj中随机抽取，可在JAGS模型内部通过离散分布实现：
  - 定义一个索引变量，从观测年份的索引中均匀抽样；
  - 用该索引提取对应的sj值作为预测期的生存参数。
方案2的模型关键修改示例：

# 替换原模型中sj.f的定义部分
# Derived parameters
# Total population size
for (t in 1:(n.occasions+K)){
  Ntot[t] <- N[1,t] + N[2,t]
}

# Bootstrap抽样观测年份的sj
idx ~ dcat(rep(1/(n.occasions-1), n.occasions-1))
sj.f <- sj[idx]

# Check whether the population is extinct in the future
for (t in 1:K){
  extinct[t] <- equals(Ntot[n.occasions+t], 0)  
}

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

在JAGS模型内从sims.list创建派生参数的可行性咨询

报错信息 #

简化示例代码 #

报错信息

简化示例代码