r - 如何从 ZIP 或 ZINB 模型获取贝叶斯 p 值的新样本

Question

希望有人可以帮助我解决这个问题，因为我真的很困惑并且没有发现我的编码错误!

我正在 JAGS(使用 R2Jags)中拟合零膨胀泊松/负二项式 GLM(无随机效应)，参数估计、先验、初始值和链收敛一切都很好。所有结果都完全符合，例如 pscl-package 的估计，包括我对模型中 PIL 逊残差的计算...

我唯一无法开始工作的是从模型中采样一个新样本以获得用于评估模型拟合的贝叶斯 p 值。我之前拟合的“正常”泊松和负二项式模型都给出了预期的重复样本，并且没有出现问题。

这是到目前为止我的代码，但重要的部分是“#New Samples”:

model{
# 1. Priors
beta  ~ dmnorm(b0[], B0[,])   
aB    ~ dnorm(0.001, 1)

    #2. Likelihood function
    for (i in 1:N){  

    # Logistic part
    W[i]           ~ dbern(psi.min1[i])
    psi.min1[i]   <- 1 - psi[i]
    eta.psi[i]    <- aB
    logit(psi[i]) <- eta.psi[i]

    # Poisson part
    Y[i]           ~ dpois(mu.eff[i])
    mu.eff[i]     <- W[i] * mu[i]
    log(mu[i])    <- max(-20, min(20, eta.mu[i]))
    eta.mu[i]     <- inprod(beta[], X[i,])

    # Discrepancy measures:
    ExpY[i]       <- mu [i] * (1 - psi[i])
    VarY[i]       <- (1- psi[i]) * (mu[i] + psi[i] * pow(mu[i], 2))
    PRes[i]       <- (Y[i] - ExpY[i]) / sqrt(VarY[i])
    D[i]          <- pow(PRes[i], 2)

    # New Samples:
    YNew[i]        ~ dpois(mu.eff[i])
    PResNew[i]    <- (YNew[i] - ExpY[i]) / sqrt(VarY[i])
    DNew[i]       <- pow(PResNew[i], 2)
    } 
Fit         <- sum(D[1:N])
FitNew      <- sum(DNew[1:N])
}

最大的问题是，我确实尝试了我认为可以/应该起作用的所有组合和更改，但是当我查看模拟样本时，我在这里得到了这个:

> all.equal( Jags1$BUGSoutput$sims.list$YNew, Jags1$BUGSoutput$sims.list$Y )

[1] TRUE

而且，更奇怪的是，当使用 Fit 和 FitNew 方法时:

> Jags1$BUGSoutput$mean$Fit
[1] 109.7883
> Jags1$BUGSoutput$mean$FitNew
[1] 119.2111

有人知道我做错了什么吗？任何帮助将不胜感激!

亲切的问候，乌尔夫

Answer 1

我怀疑情况并非如此，但我怀疑 Y[i] 和 YNew[i] 始终相同的唯一明显原因是 mu.eff[i] 是否为零，或者因为 W[i] ] 为 0 或 mu[i] 接近于零。这意味着 Y[] 始终为零，这很容易从您的数据中检查，但正如我所说，您尝试对此进行建模似乎很奇怪......否则，我不确定发生了什么。 ..尝试简化代码以查看是否可以解决问题，然后重新添加内容，直到它再次出现问题。其他一些建议:

• 查看 Y 和 YNew 的绝对值而不仅仅是 Y==YNew 可能有助于调试

• 如果您想要负二项式(= gamma-Poisson)，请尝试从 gamma 分布中采样 mu[i] - 我已在 ZINB 模型中广泛使用此公式，因此我确信它有效

  <

• 您对 aB 的先验对我来说看起来很奇怪 - 它给出了 12-88% 左右的零通胀的先验 95% CI - 这是您的意图吗？为什么平均值是 0.001 而不是 0？如果您没有预测变量，那么 psi.min 的 beta 先验似乎更自然 - 如果您没有有用的先验信息，beta(1,1) 先验将是一个明显的选择。

• 小问题，但您正在 for 循环内计算 aB 的许多确定性函数 - 这会减慢您的模型...

希望有帮助，

马特