Я пытаюсь реализовать несколько различных моделей подсчета в JAGS, но с зависимой переменной, принимающей большие положительные значения. В то время как я также изучаю варианты с непрерывной зависимой переменной, в настоящее время я пытаюсь хорошо смешать модель пуассона-гамма (Rhat <1,1), даже после запуска 3 параллельных цепочек для 50 000 приработок и 50 000 итераций впоследствии. Для сравнения, достижение требуемых значений Rhat при использовании dnegbin не является проблемой. По сравнению с обоими, логнормальный пуассон работает быстрее всех (даже в моем наборе эмпирических данных). </p>
Насколько я понимаю, сравнение значений DI C NB2 и логнормального Пуассона не имеет смысла из-за разницы в параметризации. В результате, использование модели Пуассона-Гамма позволит провести соответствующее сравнение моделей. Поэтому любые мысли о том, что может повлиять на вычисления, были бы очень признательны.
Я делюсь кодом для моделирования отрицательной биномиальной зависимой переменной с ее средним значением как функцией двух ковариат и предварительно определенного параметра размера из 1.2. Кроме того, описание модели для Пуассона-Гамма (PG), NB2 и Логнормального Пуассона (PLN) также предоставляется вместе с результатами моделирования.
Моя функция инициализации рандомизирует оценки коэффициентов с использованием стандартной нормальной распределение, в то время как параметры дисперсии инициализируются с помощью rexp (1). Я помню, как читал в руководстве JAGS (v4.3) о выборке гамма-распределения для значений меньше 1, иногда могут возникать проблемы, но я не получаю явных предупреждений. Я попытался использовать некоторое усечение, но это не помогло в моделировании исследования.
library(arm)
# library(rjags)
# library(runjags)
# library(R2jags)
library(jagsUI)
library(dplyr)
library(Metrics)
##Simulating Data
set.seed(456)
b<-c(3,2,1)
sig<-list(c(0.1,0,0),c(0,1,0),c(0,0,1))
X<-mvrnorm(n=200,mu=c(1,2,3),diag(c(0.01,1,1)))
X[,1]<-1
mu<-X %*%b %>% unlist() %>% exp()
ynb<-rnegbin(200,mu,theta=1.2)
X<-as.data.frame(X)
#> summary(ynb)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 55 3892 17964 225909 115116 6760084
###########Poisson Gamma Spec###########
cat("model{
## Likelihood
for(i in 1:n){
ynb[i] ~ dpois(theta[i])
log(theta[i])<-lambda[i] + log(e[i])
lambda[i] <- inprod(b,X[i,])
e[i] ~ dgamma(phi, phi)
}
## Priors
for (k in 1:K){
b[k] ~ dnorm(0,0.001)
}
phi ~ dgamma(0.01,0.01) #T(0.1,5)
}",file="pg.jag")
n<-200
K<-3
data <- list ("n", "K", "ynb", "X")
pg.inits <- function (){
list (b=rnorm(K),phi = rexp(1))
#list (b=rnorm(K), phi = runif(1,1,2))
}
pg.parameters <- c ("b","phi")
# pedexp.pg <- autojags(model.file = "pedexp.pg.jag",pedexp.data, pedexp.pg.inits,
# pedexp.pg.parameters,parallel=T,n.chains=3,
# iter.increment=2000,n.burnin=20000,Rhat.limit=1.1,
# max.iter=40000)
pg <-jags(model.file = "pg.jag",data, pg.inits,
pg.parameters,n.chains=3, n.iter=100000,n.burnin = 50000,
parallel=T,n.thin=2,verbose = T)
summary(pg)
print(pg)
###Results###
# mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
# b[1] 3.862 0.323 3.403 3.738 4.438 FALSE 1 6.019 3
# b[2] 1.563 0.134 1.348 1.646 1.691 FALSE 1 10.434 3
# b[3] 0.969 0.041 0.904 0.963 1.048 FALSE 1 3.418 4
# phi 0.977 0.139 0.729 0.984 1.238 FALSE 1 2.426 4
# deviance 2535.477 19.958 2498.423 2534.776 2576.344 FALSE 1 1.000 75000#
# pD = 199.2 and DIC = 2734.652
###########Negative Binomial Spec###########
cat("model{
## Likelihood
for(i in 1:n){
ynb[i] ~ dnegbin(p[i], size)
p[i] <- size / (size + mu[i])
log(mu[i]) <- max(-20, min(eta[i],30))
eta[i]<-inprod(b,X[i,])
}
## Priors
for (k in 1:K){
b[k] ~ dnorm(0,0.0001)
}
size ~ dgamma(0.01,0.01) #dunif(0.1, 100)
}",file="nb.jag")
nb.inits<-function(){
list (b=rnorm(K),size = rexp(1))
}
nb.parameters <- c ("b","size")
nb <-jags(model.file = "nb.jag",data, nb.inits,
nb.parameters,n.chains=3, n.iter=100000,n.burnin = 50000,
parallel=T,n.thin=2,verbose = T)
summary(nb)
print (nb)
###Results###
# mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
# b[1] 2.940 0.231 2.490 2.940 3.395 FALSE 1 1.000 17359
# b[2] 2.034 0.068 1.903 2.034 2.168 FALSE 1 1.001 4007
# b[3] 1.012 0.061 0.893 1.012 1.131 FALSE 1 1.000 34237
# size 1.265 0.115 1.051 1.261 1.499 FALSE 1 1.000 75000
# deviance 4508.483 2.850 4504.912 4507.848 4515.640 FALSE 1 1.000 22729
# pD = 4.1 and DIC = 4512.542
#########Poisson Lognormal Spec##########
cat("model{
## Likelihood
for(i in 1:n){
ynb[i] ~ dpois(theta[i])
log(theta[i])<-lambda[i] + eps[i]
lambda[i] <- inprod(b,X[i,])
eps[i] ~ dnorm(0,tau)
}
## Priors
for (k in 1:K){
b[k] ~ dnorm(0,0.001)
}
tau <- pow(sigma, -2)
sigma ~ dunif (0, 100)
}",file="pln.jag")
pln.inits<-function(){
list (b=rnorm(K),sigma = rexp(1))
}
pln.parameters <- c ("b","sigma")
pln <-jags(model.file = "pln.jag",data, pln.inits,
pln.parameters,n.chains=3, n.iter=100000,n.burnin = 50000,
parallel=T,n.thin=2,verbose = T)
summary(pln)
print (pln)
###Results###
#
# mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
# b[1] 2.557 0.281 2.002 2.557 3.108 FALSE 1 1 70919
# b[2] 2.057 0.081 1.898 2.056 2.215 FALSE 1 1 56997
# b[3] 0.975 0.076 0.826 0.975 1.125 FALSE 1 1 35692
# sigma 1.055 0.053 0.956 1.052 1.165 FALSE 1 1 58051
# deviance 2535.478 19.988 2498.335 2534.871 2576.864 FALSE 1 1 57735#
# pD = 199.8 and DIC = 2735.228