6
$\begingroup$

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.

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.
$\endgroup$
  • 1
    $\begingroup$ You could also look into Stan, which might be the next big thing for Bayesian MCMC. It generated c++ code for the models; as a consequence Stan works fine with ifs and elses. $\endgroup$ – Erik Feb 24 '14 at 11:05
4
$\begingroup$

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)
}
$\endgroup$
1
$\begingroup$

In this specific case you could set the beta2 term to zero

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

where indicator[] is a vector that is one for those data points that you want to model with beta2 and 0 otherwise. You could also use a scalar to change the model as a whole. However, this approach is probably less efficient because it still samples from an unwanted parameter distribution. I have used the gsub approach before, but I agree it's not pretty.

If the model is written and executed in an R script it can be manipulated as any other string:

modelstring <- paste("
model {  
  for ( i in 1:N ) {
",ifelse(timevariable,"
   mu[i] <- beta1 + beta2 * time[i]
","mu[i] <- beta1"),"
  }
  y[i]  ~ dnorm(mu[i], tau) 
}
")
writeLines(modelstring,con=filename)
$\endgroup$
  • $\begingroup$ +1 Thanks for that. Did you have a particular strategy for using gsub? E.g., do you use some form of ad hoc macro system? $\endgroup$ – Jeromy Anglim Feb 8 '14 at 3:41
  • 1
    $\begingroup$ I've amended my answer above. But this is getting into programming specifics. $\endgroup$ – user12719 Feb 8 '14 at 5:24
  • $\begingroup$ Thanks for that. That's a nice way of organising things. I agree there is always the dilemma of whether to post questions like this on stackoverflow or crossvalidated. Personally I think this would be relevant to both forums, and I'm glad you posted your code example. It's really helpful :-) $\endgroup$ – Jeromy Anglim Feb 8 '14 at 5:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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