2

私は本当に立ち往生していて、コーディングエラーが見つからないので、誰かがこれで私を助けてくれることを願っています!

JAGS（R2Jagsを使用）でゼロ膨張ポアソン/負の二項GLM（ランダム効果なし）をフィッティングしていますが、パラメーター推定、事前確率、初期値、およびチェーン収束ですべて問題ありません。すべての結果は、たとえば、モデルのピアソン残差の計算を含む、pscl-パッケージからの推定値と完全に一致しています...

私が作業できない唯一のことは、モデルから新しいサンプルをサンプリングして、モデルの適合性を評価するためのベイジアン 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] が常に同一であると推測できる唯一の明らかな理由は、W[i] が 0 であるため、mu.eff[i] が ~zero である場合です。または mu[i] がゼロに近い。これは、Y[] が常にゼロであることを意味します。これは、データから簡単に確認できますが、前述したように、これをモデル化しようとするのは奇妙に思えます...そうでなければ、何が起こっているのかわかりません. .. コードを単純化して問題が解決するかどうかを確認してから、再び壊れるまで何かを追加してください。その他の提案:

  • Y==YNew だけでなく、Y と YNew の絶対値を調べると、デバッグに役立つ場合があります。

  • 負の二項式 (= ガンマ ポアソン) が必要な場合は、ガンマ分布から mu[i] をサンプリングしてみてください。

  • あなたの aB の事前確率は私には奇妙に見えます - それは 12-88% 前後のゼロインフレに対する事前の 95% CI を与えます - それはあなたが意図したものですか? そして、なぜ平均が 0 ではなく 0.001 なのですか? 予測子がない場合は、psi.min のベータ事前分布がより自然に見えます。また、有用な事前情報がない場合は、ベータ (1,1) 事前分布が明白な選択になります。

  • マイナーポイントですが、forループ内でaBの多くの決定論的関数を計算しています-これはモデルの速度を低下させます...

それが役立つことを願って、

マット

于 2014-09-04T06:59:31.653 に答える
0

それで、神経衰弱になり、コーディングエラーを探しながら何度も何度も入力した後、これまでに犯した中で最もばかげたエラーを見つけました-これまでのところ：

保存するパラメーターとして「Y」を指定せず、「YNew」のみを指定しただけなので、sims.list の YNew と Y を all.equal と比較すると、思ったとおりの結果が得られませんでした。JAGS が (JAGS オブジェクトの sims.list から) Y を与える理由はまったくわかりませんが、何らかの理由で、Y を与えるように求められたときに YNew を与えるだけです。したがって、この部分は実際には正しいです:

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

だから私は誰にも大きな混乱を引き起こさなかったことを願っています...

于 2014-09-29T15:59:17.723 に答える