I'm using the package R2jags in R, which uses the parameters.to.save argument to specify parameters. I'm interested in the statistical distinction between a stochastic node (e.g., obs ~ dnorm(0, 0.01), where obs is not listed in parameters.to.save) vs a parameter (if obs were a parameter to save). Is this merely a programatic distinction, detailing which nodes should be returned to the user? Or is there a statistical distinction?

I've been told that there is no difference other than what is saved for user output. However, it seems plausible that there may be a statistical distinction, whereby parameters can have their distribution influenced by data, whereas the distribution of a stochastic node does not change with data.

I would experiment with this, but the only way I know how to track the value of a stochastic node in the model is to specify it as a parameter. You can see the conflict (although there is likely a more clever way to infer the effect that I haven't considered).

I realize this may ultimately be a programming issue, but I'm asking here because I'm really only concerned with the statistical interpretation/ distinction between a stochastic node and a parameter. I've read through the manual several times, and it's still a bit murky for me.


1 Answer 1


I tried experimenting, with the logic that I could compare the DIC of models to see if specifying a parameter mattered. Insofar as I can tell, it doesn't. In fact, the DIC doesn't change even if I don't specifying any parameters actually pertaining to observations (e.g., if I introduce a node in the model as rando~dnorm(0,1), then make "rando" the only parameter to save).

Here was the code I used for testing. It seems saving a node as a parameter only affects which values you can access afterwards, and does not affect the statistics.

x <- 0:29
beta1 <- 0.5
beta0 <- 1
y.true <- beta0 + beta1*x
y.obs <- y.true + rnorm(length(x), 0, 0.1)
y.obs[10] <- 0
y.obs.prec <- rep(1E3,length(x))
y.obs.prec[10] <- 1/1E3

data <- list(x=x, y.obs=y.obs, y.obs.prec=y.obs.prec, N=length(x))
params1 <- c("beta0","beta1","sigma2")
params2 <- c("beta0","beta1","sigma2","y.true[10]")

model <- function(){
    for(i in 1:N){
        y.obs[i] ~ dnorm(y.true[i], y.obs.prec[i])
        mu[i] <- beta0 + beta1*x[i]
        y.true[i] ~ dnorm(mu[i], 1/sigma2)
    beta0 ~ dnorm(0,0.001)
    beta1 ~ dnorm(0,0.001)
    sigma2 ~ dgamma(0.001, 0.001)


fit1 <- jags(data=data, parameters.to.save=params1, n.iter=5E3, model.file=model)
fit2 <- jags(data=data, parameters.to.save=params2, n.iter=5E3, model.file=model)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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