7
$\begingroup$

I found out that the function jags() in the R2jags package sometimes does not remove the burn in part even with the option n.burnin=##. Here is a very simple example in R (a simple linear model):

library(R2jags)

N <- 1000
y <- rnorm(N)
x <- rnorm(N)
data <- list("N", "y", "x")

inits <- function(){list(beta0=rnorm(1), beta1=rnorm(1), tau=1)}
parameters <- c("beta0", "beta1", "tau")

the model m.bug is like this:

model{
for (i in 1:N){
y[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 + beta1*x[i]
}

beta0 ~ dnorm(0, 0.00001)
beta1 ~ dnorm(0, 0.00001)
tau ~ dgamma(0.001, 0.001)
sigma2 <- 1/tau
}

using "jags()" in R2jags package like this:

m <- jags(data, inits, parameters, "m.bug", 
  n.chains=3, n.iter=2000, n.burnin=1000, n.thin=1)

My question:

In the output of m, the posterior estimates are based on the right number of interations (1000), but if we check the traceplot (using traceplot(m)), the burnin part seems still there(e.g. the first few values for "tau" are not converged). Why?

and there is only one progress bar (uaually two, one for burn-in, one for the rest).

also, if I change n.iter=2000, n.burnin=1000 to n.iter=2001000, n.burnin=2000000

the elapsed time does not change, which is "too fast" for so many iterations.

ps. I used R version 2.15.2 and R2jags version 0.03-08

$\endgroup$
  • $\begingroup$ I can not replicate that traceplot(m) shows the burn-in period. For me I only get the 1000 samples that I asked for. $\endgroup$ – Rasmus Bååth Dec 5 '12 at 16:49
  • $\begingroup$ if you check the traceplot of tau, you can see that the first few values are far from the posterior mean. It is not in your case? $\endgroup$ – Baoyue Li Dec 5 '12 at 22:54
  • $\begingroup$ No I don't get that. I also just get 1000 iteration in my plot, if the burn in period was also included. Wouldn't the plot show 2000 iterations? $\endgroup$ – Rasmus Bååth Dec 6 '12 at 10:45
  • $\begingroup$ Ok, I actually zoomed in on the plot now showing only 100 iterations and as you said it looks like it is showing the burn-in period. I don't know why, I can only recommend using rjags, the other R-jags bridge library. $\endgroup$ – Rasmus Bååth Dec 6 '12 at 11:13
  • 1
    $\begingroup$ If you don't find a solution to the problem maybe you should hand in a bug report to the maintainers R2jags? $\endgroup$ – Rasmus Bååth Dec 6 '12 at 13:44
5
$\begingroup$

I had exactly the same problem! After reading the R2jags::jags source I came to the same conclusion as @Hao Ye, I will expand on more detail:

First, I noticed that this happens only for some models:

  • if you remove the dgamma from your example and use dunif instead, i.e. if you modify your example this way:

    ...
    inits <- function(){list(beta0=rnorm(1), beta1=rnorm(1), sigma2=1)}
    parameters <- c("beta0", "beta1")
    ...
    sigma2 ~ dunif(0, 100)
    tau <- 1/sigma2
    ...
    

    then you will not observe the error

  • Error you observe is also produced by models with

    • combination of rules log(lambda) <- ... and dpois(lambda)
    • combination of rules log.lambda <- ... and dpois(exp(log.lambda)), but not for log.lambda <- ... and dpois(log.lambda), and not for log(lambda) <- ... and dnorm(lambda, tau)

Second, there is a bug in R2jags::jags (version 0.03-12) - as @Hao Ye noticed, it mistakenly takes burn-in parameter and uses it as adaptation, but for those models that don't need it (see above), the adaptation phase is skipped! See the code of R2jags::jags (version 0.03-12):

if (n.burnin > 0) {
    n.adapt <- n.burnin
}
else {
    n.adapt <- 100
}
....
m <- jags.model(model.file, data = data, inits = init.values, 
    n.chains = n.chains, n.adapt = 0)
adapt(m, n.iter = n.adapt, by = refresh, progress.bar = progress.bar, 
    end.adaptation = TRUE)
samples <- coda.samples(model = m, variable.names = parameters.to.save, 
    n.iter = (n.iter - n.burnin), thin = n.thin, by = refresh, 
    progress.bar = progress.bar)
....

The second call - call to adapt - will only do something in models that need the adaptation phase (I do not understand how this behaviour is defined in the ?adapt help - possibly it is some undefined default behaviour, or is it caused by the end.adaptation parameter? I don't know, the documentation seems insufficient. Anyway, it is like this).

Anyway, there is only one subsequent call to coda.samples which means that the real burn-in phase is missing. The adapt is run instead of burn-in, and only for some models.

I started to use package runjags instead of R2jags because of this bug.


Note: if someone decides to fix this bug it should be done in such a way that it does allow the JAGS' default length of the adaptation phase - see this problem in runjags package: https://stackoverflow.com/questions/22555421/runjags-how-to-use-the-jags-default-for-the-length-of-adaptation-phase

$\endgroup$
3
$\begingroup$

I think I'm a year too late to this party... (and I feel like this question might belong in stackoverflow, instead)

The issue seems to be that the jags() function treats n.burnin as the number of adaptation steps. Not all models require adaptation, which is why Baoyue has seen two progress bars for some models, but not for the example here.

I'm not sure if jags() has options for true burn-in (i.e. discard the initial portion of the MCMC run), so you might have to pad n.iter and manually discard. Alternatively, use rjags (written by the author of JAGS) instead.

$\endgroup$
  • $\begingroup$ +1 Very good answer, seems that you are completely right $\endgroup$ – Curious Mar 18 '14 at 18:08
2
$\begingroup$

Since I run a lot of analyses in JAGS and prefer the output format from R2jags I decided to hack together a temporary solution to this problem.

Essentially I added an argument to the jags() function in R2jags, n.adapt, to specify the number of 'adaptation' iterations and got rid of the code chunk that sets n.adapt = n.burnin. Then in between the R2jags calls to jags.model() [where the adaption occurs] and coda.samples() which generates samples from the posterior, I added a call to update() that updates the model for a number of iterations equal to the burnin period specified in the jags() call. Then, the output from coda.samples() will not contain those samples from the burnin period regardless of how JAGS decides to handle adaption.

My solution seems to work for the toy example provided in the question (i.e., the burn-in occurs and the traceplots look good)

You can download the modified R2jags package here. I've tested it on Windows and Linux (only the jags() function, though)

I'm no expert on building R packages, but maybe this will be helpful to others.

$\endgroup$
  • $\begingroup$ I've done a complete re-write of the package, hopefully fixing the issue described, and hosted it here. Apologies that it appeared to be offline for a while. $\endgroup$ – Ken Kellner Oct 2 '14 at 22:51
2
$\begingroup$

I've experienced this problem in my own models as well. A solution (i.e., work around until this bug is fixed): switch the prior from tau to sigma. Instead of

tau~dgamma(0.001,0.001)
sigma2<-1/tau

try

sigma2 ~ dgamma(0.001, 0.001)
tau<-1/sigma2

or

sigma2 ~ dgamma(0.001, 0.001)
tau <- 1/pow(sigma2,2)

This seems to fix the problem, and R2jags will NOT ignore the burn in period. Why this works, I don't know. Anyone see any problems with this approach? Obviously in my example you'll need to provide an initial value for sigma2, not tau (i.e., "sigma2 = 1").

$\endgroup$
0
$\begingroup$

Following up to Jason Hill's answer, while this fixes the problem, you shouldn't do it. If $\tau$, the precision, is the inverse of $\sigma^2$ (the variance), and $\tau$ is distributed gamma, then it follows that $\sigma^2$ will not be distributed gamma, but inverse gamma. I don't think JAGS supports that parameterization, which is why you see so many people using $\tau$ in their models.

Wikipedia has a summary of the distribution, if you are interested.

$\endgroup$
0
$\begingroup$

I encountered this same problem and tried to download Ken Kellner's fix. However, it appears his code is no longer online.

Based on his suggestions I was able to modify the R2jags::jags.R file to check if adaption is needed before running the adapt function. If it is not needed, then it simply calls update to burn in the model.

Here is the specific code to correctly burn in models that do not require adaption (modifying the code around line 157):

  # Handle situations where the model doesn't adapt
  # - Chris MacLellan
  if(.Call("is_adapting", m$ptr(), PACKAGE="rjags")){
    adapt( m,
           n.iter         = n.adapt,
           by             = refresh,
           progress.bar   = progress.bar,
           end.adaptation = TRUE )
  }
  else{
    update( m,
            n.iter        = n.burnin,
            by            = refresh,
            progress.bar  = progress.bar )
  }

If you want to just download the modified package, I have hosted it on my website here: http://christopia.net/data/R2jags_0.04-03.tar.gz

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.