Извлечение прогнозов из байесовской модели случайных эффектов (JAGS) - PullRequest
Новая
       13

Извлечение прогнозов из байесовской модели случайных эффектов (JAGS)

1 голос
/ 07 марта 2019

Я новичок в байесовском учении и пытаюсь извлечь прогнозы и достоверные интервалы для построения графиков и не могу понять, как его кодировать.В идеале, я хотел бы, чтобы как общий прогноз (например, значения прогнозирования dist_feed, сравниваемый между двумя обработками), так и значения прогнозирования явно включали другие переменные (такие как сравнение самцов и самок в рамках каждой процедуры и гнездящихся самок и гнездящихся мужчин в каждой обработке).).Будем весьма благодарны за любые подсказки.

модель: 

sink("dist_feeder.jags") 
cat("
model {

for( r in 1 : Ndata ) {
dist_feeder[r] ~ dnorm(mu[r] , tau) 
mu[r] <-  alphanobo[nobo[r]] + b0[nobo[r]]*treatment[r] + b1[nobo[r]]*sex[r] + b2[nobo[r]]*nest[r] +b3[nobo[r]]*nest[r]*sex[r] + b4[nobo[r]]*brood[r]+ b5[nobo[r]]*brood[r]*sex[r]
}


for (i in 1:Nnobo) { alphanobo[i] ~ dnorm(0, taunobo)
                      b0[i] ~ dnorm(mutreatment,tautreatment)
                      b1[i] ~dnorm(mu1G , tau1G )
                      b2[i] ~dnorm(mu2G , tau2G)
                      b3[i] ~dnorm(mu3G , tau3G )
                      b4[i] ~dnorm(mu4G , tau4G )
                      b5[i] ~dnorm(mu5G , tau5G )
                      }

munobo ~ dnorm(0,0.01)

tau <- 1/(sigma*sigma)
sigma ~ dunif(0,500) 

taunobo <- 1/(sigma.nobo*sigma.nobo)

sigma.nobo ~ dunif(0,500)

mutreatment ~ dnorm(0,0.01) 
tautreatment ~ dgamma(0.1,0.1)

mu1G ~ dnorm(0,.01)
tau1G ~ dgamma(.1,.1)

mu2G ~ dnorm(0,.01)
tau2G ~ dgamma(.1,.1)

mu3G ~ dnorm(0,.01)
tau3G ~ dgamma(.1,.1)

mu4G ~ dnorm(0,.01)
tau4G ~ dgamma(.1,.1)

mu5G ~ dnorm(0,.01)
tau5G ~ dgamma(.1,.1)

 }

",fill = TRUE)
sink()


 ni <- 55000
 nt <- 1
 nb <- 25000
 nc <- 3
 parameters <- 

 c("munobo","sigma", "alphanobo","sigma.nobo","mutreatment", 
   "tautreatment","b0","b1", "b2","b3","b4","b5","mu1G","mu2G", 
   "mu3G","mu4G","mu5G","tau1G","tau2G","tau3G","tau4G","tau5G")

mod1 <- jagsUI(jags.data, inits=NULL, parameters, "dist_feeder.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)

попытка кода извлечения прогноза: 

pred_fdist <-matrix(nrow=80000,ncol=2)

for(i in 1:80000){
 pred_fdist[i,] <- mod1$sims.list$mutreatment[i] + mod1$sims.list$munobo[i] + mod1$sims.list$mu1G[i] + mod1$sims.list$mu2G[i] + mod1$sims.list$mu3G[i] + 
mod1$sims.list$mu4G[i] + mod1$sims.list$mu5G[i]
}
...