JAGS实现嵌套分组线性AR1预测模型索引越界报错求助

问题场景 #

需要在JAGS中实现带平衡嵌套分组的线性AR1模型，完成缺失值预测：

• 共2个分组，每组包含3条时序观测
• 每组第3条观测为缺失值，需利用组内前期观测预测该缺失值
• 原始代码运行时报错：Index out of range taking subset of group.

原始代码如下：

data<-data.frame(group=c(1, 1, 1, 2, 2, 2), y=c(3, 5, NA, 10, 2, NA))

model<-function(){
    for(i in 1:6){
        y[group[i]]~dnorm(mu[group[i]], tau)
        mu[group[i]]<-beta[1] + beta[2]*y[group[i-1]]
    }
    for(i in 1:2){
        beta[i]~dnorm(0, .01)
    }
    tau~dgamma(.01, .01)
}

model.data<-list("group", "y")
model.data<-list(group=data$group, y=data$y)
model.params<-c("y")

model.fit<-jags(data=model.data, inits=NULL, model.params, n.chains=2, n.iter=10000, n.burnin=1000, model.file=model)

错误原因 #

• 索引逻辑错位：将观测位置索引和分组编号混淆，y是长度为6的观测向量，索引应为1~6的观测序号，而非1/2的分组号，原始写法y[group[i]]会重复索引分组对应的y值，完全偏离数据结构。
• 循环边界越界：当循环到i=1时，i-1=0，JAGS不支持0索引，直接触发子集索引越界报错。
• 滞后项未做分组隔离：直接使用y[i-1]作为滞后项，会将上一分组的末尾观测值，错误作为下一分组首条观测的前置值，不符合组内独立的嵌套结构要求。
• 存在冗余代码：连续两次对model.data赋值，第一次赋值为无效代码。

修正代码 #

首先补充组内时序位置标识，明确每条观测在所属组内的时间顺序，再调整模型索引逻辑，将AR1滞后计算严格限制在组内：

library(R2jags)

# 补充组内时间点索引t
data <- data.frame(
  group = c(1, 1, 1, 2, 2, 2),
  t = c(1,2,3,1,2,3),
  y = c(3, 5, NA, 10, 2, NA)
)

# 修正后的嵌套AR1模型
model <- function(){
  # 观测层似然
  for(i in 1:N){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta[1] + beta[2] * y_lag[i]
  }

  # 构造组内滞后项，避免跨组值污染
  for(i in 1:N){
    # 仅组内时间点大于1的观测使用前一期值，首条观测的滞后值无实际计算意义，填充0即可
    y_lag[i] <- ifelse(t[i] > 1, y[i-1], 0)
  }

  # 参数先验
  for(j in 1:2){
    beta[j] ~ dnorm(0, 0.01)
  }
  tau ~ dgamma(0.01, 0.01)
}

# 模型输入数据
model.data <- list(
  N = nrow(data),
  t = data$t,
  y = data$y
)
# 待监控参数，缺失y值会被JAGS自动采样预测
model.params <- c("beta", "tau", "y")

# 拟合模型
model.fit <- jags(
  data = model.data,
  inits = NULL,
  parameters.to.save = model.params,
  n.chains = 2,
  n.iter = 10000,
  n.burnin = 1000,
  model.file = model
)

# 查看结果，包含两组缺失值的后验预测分布
print(model.fit)

修正后模型不会再触发索引错误，AR1自回归计算被严格限制在组内，可直接输出两组第3条缺失观测的预测结果。

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