Странный вывод пуассоновской GLM со случайным эффектом iid в r - PullRequest
Новая
       33

Странный вывод пуассоновской GLM со случайным эффектом iid в r

0 голосов
/ 20 октября 2018

Я пытаюсь запустить rjags в R (через Rstudio), чтобы оценить параметры alpha & beta и гиперпараметр tau.nu модели, следующей: 

y_i|x_i~pois(eta_i),
eta_i=exp(alpha + beta*x_i + nu_i),
nu_i~N(0,tau.nu)

есть мой код: 

#generating data
N = 1000
x = rnorm(N, mean=3,sd=1) 
nu = rnorm(N,0,0.01)
eta = exp(1 + 2*x + nu)
y = rpois(N,eta) 
data=data.frame(y=y,x=x)
###MCMC
library(rjags)
library(coda)
mod_string= "model {  
  for(i in 1:1000) {
    y[i]~dpois(eta[i])
    eta[i]=exp(alpha+beta*x[i]+nu[i])
    nu[i]~dnorm(0,tau.nu)
  }
  alpha  ~ dnorm(0,0.001)
  beta  ~ dnorm(0,0.001) 
  tau.nu ~ dgamma(0.01,0.01) 
}"

params = c("alpha","beta","tau.nu")

inits = function() {
  inits = list("alpha"=rnorm(1,0,100),"beta"=rnorm(1,0,80),"tau.nu"=rgamma(1,1,1))
}
mod = jags.model(textConnection(mod_string), data=data, inits=inits, n.chains =3)
update(mod,5000)
mod_sim = coda.samples(model=mod,
                       variable.names=params,
                       n.iter=2e4)
mod_csim = as.mcmc(do.call(rbind, mod_sim)) 
plot(mod_csim)

я получаю странный вывод, я не знаю, где я ошибаюсь. MCMC не работает в этой модели? Или я просто делаю что-то не так в кодировании?

enter image description here

1 Ответ

0 голосов
/ 22 октября 2018

Эта модель не сходится с использованием стандартных сэмплеров.Это происходит, если вы используете сэмплеры в модуле glm.(но это не всегда так [1] )

Без загруженного модуля glm 

library(rjags)

mod_sim1 <- jagsFUN(dat)
plot(mod_sim1)

enter image description hereПосле загрузки 

load.module("glm")
mod_sim2 <- jagsFUN(dat)
plot(mod_sim2)

enter image description here

# function and data
# generate data
set.seed(1)
N = 50 # reduced so could run example quickly
x = rnorm(N, mean=3,sd=1) 
nu = rnorm(N,0,0.01)
eta = exp(1 + 2*x + nu)
y = rpois(N,eta) 
dat = data.frame(y=y,x=x)

# jags model
jagsFUN <- function(data) {
  mod_string= "model {  
    for(i in 1:N) {
      y[i] ~ dpois(eta[i])
      log(eta[i]) = alpha + beta* x[i] + nu[i]
    }

    # moved prior outside the likelihood
    for(i in 1:N){
        nu[i] ~ dnorm(0,tau.nu)
    }
    alpha  ~ dnorm(0,0.001)
    beta  ~ dnorm(0,0.001) 
    tau.nu ~ dgamma(0.001,0.001) 
    # return on variance scale
    sig2 = 1 / tau.nu
  }"

  mod = jags.model(textConnection(mod_string), 
                   data=c(as.list(data),list(N=nrow(data))), 
                   n.chains = 3)
  update(mod,1000)
  mod_sim = coda.samples(model=mod,
                         variable.names=c("alpha","beta","sig2"),
                         n.iter=1e4)
  return(mod_sim)
}
...