How to dynamically create an openbugs model file (R, bash) given a set of model parameters

Question

I am fitting a set of models using r2openbugs. This requires me to create a new model file for each model I try. I have a set of combinations that I want to test (different covariate sets, whether to include random effects, whether to calculate predictions etc). I would like to create a script that takes a set of parameters (list of covariates to include, T/F for random effects and predictions) and outputs a valid openbugs model file.

Current attempt at the problem: I am currently writing out each model I wish to use and selecting between them using R closures. However, this is obviously duplicating a lot of code, and takes a long time. I'm sure it is possible to generate an appropriate model file by adding/ modifying text for a base model depending on some parameters, but I don't really know where to begin.

I would like to use R to do this, but would be happy to use bash if necessary.

Example models:

BASIC GLM MODEL

create_model(random_effects = F, 
             predictions = F, 
             cov_names = c(cov1, cov2, cov3))

would return

for (k in 1:N_data) {
  BetaX_data[k] <- beta[1] + 
                   beta[2] * cov1[k] + 
                   beta[3] * cov2[k] + 
                   beta[4] * cov3[k]

  logit(pi_data[k]) <- BetaX_data[k]
  m[k] ~ dbin(pi_data[k], n[k])
}
## PRIORS ##
for (g in 1:N_cov) {
  beta[g] ~ dnorm(0, 0.00001) 
}

GLMM MODEL (same as previous model, but add lines for random effect)

create_model(random_effects = T, 
             predictions = F, 
             cov_names = c(cov1, cov2, cov3))

would return

for (i in 1:N_loc) { ## additional lines
  U[i] ~ dnorm(0, tau_iid)
}
for (k in 1:N_data) {
  BetaX_data[k] <- beta[1] + 
                   beta[2] * cov1[k] + 
                   beta[3] * cov2[k] + 
                   beta[4] * cov3[k]

  logit(pi_data[k]) <- BetaX_data[k] + U[location_id[k]]
  m[k] ~ dbin(pi_data[k], n[k])
}
## PRIORS ##
for (g in 1:N_cov) {
  beta[g] ~ dnorm(0, 0.00001) 
}
tau_iid ~ dgamma(0.001,0.001) ## additional lines

GLM MODEL (with cov2 removed)

create_model(random_effects = F, 
             predictions = F, 
             cov_names = c(cov1, cov3))

would return

for (k in 1:N_data) {
  BetaX_data[k] <- beta[1] + 
                   beta[2] * cov1[k] + 
                   beta[3] * cov3[k]  ## lines changed

  logit(pi_data[k]) <- BetaX_data[k]
  m[k] ~ dbin(pi_data[k], n[k])
}
## PRIORS ##
for (g in 1:N_cov) {
  beta[g] ~ dnorm(0, 0.00001) 
}

GLM MODEL (with predictions)

create_model(random_effects = F, 
             predictions = T, 
             cov_names = c(cov1, cov2, cov3))

would return

for (k in 1:N_data) {
  BetaX_data[k] <- beta[1] + 
                   beta[2] * cov1[k] + 
                   beta[3] * cov2[k] + 
                   beta[4] * cov3[k]

  logit(pi_data[k]) <- BetaX_data[k]
  m[k] ~ dbin(pi_data[k], n[k])
}
## PRIORS ##
for (g in 1:N_cov) {
  beta[g] ~ dnorm(0, 0.00001) 
}
## PREDICTIONS ##
for(l in 1:N_pred) { # additional lines which also 
                     # change depending on covariate set
  BetaX_pred[l] <- beta[1] + 
                   beta[2] * cov1[l] + 
                   beta[3] * cov2[l] + 
                   beta[4] * cov3[l]

  logit(pi_pred[l]) <- BetaX_pred[l]
}

Answer 1

I don't know of a way to do this specifically for OpenBUGS, but you could check out the template.jags function in the runjags package (https://www.rdocumentation.org/packages/runjags/versions/2.0.4-2/topics/template.jags) which does pretty exactly what you ask for JAGS models but based on an lme4-style model syntax. See also the 'Generating a GLMM template' section from here:

http://runjags.sourceforge.net/quickjags.html

If you used write.data=FALSE and write.inits=FALSE then it just creates the model, which should also be valid OpenBUGS code. Or you could just install JAGS and use that instead :)

Hope that helps.

Matt