RUNJAGS-set seed without prior for simulations

I am sampling some data from a normal distribution using runjags. I don't have any prior for the parameters that I used for my simulations. It seems that runjages dos not use the argument to fix the seed: list(".RNG.name"="base::Super-Duper", ".RNG.seed"=1). I changed the argument to list(muOfClustsim=rep(1, npop), ".RNG.name"="base::Super-Duper", ".RNG.seed"=1) but it does not work either. Is-there a way to fix the seed for such models in runjags?

Here is a minimal reproducible example:

library(runjags)

npop=3
nrep=10
sdpop=7
sigma=5
seed=4

set.seed(seed)
N = npop*nrep # nb of observations

## Population identity of each individual used to sample genotypes but not used for common garden test
pop <- rep(1:npop, each=nrep)

muOfClustsim <- rnorm(npop, 0, sdpop) # vector of population means
(tausim <- 1/(sigma*sigma)) # precision of random individual error

# parameters are treated as data for the simulation step
data <- list(N=N, pop=pop, muOfClustsim=muOfClustsim, tausim=tausim)

## JAG model
txtstring <- "
data{
  # Likelihood:
  for (i in 1:N){
    ysim[i] ~ dnorm(eta[i], tausim) # tau is precision (1 / variance)
    eta[i] <- muOfClustsim[pop[i]]
  }

}

model{
fake <- 0
}
"
## Initial values with seed for reproducibility
initssim <- list(".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
##initssim <- list(muOfClustsim=rep(1, npop), ".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
## Simulate with jags
set.seed(seed)
out <- run.jags(txtstring, data = data, monitor=c("ysim"), sample=1, n.chains=1, inits=initssim, summarise=FALSE)
## reformat the outputs
(ysim1 <- coda::as.mcmc(out)[1:N])

set.seed(seed)
out <- run.jags(txtstring, data = data, monitor=c("ysim"), sample=1, n.chains=1, inits=initssim, summarise=FALSE)
## reformat the outputs
(ysim2 <- coda::as.mcmc(out)[1:N])

identical(ysim1, ysim2)

Solution

The .RNG.name / .RNG.seed / .RNG.state initial values apply to the model (or more specifically, a chain within the model), and are not used within the data block. This means that (as far as I know) there is no way to make any stochastic representations within the data block reproducible in JAGS <= 4.3. This is something that could be added for a future version of JAGS, but it would be a low priority I'm afraid as it is always possible (and usually better) to simulate data in R before passing it to JAGS.

In your case the answer (assuming you want to use JAGS) is to do the simulation in the model block rather than the data block:

txtstring <- "
model{
  # Likelihood:
  for (i in 1:N){
    ysim[i] ~ dnorm(eta[i], tausim) # tau is precision (1 / variance)
    eta[i] <- muOfClustsim[pop[i]]
  }
}
"

The rest of your code then runs as expected §. It is worth noting that data simulation is a task generally better suited to R than JAGS, but I am assuming that there is a specific reason that you want to use JAGS in this case...

Matt

§ Although you should not generally expect strict equality of doubles, for example:

    identical(0.1+0.2, 0.3)

But:

    abs(0.3 - (0.1+0.2)) < sqrt(.Machine$double.eps)

Or even better:

    isTRUE(all.equal(0.1+0.2, 0.3))

This is well worth watching: https://www.youtube.com/watch?v=3Bu7QUxzIbA&t=1s