首页 文章

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

提问于
浏览
2

希望有人可以帮助我这个,因为我真的卡住了,没有找到我的编码错误!

我在JAGS中使用零膨胀泊松/负二项式GLM(没有随机效应)(使用R2Jags),并且参数估计,先验,初始值和链收敛都很好 . 所有结果都完全符合例如pscl-package的估计值,包括我在模型中计算的皮尔逊残差...

我唯一无法工作的是从模型中采样一个新样本以获得贝叶斯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

有谁知道我做错了什么?任何帮助将深表感谢!

亲切的问候,Ulf

2 回答

  • 0

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

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

    • 如果你想要一个负二项式(= gamma-Poisson)尝试从一个伽马分布中采样mu [i] - 我已经广泛使用了这个公式用于ZINB模型,所以我相信它是有效的

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

    • 小点,但你在for循环中计算aB的许多确定性函数 - 这会减慢你的模型...

    希望有所帮助,

    马特

  • 2

    因此,在搜索到我的编码错误后,在经历了一次神经衰弱并一次又一次地输入所有内容之后,我发现了我曾经做过的最愚蠢的错误 - 到目前为止:

    我只是没有指定“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

    所以我希望我没有对任何人造成重大困惑......

相关问题