RJags报错：计算日巢存活率时x索引越界问题求助

背景 #

正在构建覆盖4年的日巢存活率零模型，仅拟合截距项，采用JAGS工具通过贝叶斯定理与MCMC生成参数后验分布。模型参考Kery和Royle《模型构建策略》第5.15章最简版本，基于单年度数据分析代码修改。运行模型时触发索引越界错误，已核对first、lastCheck与矩阵x的长度和数值范围，仍无法定位原因。

索引定义 #

first <- nest_filtered$first 
# 计算方式：dt_found - earliest_dt_found + 1，共366个值，均小于42

lastCheck <- nest_filtered$lastCheck 
# 计算方式：dt_lastCheck - earliest_dt_found，共366个值，均小于42

x <- encounter_matrix 
# 366行42列的矩阵，记录巢穴每日存活状态（1=存活，0=死亡）

模型代码 #

R调用代码 #

null <- list(x=x, nest_row=nrow(x), first=first, lastCheck=lastCheck)

# 写入JAGS模型文件
writeLines("model {  
  # 固定效应先验
  intercept ~ dunif(-5, 5)  # 截距项的先验分布
  
  # 模型主体
  for (i in 1:nest_row) {  # 遍历每个巢穴
    for (j in (first[i]):lastCheck[i]) {  # 遍历巢穴被观测的每一天
      lq[i, j] <- intercept  # 仅截距项模型
      q[i, j] <- 1 / (1 + exp(-lq[i, j]))  # logit转换为日存活概率
      p[i, j] <- x[i, j-1] * q[i, j]  # 当日观测的存活概率（依赖前一日状态）
      x[i, j] ~ dbern(p[i, j])  # 观测数据服从伯努利分布
    }
  }
  
  # 衍生变量
  realsurvival <- 1 / (1 + exp(-intercept))  # 实际日存活概率
  realcumsurvival <- pow(realsurvival, 29) # 29天累计存活概率（从初始到孵化的平均时长）
}", con = "model.txt")

mNULL <- 'model.txt'

# 初始值函数
initsN <- function(){  
  list(intercept=rnorm(1,0,0.2)) # 从正态分布(0,0.2)中随机抽取截距初始值
}

# 需要估计的参数
parametersN <- c("intercept","realsurvival","realcumsurvival")

# 运行JAGS模型
outN <- jags(null,            
            initsN,            
            parameters.to.save=parametersN,            
            model.file=mNULL,            
            n.thin=1,           
            n.chains=2,           
            n.burnin=1000,           
            n.iter=22000)

报错信息 #

Error in rjags::jags.model(file = model.file, data = data, inits = inits,  :  
  RUNTIME ERROR:Compilation error on line 11.Index out of range taking subset of  x

报错对应JAGS模型的第11行：

p[i, j] <- x[i, j-1] * q[i, j]

错误原因与解决方案 #

核心原因 #

JAGS使用1-based索引，当first[i] = 1时，内层循环的起始j=1，此时j-1=0，x[i,0]超出矩阵的列索引范围（矩阵x的列索引从1到42）。

另外，巢存活模型的逻辑中，首次发现巢穴的当日（j=first[i]），巢穴必然是存活的，这一状态是已知的，不需要纳入模型拟合。

修复步骤 #

1. 调整内层循环的起始值：将JAGS模型中的内层循环从(first[i]):lastCheck[i]改为(first[i]+1):lastCheck[i]，确保j-1的最小值为first[i]（合法索引）。修改后的模型片段：

for (i in 1:nest_row) {  
  for (j in (first[i]+1):lastCheck[i]) {  # 从首次观测的次日开始拟合
    lq[i, j] <- intercept  
    q[i, j] <- 1 / (1 + exp(-lq[i, j]))  
    p[i, j] <- x[i, j-1] * q[i, j]  
    x[i, j] ~ dbern(p[i, j])  
  }
}

2. 过滤无效数据（可选）：如果存在first[i] == lastCheck[i]的巢穴（仅被观测一次，无后续存活状态数据），可提前从数据中过滤，避免空循环：

nest_filtered <- nest_filtered[nest_filtered$first < nest_filtered$lastCheck, ]
# 同步更新x、first、lastCheck
x <- x[rownames(nest_filtered), ]
first <- nest_filtered$first
lastCheck <- nest_filtered$lastCheck

3. 验证索引范围：再次确认first的最小值，确保first[i]+1 <= lastCheck[i]对所有巢穴成立，避免循环起始值大于终止值的情况。

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