Weird output of poisson GLM with an iid random effect in r

Question

I am trying run rjags in R (via Rstudio) to estimate parameters alpha&beta and hyperparameter tau.nu of the model following:

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

there is my code:

#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)

the I get weird output,I don't konw where I get wrong.Does MCMC not work in this model?Or I just do something wrong in coding?

Fluoroscopy

I made a couple of changes: moved the prior for nu outside the likelihood, increased burnin and iterations -- these made no difference. After loading the glm module (load.module("glm")) it converged as expected. (caveat, I used N=100 as it would take too long otherwise) – Cotemporary

Thx! www.I tried glmmodule that you suggested,it converged,but i still have problems:I set iterations as 5e4,why the range of x-axis is 0~1.5e5 ?And the value of estimation of tau.mu is far away from true value 0.01. – Fluoroscopy

tau.nu is the precision. to return the variance add a line under the priors, sig2 = 1 / tau.nu, Then follow sig2 instead of tau.nu. From your example, the iterations on the x-axis depend the n.adapt (1000) from jags.model + 5000 burnin + n.iter (2e4) from coda.samples, so I'd expect the x-axis to start at 6000, and go to 26000 in plot(mod_sim) - which is what I think you really want. As you are combining the chains you get 3*5e4 = 1.5e5 -- so one after the other. – Cotemporary

Hi,www.I don't know how to continue asking my problem about this case,so I try to comment who answered me.Hope you don't mind.I still wondering whether there is a way to plot interval estimates and hist from this simulation mcmc draws ? I tried bayesplot package ,it made problems when I try to add output of other approach.So I tried jags.hist function in R2jags package, it errors with 'not found jags.hist function',even I run the example case 1.of course,I did install R2jags – Fluoroscopy

Yuki, it was me not www who answered you in the comments (you can see who left the comments by username given at the end of the comment). ps to notify me you can use @Cotemporary -- its the @ that notifies users – Cotemporary

Thank you,you're so nice.And the ggplot did help though I still met some problem. – Fluoroscopy

here is some code: library(data.table) ; chn <- rbindlist(lapply(mod_sim2, as.data.table), idcol="chain") ; mchn <- melt(chn, id.vars = "chain") ; library(ggplot2) ; ggplot(mchn, aes(value)) + geom_histogram() + facet_wrap(~variable, scales="free") (ps density plots are definitely the better option compared to histograms) – Cotemporary

Hi,user20650.I work out with the histgram. I post my problem inlink .It would be nice if you take look at it. – Fluoroscopy

About hist thing.rjags.mat = as.matrix(mod_csim);rjags.dat = as.data.frame(rjags.mat) ;hist.sig2=ggplot() + geom_histogram(data=rjags.dat,aes(x=sig2,y=..density..), binwidth=.1, colour="gray", fill="white") + geom_density(data = dens.inla,aes(x=sig2),col='black') I tired like this.And thank you sooooooo much. – Fluoroscopy

you're welcome. Your code is fine , although I needed to reduce the binwidth. – Cotemporary

Answer 1

This model doesn't converge using the standard samplers. It does if you use the the samplers in the glm module. (but this may not always be the case [1])

Without the glm module loaded

library(rjags)

mod_sim1 <- jagsFUN(dat)
plot(mod_sim1)

After loading

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

---

# 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)
}

Cotemporary

Hi,I still wondering whether there is a way to plot interval estimates and hist from this simulation mcmc draws ? I tried bayesplot package ,it made problems when I try to add output of other approach.So I tried jags.hist function in R2jags package, it errors with 'not found jags.hist function',even I run the example case1,of course I did install the R2jags package. – Fluoroscopy

Hi, I dont know what plot intervals you want. Re the histogram, you can just pull out the relevant columns from the mcmc output and plot i.e. hist(mod_sim2[[1]][,1]). I suppose you could reshape the mcmc output and plot all parameters at once using ggplot – Cotemporary