Как получить новые образцы из ZIP или ZINB-модели для байесовского p-значения

Надеюсь, кто-нибудь сможет мне помочь с этим, потому что я действительно застрял и не нашел свою ошибку кодирования!

Я подбираю нулевые пуассоновские / отрицательные биномиальные GLM (без случайных эффектов) в JAGS (с R2Jags), и все в порядке с оценками параметров, априорными значениями, начальными значениями и сходимостью цепочек. Все результаты полностью соответствуют, например, оценкам из пакета 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

Кто-нибудь знает, что я делаю не так? Любая помощь будет принята с благодарностью!

С уважением, Ульф


r bayesian poisson glm jags
person Woosah    schedule 20.08.2014    source источник
comment
вы используете Jags1$BUGSoutput$sims.list$Y, но это не имеет смысла, поскольку Y - это данные, а не выходной параметр, поэтому не должно быть никакого моделирования.   -  person Tomas    schedule 15.03.2015

Ответы (2)


arrow_upward
2
arrow_downward

Я подозреваю, что это не так, но единственная очевидная причина, по которой я могу подозревать, что Y [i] и YNew [i] всегда идентичны, - это то, что mu.eff [i] равно ~ нулю, либо потому, что W [i] равно 0 или mu [i] близко к нулю. Это означает, что Y [] всегда равно нулю, что легко проверить по вашим данным, но, как я уже сказал, кажется странным, что вы пытаетесь смоделировать это ... В противном случае я не уверен, что происходит. .. попробуйте упростить код, чтобы увидеть, решит ли это проблему, а затем добавляйте вещи обратно, пока он снова не сломается. Некоторые другие предложения:

  • Для отладки может быть полезно посмотреть на абсолютные значения Y и YNew, а не только на Y == YNew

  • Если вам нужен отрицательный бином (= гамма-Пуассон), попробуйте выборку mu [i] из гамма-распределения - я широко использовал эту формулировку для моделей ZINB, поэтому уверен, что она работает

  • Ваш априор для aB выглядит для меня странным - он дает предварительный 95% доверительный интервал для нулевой инфляции около 12-88% - это то, что вы планировали? И почему среднее значение 0,001, а не 0? Если у вас нет предикторов, то предварительная бета-версия для psi.min кажется более естественной, а если у вас нет полезной предварительной информации, очевидным выбором будет предшествующая бета (1,1).

  • Незначительный момент, но вы вычисляете множество детерминированных функций aB в цикле for - это замедлит вашу модель ...

Надеюсь, это поможет,

Мэтт

person Matt Denwood    schedule 04.09.2014
comment
Спасибо за ваш ответ! Я сравнил абсолютные значения Y и YNew, как вы предложили, что привело меня к ответу, потому что значения Y, которые я видел, не были целыми числами, поскольку они должны быть исходными данными ... - person Woosah; 29.09.2014
comment
И второе спасибо, потому что вы указали мне на слишком узкую приору для aB! Помещение aB в цикл просто вытекает из плана по добавлению предикторов позже, а с aB внутри цикла мне не нужно больше ничего менять, чтобы запустить расширенную модель ...! - person Woosah; 29.09.2014

arrow_upward
0
arrow_downward

Итак, после нервного срыва и набора текста снова и снова в поисках своей ошибки кодирования я обнаружил самую глупую ошибку, которую когда-либо делал - пока что:

Я просто не указал «Y» в качестве параметра для сохранения, только «YNew», поэтому, когда я сравнил YNew и Y из sims.list с all.equal, я не получил того, что, как я думал, должен. Я не знаю, почему JAGS вообще дает мне Y (из sims.list объекта JAGS), но по какой-то причине он просто дает мне YNew, когда его просят дать Y. Так что эта часть на самом деле правильная:

Jags1 $ BUGSoutput $ mean $ Fit [1] 109,7883 Jags1 $ BUGSoutput $ mean $ FitNew [1] 119,2111

Так что я надеюсь, что никого не запутал ...

person Woosah    schedule 29.09.2014
comment
Рад, что ответ помог. Обратите внимание, что на самом деле это не JAGS или rjags, которые возвращают 'Y' для YNew - это пакет r2jags (который, как я предполагаю, вы используете), находящийся между вами и JAGS, вызвал это! Это произошло из-за того, что некоторые (но не все) методы извлечения элементов списков в R выполняют частичное сопоставление - так в основном это произошло: mylist ‹- list (YNew = 'any'); mylist $ Y - с другой стороны, попытка доступа к mylist ['Y'] приводит к возвращению ожидаемого NULL. - person Matt Denwood; 01.10.2014