Monday, 6 August 2012

A bunch of R (and JAGS) scripts

I finally (nearly) got around to prepare the R code to replicate the examples in the book. I divided the examples by chapter and then linked to the R scripts and, for those involving Bayesian analysis, the associated JAGS models.

At the moment, the scripts basically cover 3 running examples (discussed in several parts of the book):
  1. MCMC.R. A Gibbs sampler for the very simple case of a semi-conjugated Normal model \begin{eqnarray*}y_i &\sim& \mbox{Normal}(\mu,\sigma^2)\\ \mu\mid\sigma^2 & \sim & \mbox{Normal}(\mu_0,\sigma^2_0) \\ \tau=1/\sigma^2 &\sim& \mbox{Gamma}(\alpha_0,\beta_0)\end{eqnarray*} to show the basics of the method. That's all written in R and I also provide some additional code to make some plots, showing convergence varying the number of iterations, something like this:

  2. normalModel.R. A linear regression model \begin{eqnarray*} y_i & \sim & \mbox{Normal}(\mu_i,\sigma^2) \\ \mu_i & = & \alpha + \beta x_i \\ \alpha,\beta & \sim & \mbox{Normal}(0,h^2) \\ \log(\sigma) & \sim & \mbox{Uniform}(-k,k) \end{eqnarray*} This has several slightly different specifications to show how things change when centring a covariate, or selecting a longer burn-in for the MCMC, or using thinning. There are actually two versions of this script and this modifies the original model to include blocking to improve convergence and estimate the predictive distribution. Both versions also include the JAGS code to run the different models.
  3. HEexample.R. This runs a (reasonably simple) health economic model whose aim is to estimate the cost-effectiveness of a new chemotherapy drug. This actually consists of two steps. The first one specifies the model in terms of a set of relevant parameters. These are estimated through MCMC (and the JAGS code is provided). Then, using BCEA, the actual cost-effectiveness analysis is run. The script provides the code to run the basic analysis, as well as more advanced ones (including the computation of the Expected Value of Partial Information, model average to account for uncertainty in model specification, decision-maker's risk aversion and the mixed analysis to consider non-optimal market configurations $-$ I've discussed this here).
Soon(-ish), I'll add the code for the examples in chapter 5, which are specific health economic evaluations. These are also a combination of R and JAGS codes that basically set up and run the Bayesian model via R2jags and then post-process the results to produce the health economic evaluation.


  1. So why back to JAGS and not R-INLA for those models ?

  2. The models we do in health economics are generally a bit too complex for JAGS. Most of the time they are decision models, including several related variables (often some are functions or random parameters). So, all in all JAGS really is the best bet in this case. If you take a look at the codes you can see some examples.