如何在R中无需重新拟合INLA模型即可获取预测结果？

问题背景 #

当用大型数据集拟合INLA模型后，希望基于更新的协变量定期计算预测结果。常规方法是将预测用协变量加入原数据集、把结果列设为NA后重新拟合，但我们期望仅拟合一次模型并保存，后续直接调用模型进行预测（类似R中lm()搭配predict()的用法）。

实现方法 #

INLA没有原生的predict()函数，但可以通过提取已拟合模型的后验参数，结合新协变量手动计算预测分布。核心思路是：先拟合模型并保留后验样本信息，再用新协变量构造线性预测器，最后基于参数后验样本生成预测值。

可复现代码 #

library(INLA)

# 模拟数据集
n = 100; a = 1; b = 1; tau = 100
z = rnorm(n)
eta = a + b*z

scale = exp(rnorm(n))
prec = scale*tau
y = rnorm(n, mean = eta, sd = 1/sqrt(prec))
plot(z,y)

# 拟合INLA模型（保留后验配置信息）
data = list(y=y, z=z)
formula = y ~ 1+z
# 加入control.compute=list(config = TRUE)以保留生成后验样本的配置
result = inla(formula, family = "gaussian", data = data,
              control.compute=list(config = TRUE))
summary(result)

# 定义新的预测用协变量（无需加入原数据集）
new_z = seq(2,4,length.out=100)

# 提取模型的后验参数样本
post_samples <- inla.posterior.sample(n = 1000, result = result)

# 基于后验样本计算新协变量的预测值
pred_list <- lapply(post_samples, function(sample) {
  # 提取截距和z的系数
  intercept = sample$latent[rownames(sample$latent) == "(Intercept)"]
  z_coef = sample$latent[rownames(sample$latent) == "z"]
  # 计算线性预测器
  eta_pred = intercept + z_coef * new_z
  # 基于高斯分布生成预测值（考虑精度参数）
  prec_pred = sample$hyperpar[rownames(sample$hyperpar) == "Precision for the Gaussian observations"]
  rnorm(length(new_z), mean = eta_pred, sd = 1/sqrt(prec_pred))
})

# 将预测样本整理成矩阵
pred_matrix <- do.call(cbind, pred_list)
# 计算每个新协变量的预测均值和95%置信区间
pred_mean = rowMeans(pred_matrix)
pred_lower = apply(pred_matrix, 1, quantile, 0.025)
pred_upper = apply(pred_matrix, 1, quantile, 0.975)

# 可视化预测结果
plot(z, y, main = "INLA模型预测结果", xlab = "z", ylab = "y")
lines(new_z, pred_mean, col = "red", lwd = 2)
lines(new_z, pred_lower, col = "blue", lty = 2)
lines(new_z, pred_upper, col = "blue", lty = 2)
legend("topleft", legend = c("观测值", "预测均值", "95%置信区间"),
       col = c("black", "red", "blue"), lty = c(NA, 1, 2), pch = c(1, NA, NA))

关键说明 #

• 拟合模型时必须设置control.compute=list(config = TRUE)，这样才能生成后验样本。
• 无需将新预测数据加入原数据集，直接用新协变量结合后验参数计算预测值，避免重复拟合模型。
• 针对不同的分布族（如泊松、负二项等），需要调整预测值的生成逻辑（比如高斯用rnorm，泊松用rpois）。

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