基于R2jags包的JAGS模型利用GIS新数据预测新观测值的问题

我用R2jags包构建了一个JAGS模型，拟合经验数据的二次函数并估算参数，现在想给栅格图层转成的新数据框做预测，但遇到了问题：

第一次尝试的问题 #

运行「第一次尝试」的代码后，结果不符合预期，预期数据会有更多噪声波动，但实际没有。

第二次尝试的问题 #

把新数据框的响应变量y列设为NA后重新运行模型，出现如下错误：

Error in all$sims.array[, , "deviance", drop = FALSE] :
subscript out of bounds
In addition: Warning message:
In FUN(X[[i]], ...) : Failed to set trace monitor for deviance
There are no observed stochastic nodes

这是我第一次用JAGS做预测，求帮助。

可复现代码 #

library(geodata)
library(R2jags)

raw_data <- data.frame(T = c(-24, -22, -20, -19, -18, -16, -15, -12, -10 , -5, -2, -1, seq(1, 34, by = 1), 35, 37, 39, 41, 42, 43, 45, 46, 47), 
                       dev = c(0, 0, 0, 0, 0, 0, 0, 0.06, 0.6, 0.9, 0.9, 1, rep(1, length(seq(1, 34, by = 1))), c(0.8, 0.8, 0.4, 0, 0, 0, 0, 0, 0)))
## plot(raw_data$T, raw_data$dev)

## Define a JAGS model
mod <- function(){
  
  ## Priors
  q ~ dunif(0, 1)
  T0 ~ dunif(-24, -2)
  Tn ~ dunif(35, 47)
  sigma ~ dunif(0, 1000)
  tau <- 1 / (sigma * sigma)
  
  ## Likelihood
  for(i in 1:n){
    mu[i] <- -1 * q * (T[i] - T0) * (T[i] - Tn) * (Tn > T[i]) * (T0 < T[i]) * (-1 * q * (T[i] - T0) * (T[i] - Tn) < 1) + (-1 * q * (T[i] - T0) * (T[i] - Tn) > 1)
    y[i] ~ dnorm(mu[i], tau)
  }
  
}

## Group the data
data <- list(y = raw_data$dev, n = dim(raw_data)[1], T = raw_data$T)

## Define the parameters of interest to be estimated
param <- c("q", "T0", "Tn", "sigma", "mu")

## Define the initial values for the parameters of interest to be estimated
inits <- function(){list(q = 0.01, T0 = -20, Tn = 40, sigma = rlnorm(1))}

## Run the JAGS model
jags_mod <- R2jags::jags(data = data, inits = inits, parameters.to.save = param, n.thin = 8, n.chains = 3, n.burnin = 5000, n.iter = 25000, model.file = mod)

第一次尝试代码 #

## Import the raster layer
bio_r <- worldclim_country(country = "France", var="bio", res = 10, path = tempdir())
## plot(bio_r$wc2.1_30s_bio_1)

## Convert the raster layer to a data fame
T_df <- terra::as.data.frame(bio_r$wc2.1_30s_bio_1, xy = TRUE, na.rm = TRUE)
## summary(T_df)
colnames(T_df) <- c("x", "y", "T")

## Apply the quadratic function to the raster layer
q <- jags_mod$BUGSoutput$summary[row.names(jags_mod$BUGSoutput$summary) %in% "q", "mean"]
T0 <- jags_mod$BUGSoutput$summary[row.names(jags_mod$BUGSoutput$summary) %in% "T0", "mean"]
Tn <- jags_mod$BUGSoutput$summary[row.names(jags_mod$BUGSoutput$summary) %in% "Tn", "mean"]
T_df$dev <- -1 * q * (T_df$T - T0) * (T_df$T - Tn) * (Tn > T_df$T) * (T0 < T_df$T) * (-1 * q * (T_df$T - T0) * (T_df$T - Tn) < 1) + (-1 * q * (T_df$T - T0) * (T_df$T - Tn) > 1)
## summary(T_df)
plot(T_df$T, T_df$dev, pch = 16)

第二次尝试代码 #

T_df$dev <- NA
new_data <- list(y = T_df$dev, n = dim(T_df)[1], T = T_df$T)
jags_mod <- R2jags::jags(data = new_data, inits = inits, parameters.to.save = param, n.thin = 8, n.chains = 3, n.burnin = 5000, n.iter = 25000, model.file = mod)

解决方案 #

问题1：第一次尝试无噪声的原因 #

你直接用参数的均值计算了确定性的dev值，没有加入模型拟合得到的随机误差（即sigma对应的正态噪声）。要生成带噪声的预测，需要基于每个预测点的mu，从dnorm(mu, tau)中抽样。

问题2：第二次尝试报错的原因 #

当所有y都是NA时，JAGS没有观测节点可以计算似然，自然无法追踪deviance（偏差），同时模型没有可拟合的数据，参数无法被约束，导致报错。正确的做法是保留原始拟合好的参数后验分布，基于此对新数据做预测，而不是重新运行无观测数据的模型。

正确的预测步骤 #

1. 提取原始模型中所有参数的后验抽样（不只是均值）
2. 对新数据的每个点，用每一组后验参数计算mu，然后从正态分布中抽样得到带噪声的预测值
3. 可以计算预测值的均值、置信区间等统计量

示例代码：

# 提取参数的后验抽样
post_samples <- jags_mod$BUGSoutput$sims.list

# 获取新数据的T值
new_T <- T_df$T
n_new <- length(new_T)
n_samples <- nrow(post_samples$q)

# 初始化预测矩阵
pred_dev <- matrix(NA, nrow = n_samples, ncol = n_new)

# 循环每个后验样本，计算预测值
for (s in 1:n_samples) {
  q_s <- post_samples$q[s]
  T0_s <- post_samples$T0[s]
  Tn_s <- post_samples$Tn[s]
  sigma_s <- post_samples$sigma[s]
  tau_s <- 1/(sigma_s^2)
  
  # 计算每个新数据点的mu
  mu_s <- -1 * q_s * (new_T - T0_s) * (new_T - Tn_s) * 
    (Tn_s > new_T) * (T0_s < new_T) * 
    (-1 * q_s * (new_T - T0_s) * (new_T - Tn_s) < 1) + 
    (-1 * q_s * (new_T - T0_s) * (new_T - Tn_s) > 1)
  
  # 从正态分布抽样得到带噪声的预测
  pred_dev[s, ] <- rnorm(n_new, mean = mu_s, sd = sigma_s)
}

# 计算每个预测点的均值、95%置信区间
T_df$dev_mean <- colMeans(pred_dev)
T_df$dev_lower <- apply(pred_dev, 2, quantile, 0.025)
T_df$dev_upper <- apply(pred_dev, 2, quantile, 0.975)

# 绘制带置信区间的预测图
plot(T_df$T, T_df$dev_mean, pch = 16, ylim = range(pred_dev), main = "带噪声的预测结果")
points(T_df$T, pred_dev[sample(n_samples, 1), ], col = "gray", pch = 16) # 随机展示一组抽样预测
segments(T_df$T, T_df$dev_lower, T_df$T, T_df$dev_upper, col = "blue")

补充说明 #

• 如果你不需要带噪声的预测，只需要拟合的趋势（即mu的后验均值），那第一次尝试的代码是对的，但你预期的噪声其实是模型的随机误差项，需要额外抽样加入。
• 不要用全NA的y重新运行JAGS模型，正确的预测应该基于已拟合好的后验分布来生成，这样既高效又符合贝叶斯预测的逻辑。

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