2

希望有人能帮我解决这个问题，因为我真的被卡住了，没有发现我的编码错误！

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

我唯一不能开始工作的是从模型中抽取一个新样本以获得贝叶斯 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

有谁知道我做错了什么？任何帮助将不胜感激！

亲切的问候，乌尔夫

4

2 回答 2

2

我怀疑情况并非如此，但我可以怀疑 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 的许多确定性函数-这会减慢您的模型速度...

希望有帮助，

马特

于 2014-09-04T06:59:31.653 回答
0

因此，在精神崩溃并在搜索我的编码错误时一次又一次地输入所有内容之后，我发现了我犯过的最愚蠢的错误 - 到目前为止：

我只是没有指定“Y”作为要保存的参数，只有“YNew”，所以当我将 sims.list 中的 YNew 和 Y 与 all.equal 进行比较时，我没有得到我认为应该的结果。我不知道为什么 JAGS 给了我 Y（来自 JAGS 对象的 sims.list），但由于某种原因，当被要求给 Y 时它只是给了我 YNew。所以这部分实际上是正确的：

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

所以我希望我没有给任何人造成严重的困惑......

于 2014-09-29T15:59:17.723 回答