如何在R的分层JAGS模型中定义残差？

分层模型的残差分为个体水平残差和组水平残差，以下是具体的定义方法和实现步骤：

1. 个体水平残差定义 #

个体残差是观测值与模型对该个体的预测值之差，直接在JAGS模型的循环块中添加定义即可。

假设你的分层模型基础结构如下：

for (i in 1:N) {
  y[i] ~ dnorm(mu[i], tau_y)
  mu[i] <- beta0 + b[group[i]] + beta1 * x[i]
}
# 组随机效应
for (j in 1:J) {
  b[j] ~ dnorm(0, tau_b)
}
# 先验分布
beta0 ~ dnorm(0, 1e-6)
beta1 ~ dnorm(0, 1e-6)
tau_y ~ dgamma(0.01, 0.01)
tau_b ~ dgamma(0.01, 0.01)

添加个体残差的定义，只需在个体循环中加入一行：

for (i in 1:N) {
  y[i] ~ dnorm(mu[i], tau_y)
  mu[i] <- beta0 + b[group[i]] + beta1 * x[i]
  # 定义个体水平残差
  res_ind[i] <- y[i] - mu[i]
}

2. 组水平残差定义 #

组水平残差反映组层面的偏离情况，通常是组随机效应与组水平预测值的差。

如果你的组随机效应是围绕总体均值的偏差（如上面的示例），组残差可以直接用随机效应本身：

for (j in 1:J) {
  b[j] ~ dnorm(0, tau_b)
  # 组水平残差（组偏离总体截距的部分）
  res_group[j] <- b[j]
}

如果组层面有协变量（比如组特征z[j]），则需要计算组随机效应与组预测值的差：

for (j in 1:J) {
  b[j] ~ dnorm(mu_b[j], tau_b)
  mu_b[j] <- gamma0 + gamma1 * z[j]
  # 定义组水平残差
  res_group[j] <- b[j] - mu_b[j]
}

3. 在R中提取残差 #

运行JAGS模型时，需要将残差变量加入待保存的参数列表：

library(rjags)
# 加载模型
model <- jags.model("your_model.jags", data = your_data_list, n.chains = 3)
# 烧录迭代
update(model, n.iter = 1000)
# 抽取样本，包含残差
samples <- coda.samples(model, 
                        variable.names = c("beta0", "beta1", "res_ind", "res_group"), 
                        n.iter = 5000)

提取残差的后验均值用于后续分析：

# 提取个体残差的后验均值
res_ind_mean <- apply(samples[[1]][, grep("res_ind", colnames(samples[[1]]))], 2, mean)
# 提取组残差的后验均值
res_group_mean <- apply(samples[[1]][, grep("res_group", colnames(samples[[1]]))], 2, mean)

4. 广义分层模型的残差注意事项 #

如果是广义线性分层模型（如泊松、logistic模型），不要用简单的观测减预测，建议使用皮尔逊残差或偏差残差：

• 泊松模型的皮尔逊残差：

  res_pearson[i] <- (y[i] - mu[i]) / sqrt(mu[i])
  

• Logistic模型的皮尔逊残差：

  res_pearson[i] <- (y[i] - mu[i]) / sqrt(mu[i] * (1 - mu[i]))
  

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