如何在rjags中实现三层嵌套线性混合效应模型？
你遇到的问题核心是混淆了嵌套模型和交叉模型的参数化方式，导致你的JAGS代码尝试估计不可识别的参数，最终出现个体随机效应被过度惩罚的情况。让我们一步步理清并修正这个问题：
问题根源：嵌套模型 vs 交叉模型
你提到的lme4代码lmer(yN ~ x + (1 |group/indiv), data=qq)等价于lmer(yN ~ x + (1|group) + (1|group:indiv), data=qq)，这是真正的嵌套模型：个体的随机效应是嵌套在组内的，即每个个体的截距是围绕其所属组的截距波动的。
而你写的JAGS代码对应的是lmer(yN ~ x + (1|group) + (1|indiv), data=qq)，这是交叉模型——当每个个体只属于一个组时，这个模型是不可识别的：个体的随机效应已经包含了组的效应，导致模型无法区分组和个体的独立贡献，最终出现参数被过度压缩（惩罚）的情况。
正确的JAGS嵌套模型代码
我们可以直接实现嵌套结构，不需要创建额外的交互变量，核心是让个体的随机效应以其所属组的随机效应为均值，而不是直接围绕总体均值波动。
完整JAGS模型代码
st <- " model { # 固定效应 beta0 ~ dnorm(0, 0.0001) # 总体截距 beta1 ~ dnorm(0, 0.0001) # x的斜率 # 残差精度 tau_resid ~ dgamma(0.01, 0.01) sigma_resid <- sqrt(1/tau_resid) # 组水平随机效应：每个组的截距围绕总体均值波动 for (k in 1:nGrp) { b_group[k] ~ dnorm(0, tau_group) } tau_group ~ dgamma(0.001, 0.001) sigma_group <- sqrt(1/tau_group) # 个体水平随机效应：每个个体的截距围绕其所属组的截距波动 for (j in 1:nInd) { b_indiv[j] ~ dnorm(b_group[grp_of_ind[j]], tau_indiv) } tau_indiv ~ dgamma(0.001, 0.001) sigma_indiv <- sqrt(1/tau_indiv) # 观测模型 for (i in 1:n) { mu[i] <- beta0 + beta1 * x[i] + b_indiv[ind[i]] y[i] ~ dnorm(mu[i], tau_resid) } } "
关键修改点
- 嵌套结构的参数化：个体随机效应
b_indiv[j]的先验均值不再是0，而是其所属组的随机效应b_group[grp_of_ind[j]]，完美对应(1|group/indiv)的嵌套逻辑。 - 新增映射关系：我们需要一个向量
grp_of_ind，用来存储每个个体对应的组索引（比如个体a属于组1，个体b属于组1，等等）。
运行模型的准备代码
在运行JAGS之前，先创建grp_of_ind向量：
library(rjags) library(lme4) # 处理数据集，创建个体到组的映射 qq$ind_num <- as.integer(qq$indiv) qq$group_num <- as.integer(qq$group) # 得到每个个体对应的组：取每个个体的第一个观测的组编号 grp_of_ind <- tapply(qq$group_num, qq$ind_num, function(x) unique(x)) grp_of_ind <- unname(grp_of_ind) # 整理JAGS输入数据 jags_data <- list( y = qq$yN, x = qq$x, ind = qq$ind_num, n = nrow(qq), nInd = length(unique(qq$ind_num)), nGrp = length(unique(qq$group_num)), grp_of_ind = grp_of_ind ) # 初始化并运行模型 mod <- jags.model( textConnection(st), data = jags_data, n.adapt = 100000, # 不需要1e6这么大，1e5足够 inits = list(.RNG.seed=1, .RNG.name="base::Wichmann-Hill") ) update(mod, n.iter=50000) # 额外的 burn-in mod_samples <- coda.samples( mod, variable.names=c("beta0", "beta1", "sigma_group", "sigma_indiv", "sigma_resid"), n.iter=500000, thin=10 ) # 查看结果 summary(mod_samples)
对应Pastes数据集的实现
对于Pastes数据集的嵌套模型lmer(strength ~ 1 + (1|batch/cask), data=Pastes)，我们可以用完全相同的逻辑实现，不需要创建batch:cask交互变量：
JAGS模型代码
pastes_st <- " model { # 总体截距 beta0 ~ dnorm(0, 0.0001) # 残差精度 tau_resid ~ dgamma(0.01, 0.01) sigma_resid <- sqrt(1/tau_resid) # 批次（batch）水平随机效应 for (b in 1:nBatch) { b_batch[b] ~ dnorm(0, tau_batch) } tau_batch ~ dgamma(0.001, 0.001) sigma_batch <- sqrt(1/tau_batch) # 桶（cask）水平随机效应：嵌套在批次中 for (c in 1:nCask) { b_cask[c] ~ dnorm(b_batch[batch_of_cask[c]], tau_cask) } tau_cask ~ dgamma(0.001, 0.001) sigma_cask <- sqrt(1/tau_cask) # 观测模型 for (i in 1:n) { mu[i] <- beta0 + b_cask[cask[i]] strength[i] ~ dnorm(mu[i], tau_resid) } } "
运行代码
data(Pastes, package="lme4") Pastes$batch_num <- as.integer(Pastes$batch) Pastes$cask_num <- as.integer(Pastes$cask) # 创建桶到批次的映射 batch_of_cask <- tapply(Pastes$batch_num, Pastes$cask_num, function(x) unique(x)) batch_of_cask <- unname(batch_of_cask) pastes_data <- list( strength = Pastes$strength, cask = Pastes$cask_num, n = nrow(Pastes), nBatch = length(unique(Pastes$batch_num)), nCask = length(unique(Pastes$cask_num)), batch_of_cask = batch_of_cask ) pastes_mod <- jags.model( textConnection(pastes_st), data = pastes_data, n.adapt = 10000, inits = list(.RNG.seed=2) ) update(pastes_mod, n.iter=5000) pastes_samples <- coda.samples( pastes_mod, variable.names=c("beta0", "sigma_batch", "sigma_cask", "sigma_resid"), n.iter=100000, thin=5 ) summary(pastes_samples)
为什么你的原代码会出现惩罚过度？
你的原代码同时估计组和个体相对于总体均值的随机效应，当个体完全嵌套在组内时，这两个效应是共线性的——模型无法区分“组带来的差异”和“个体带来的差异”，因此JAGS会自动压缩其中一个效应的方差（通常是个体效应，因为组的水平更少），最终出现你看到的“个体水平随机效应惩罚过度”的情况。
通过改用嵌套参数化（个体效应围绕组效应波动），我们让模型的参数变得可识别，就能得到和lme4一致的合理估计结果。
内容的提问来源于stack exchange，提问作者user2957945
