Can we enforce that the user can only set scalar parameters?

[joethorley]

If we can then it will be much easier to set parameters in a jags model.

Note we shouldn't do it by replacing say alpha = 0.5 everywhere because

 alpha ~ dunif(0, 1)

will produce the following invalid code

0.5 ~ dun(0, 1)

instead we should search for alpha ~ ... and replace the ~ and everything after with 0.5.

[joethorley]

A second benefit of enforcing scalars is that the parameter set can be given by a uniquely named numeric vector ie

c(alpha = 1, beta = 3.1)

also multiple combinations could be set by list(alpha = c(1,2), beta = 3.1) which would automatically be processed to produce

list(c(alpha = 1, beta = 3.1), c(alpha = 2, beta = 3.1))

[audrey-b]

If we can then it will be much easier to set parameters in a jags model. Note we shouldn't do it by replacing say alpha = 0.5 everywhere because alpha ~ dunif(0, 1)

will produce the following invalid code 0.5 ~ dun(0, 1)

instead we should search for alpha ~ ... and replace the ~ and everything after with 0.5.

Or we could simply erase the whole alpha ~ … line and pass alpha = 0.5 as data?

[audrey-b]

What would someone do if they need to pass values for alpha[i,j] for all combinations of i and j? E.g. in capture-recapture with 30 years of data, need to pass phi[1] to phi[30], etc.

[audrey-b]

Joe's idea (need to think about) : could we pass the parameters as inits in a model block? I think that what we have to do to test this is run the chain for 1 iteration and see if the parameter values have changed in the results. If they have changed, it means that jags has reported the next iteration and this idea won't work. If they haven't changed then maybe it's doing what we want?

[audrey-b]

Joe's idea (need to think about) : could we pass the parameters as inits in a model block? I think that what we have to do to test this is run the chain for 1 iteration and see if the parameter values have changed in the results. If they have changed, it means that jags has reported the next iteration and this idea won't work. If they haven't changed then maybe it's doing what we want?

I tried the code below and the first iteration reported for the mcmc algorithm has updated parameter values, so I don't think we can use a model block to generate data:

library(R2jags) library(runjags) library(mcmcplots)

likelihood <- "

Likelihood:

for (i in 1:N){ y[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance) mu[i] <- alpha + beta * x[i] }"

parameters for simulations

N = 30 # nb of observations x <- 1:N # predictor alpha = 0.5 # intercept beta = 1 # slope sigma <- .1 # residual sd tau <- 1/(sigma*sigma) # precision

priors <- "

Priors:

alpha ~ dnorm(0, 0.01) # intercept beta ~ dnorm(0, 0.01) # slope sigma ~ dunif(0, 100) # standard deviation tau <- 1 / (sigma * sigma) "

modell <- " model { %s %s } "

specify model in BUGS language

model <- paste(sprintf(modell, likelihood, priors))

data

jags.data <- list(y = rep(NA, N), N = N, x = x)

initial values

inits <- function(){list(alpha = alpha, beta = beta, sigma = sigma)}

parameters monitored

parameters <- c("alpha", "beta", "sigma")

MCMC settings

ni <- 2 nt <- 1 nb <- 0 nc <- 1

call JAGS from R

res <- jags(jags.data, inits, parameters, model=textConnection(model), n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, working.directory = getwd(), DIC=FALSE)

as.mcmc(res)

try runjags instead

res2 <- run.jags(data=jags.data, inits=inits, monitor=parameters, model=model, n.chains = nc, thin = nt, sample = ni, burnin = nb, adapt=0, method="simple")

names(res2) res2$mcmc

[joethorley]

OK I'll check it out!

[joethorley]

@audrey-b - I concur - so we definitely have to use a data block.

[joethorley]

I like you proposition that if we have alpha = 0.5 in the data then we simply delete

alpha ~ dunif(0, 1) in the model block.

as you suggest this should also work with vectors or arrays

beta <- c(1, 3)

for(i in 1:2) { beta[i] ~ dunif(0,2) }

After deleting all such nodes we could then run a function to delete empty for loops.

With this approach I'm thinking there is no need to move anything to a data block.

[audrey-b]

Sounds good to me. So the user will specify a model using model{}. Should we print a message to say that the prior distribution on beta was not used to simulate the data? There is a risk that a user forgets to specify the value for a parameter and gets incorrect results because we would be sampling from the prior...

[audrey-b]

What if the user is using a hierarchical prior e.g. for(i in 1:2) { beta[i] ~ dnorm(mu,2) } mu ~ dunif(-3,3)

but then decides to input values for beta, but not mu? The beta line would be erased but mu ~ dunif(-3,3) would remain in the code. I guess that's not a big deal as they wouldn't care about mu in this case.

[audrey-b]

Having the likelihood and prior separate as I had originally insures the user thinks through the model they want to use to generate the data. But the drawback is that if they already have a bugs model defined, they have to go through it to delete the lines of code that contain priors.

[joethorley]

we could have an argument called something like uncertain_priors = FALSE (that takes values of FALSE (errors), TRUE (no problem) and NA (warning)) if there are any priors without fixed values. By default it would be FALSE but the user could relax it if they want to have uncertainty in a particular parameter for a set of simulations? also it seems to me that this approach would still allows the user to edit a bugs model and use the data block if they wished.

[audrey-b]

It would be hard to identify priors without fixed values because of the case: lb <- 0 ub <- 1 beta ~ dunif(lb,ub) Here beta is a prior but it's not obvious because the numeric values are specified through variables. I think that it would be nice to have a chat about this as I am leaning more towards specifying the likelihood and priors separately but would like to understand your perspective better.

[joethorley]

allow user to pass matrices and arrays