# Effect of “parameters.to.save” in R2jags/ JAGS

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.

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 <- 0
y.obs.prec <- rep(1E3,length(x))
y.obs.prec <- 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")

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)