Solved – How to conditionally run element of JAGS script based on user supplied variable

Background: I often find that when working with JAGS that I end up with a lot of JAGS scripts as I explore different models. After a while I might settle on a set of models that I'm going to report, but then there is the issue of keeping settings consistent across these models that are meant to be consistent. For example, the priors for parameters that are common should be consistent, or the variable names for parameters should be consistent. While consistency is not that difficult to achieve, it's often a little annoying. For example, I might decide that I want to call something theta rather than beta, and now all the scripts need updating. So I'm trying to work out how I can have one script that can serve multiple purposes.

Example:
Here's a simple example. Assume I have two models, one with time as a predictor and one without. I supply to jags a logical variable timevariable to indicate whether I want to model time or not. Thus, I'd like to be able to write something like this.

for (i in 1:length(y)) {     if (timevariable) {         mu[i] <- beta1 + beta2 * time[i]     } else {         mu[i] <- beta1     }     y[i]  ~ dnorm(mu[i], tau) }

However, if-then is not part of the lanugage. Note that I'm not trying to incorporate an if-then based on estimated values of parameters. Rather I just want to have a script that can take flexible input.

Contents hide

1 Question

2 Initial thoughts

2.1 Best Answer

3 Explanation of the system

4 Example

5 Demonstration

6 Similar Posts:

Question

What is a good way to conditionally run parts of a JAGS script based on a user supplied variable?

Initial thoughts

Put a macro name in the text and use a text replacement function in R like gsub to replace the macro name with the desired text. I guess this would work, but it seems like a bit of a hack. It would be nice to have everything inside the script.

Best Answer

Here is one example of implementing a basic macro substitution system for JAGS scripts.

Explanation of the system

Define a function that takes as arguments any optional elements of the script.
For any aspects of the script that vary across argument values, record a macro token. This should be some unique text. Starting and ending with some symbols may assist to make this unambiguous.
For each macro token include code for what the replacement string should be under all possible combinations of the function arguments.
Replace the macro tokens with appropriate placement text based on arguments.
The code below provides one example of how macros could be specified and how to apply the macros to the raw script (suggestions for doing this more elegantly are welcome).
The function returns a replaced model string that could if necessary be passed to a textConnection function for use with rjags.

I like this system for a few reasons:

the raw script is easy to read
the resulting script has proper indentation
you don't have to worry about errors related to commands appearing on the same line

Example

Specifically, this example, aims to allow the user to fit a particular type of multilevel nonlinear model. It is designed to allow for three functional forms: a two parameter power, three parameter power, and a three parameter exponential.

The macro section structures macros as a list of lists. The top level list contains one element for each macro token. For each macro token, there is the macro token text, and the conditional replacement text.

Finally, a for loop applies all the macro replacements to the raw script.

See below (scrolling is required):

jags_model <- function (f=c('power2', 'power3', 'exp3')) {     f <- match.arg(f)      # raw script     script <-  "model {     # Model     for (i in 1:length(y)) {         $FUNCTION         y[i]  ~ dnorm(mu[i], tau[subject[i]])     }      # Random coefficients     for (i in 1:N) {             theta1[i] ~ dnorm(theta1.mu, theta1.tau)T(0, 1000)         theta2[i] ~ dnorm(theta2.mu, theta2.tau)T(-10,  0)         $THETA3DISTRIBUTION         sigma[i] ~ dnorm(sigma.mu, sigma.tau)T(0, 100);         tau[i] <- 1/(sigma[i]^2)     }       theta1.mu  ~ dunif(0, 100)     theta2.mu   ~ dunif(-2, 0)     $THETA3PRIOR.MU     sigma.mu ~ dunif(0, 20)      theta1.sigma ~ dunif(0, 100)     theta2.sigma ~ dunif(0, 2)     $THETA3PRIOR.SIGMA     sigma.sigma ~ dunif(0, 10)      # Transformations     theta1.tau <- 1/(theta1.sigma^2)     theta2.tau <- 1/(theta2.sigma^2)     $THETA3.TAU      sigma.tau <- 1/(sigma.sigma^2) }"       # define macros     macros <- list(list("$FUNCTION",                 switch(f,                      power2="mu[i] <- theta1[subject[i]] * pow(trial[i], theta2[subject[i]])",                       power3="mu[i] <- theta1[subject[i]] * pow(trial[i], theta2[subject[i]]) + theta3[subject[i]];",                      exp3="mu[i] <- theta1[subject[i]] * exp(theta2[subject[i]] * (trial[i] - 1)) + theta3[subject[i]];") ),              list("$THETA3DISTRIBUTION",             switch(f,                   power3=, exp3= "theta3[i] ~ dnorm(theta3.mu, theta3.tau)T(0, 1000)",                   power2="") ),          list("$THETA3PRIOR.MU",                 switch(f,                       power3=, exp3= "theta3.mu  ~ dunif(0, 100)",                       power2="") ),              list("$THETA3PRIOR.SIGMA",             switch(f,                   power3=, exp3= "theta3.sigma ~ dunif(0, 100)",                   power2="") ),          list("$THETA3.TAU",             switch(f,                   power3=, exp3= "theta3.tau <- 1/(theta3.sigma^2)",                   power2="") )         )      # apply macros     for (m in seq(macros)) {         script <- gsub(macros[[m]][1], macros[[m]][2], script, fixed=TRUE)     }     script }

Demonstration

Thus, we can produce the processed JAGS model with

x <- jags_model(f='power3')

And if we want to view the resulting model, we can do

cat(x)

which results in

model {     # Model     for (i in 1:length(y)) {         mu[i] <- theta1[subject[i]] * pow(trial[i], theta2[subject[i]]) + theta3[subject[i]];         y[i]  ~ dnorm(mu[i], tau[subject[i]])     }      # Random coefficients     for (i in 1:N) {             theta1[i] ~ dnorm(theta1.mu, theta1.tau)T(0, 1000)         theta2[i] ~ dnorm(theta2.mu, theta2.tau)T(-10,  0)         theta3[i] ~ dnorm(theta3.mu, theta3.tau)T(0, 1000)         sigma[i] ~ dnorm(sigma.mu, sigma.tau)T(0, 100);         tau[i] <- 1/(sigma[i]^2)     }       theta1.mu  ~ dunif(0, 100)     theta2.mu   ~ dunif(-2, 0)     theta3.mu  ~ dunif(0, 100)     sigma.mu ~ dunif(0, 20)      theta1.sigma ~ dunif(0, 100)     theta2.sigma ~ dunif(0, 2)     theta3.sigma ~ dunif(0, 100)     sigma.sigma ~ dunif(0, 10)      # Transformations     theta1.tau <- 1/(theta1.sigma^2)     theta2.tau <- 1/(theta2.sigma^2)     theta3.tau <- 1/(theta3.sigma^2)      sigma.tau <- 1/(sigma.sigma^2) }