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?

enter image description here

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) — user20650, Oct 20 '18 at 15:22
Thx! www.I tried `glm`module 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. — Yuki, Oct 21 '18 at 15:06
`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. — user20650, Oct 21 '18 at 16:23
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](https://rdrr.io/github/mbtyers/jagsplot/man/jags.hist.html).of course,I did install `R2jags` — Yuki, Oct 25 '18 at 05:34
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 @user20650 -- its the @ that notifies users — user20650, Oct 25 '18 at 12:29
Thank you,you're so nice.And the ggplot did help though I still met some problem. — Yuki, Oct 27 '18 at 08:51
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) — user20650, Oct 27 '18 at 10:12
Hi,user20650.I work out with the histgram. I post my problem in[link](https://stackoverflow.com/questions/53020720/credibility-interval-with-respect-two-factors-using-ggplot2-in-r) .It would be nice if you take look at it. — Yuki, Oct 27 '18 at 10:16
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. — Yuki, Oct 27 '18 at 10:20
you're welcome. Your code is fine , although I needed to reduce the binwidth. — user20650, Oct 27 '18 at 10:26

score 6 · Accepted Answer · answered Oct 22 '18 at 19:02

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

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 case[1](https://rdrr.io/github/mbtyers/jagsplot/man/jags.hist.html),of course I did install the `R2jags` package. — Yuki, Oct 25 '18 at 05:18
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 — user20650, Oct 25 '18 at 09:44

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

1 Answers1

Linked