Linear regression in OpenBUGS

I always wondered why is it so difficult to find an OpenBUGS example of simple linear regression on the Web. Curiously, such example is even missing in the OpenBUGS help. The only nice example so far is in the book by Marc Kéry. I have simplified the code. You need to have OpenBUGS (or WinBUGS) installed, together with the R2OpenBUGS package. Then just copy-paste the code into R.
Good luck!

# generation of artificial data
x <- sort(rnorm(100))
y <- 5 + 3*x + rnorm(100)
plot(x,y)

# putting the data in the OpenBugs-friendly format
my.data <- list(x=x, y=y, N=100)

# loading the R2OpenBUGS package
library(R2OpenBUGS)

# sinking the model into a file
sink("regression.txt")
cat(" 
model
{
  # priors
  beta0 ~ dnorm(0,0.001)
  beta1 ~ dnorm(0,0.001)
  sigma ~ dunif(0,100)

  tau <- 1/(sigma*sigma)

  # likelihood
  for(i in 1:N)
  {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta0 + beta1*x[i]
  }
}
")
sink()

# setting the parameters to be monitored
params <- c("beta0", "beta1")

# specifying the initial MCMC values
inits <- function()
{
  list (beta0=rnorm(1,0,10000),
        beta1=rnorm(1,0,10000))
}

# calling OpenBugs to sample from the posterior distributions
# (you may need to change the OpenBUGS.pgm directory)
res <- bugs(data=my.data, inits=inits, parameters.to.save=params,
            n.iter=2000, n.chains=3, n.burnin=1000,
            model.file="regression.txt", debug=FALSE, codaPkg=TRUE,
            OpenBUGS.pgm="C:/Program Files (x86)/OpenBUGS/OpenBUGS321/OpenBUGS.exe",
            working.directory=getwd())

# and, finally, getting the results
br <- read.bugs(res)
plot(br, ask=T)
summary(br)

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>