如何在R2jags中实现固定Power prior模型（伯努利分布）

我们要实现一个固定Power Prior模型，事件服从参数为$p_i$的伯努利分布，需要将伯努利似然函数取$w$次幂。尝试直接在JAGS中写dbern(p[i])^w但无法识别，现已有可行的Stan代码，需要对应的JAGS（R2jags）实现。

错误的JAGS尝试代码 #

mod1<-function(){  
  for(i in 1 : N) {
    y[i] ~ dbern(p[i])^w # 此处JAGS无法识别语法
    logit(p[i]) <- a[group[i]] + t[group[i]]*d[i]
  }  
  for (j in 1:3) {
    a[j] ~ dnorm(mu.a, pow(sd.a,-2))
    t[j] ~ dnorm(mu.t, pow(sd.t,-2))
  }  
  mu.a~dnorm(0,0.01)
  mu.t~dnorm(0,0.01)
  sd.a ~ dunif(0,10)
  sd.t ~ dunif(0,10)
}

可行的Stan代码 #

mod1 <- '
data {
  int<lower = 1> N;
  real w;
  int y[N];
  vector[N] d;
  int group[N];
}
parameters {
  vector[3] a;
  vector<lower=0>[3] t;
  real mua;
  real mut;
  real sigmaa;
  real sigmat;
}
model {
  real p[N];
  // 先验分布
  mua~normal(0,10);
  mut~normal(0,10);
  sigmaa ~ uniform(0,10);
  sigmat ~ uniform(0,10);
  for (j in 1:3) {
    a[j] ~ normal(mua, sigmaa);
    t[j] ~ normal(mut, sigmat);
  }
  
  // 似然函数（取w次幂）
  for (i in 1:N) {
    p[i]=inv_logit(a[group[i]] + t[group[i]]*d[i]);
    target += bernoulli_lpmf(y[i]|(p[i]))*w;
  }
}
'

测试样本数据 #

y <- rbinom(160, 1, 0.3)
d <- runif(160, 0, 41)
group <- sample(1:3, 160, replace = TRUE)
w <- rep(0.5, 160)
N=160
mydata <- data.frame(y, d, group, w, N) # 移除原代码中未定义的x4

JAGS等效实现代码 #

JAGS不支持直接对分布取幂，但可以通过手动构造对数似然并累加到后验分布的方式实现，逻辑和Stan的target +=一致：

mod1 <- function(){
  # 初始化对数似然累加器
  loglik <- 0
  
  for(i in 1 : N) {
    logit(p[i]) <- a[group[i]] + t[group[i]]*d[i]
    # 手动计算伯努利对数似然，乘以权重w[i]后累加
    loglik <- loglik + w[i] * (y[i] * log(p[i]) + (1 - y[i]) * log(1 - p[i]))
  }
  
  # 将构造的对数似然加入后验（dflat为平坦分布，不改变后验结构，仅实现累加效果）
  loglik ~ dflat()
  
  # 先验分布部分
  for (j in 1:3) {
    a[j] ~ dnorm(mu.a, pow(sd.a, -2))
    t[j] ~ dnorm(mu.t, pow(sd.t, -2))
  }
  mu.a ~ dnorm(0, 0.01)
  mu.t ~ dnorm(0, 0.01)
  sd.a ~ dunif(0, 10)
  sd.t ~ dunif(0, 10)
}

关键说明 #

1. JAGS没有类似Stan的target机制，因此手动计算每个样本的对数似然，乘以对应权重后累加，再通过loglik ~ dflat()将累加值注入模型的总对数后验中。
2. 测试数据中w是长度为N的向量，代码中用w[i]匹配每个样本的权重。
3. 调用R2jags时，需确保传入的数据包含w变量，即使用整理后的mydata。

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