Simple linear regression in WinBUGS

---

We observed the data:
X:   1   2   3   4   5   6   7   8   9 10
Y:   6   7   8   4   9 11 12 14 15 19

Here is WinBUGS code to perform simple linear regression

#The model:
model{
   for(i in 1:n){
      y[i]~dnorm(mu[i],taue) #taue is a precision, not a variance!
      mu[i] <- alpha + beta*x[i]
   }
   taue~dgamma(0.01,0.01)
   alpha~dnorm(0,0.01)
   beta~dnorm(0,0.01)
}

#The data
list(n=10,x=c(1,2,3,4,5,6,7,8,9,10),y=c(6,7,8,4,9,11,12,14,15,19))

#The initial values
list(taue=1,alpha=0,beta=1)



Running WinBUGS

• Highlight the word `model' in the first line of the model code.
  
• Go to `model/specification' and click `check model'. `Model is syntactically correct' should pop up on the bottom of the WinBUGS window.
  
• Highlight `list' in the data list and click `load data'. `data loaded' should pop up on the bottom of the WinBUGS window.
  
• Highlight `list' in the initial values list and click `load inits'. Click `gen inits' to have BUGS randomly generate initial values.
  
• For each variable you'd like to analyze, go to `inference/sample', enter the name of the variable in the `nodes' box, and click `set'. Samples of variables you don't request to record are automatically discarded.
  
• To begin sampling, go to model/update and click `update'.
  
• When the MCMC sampling is completed, go back to `inference/sample' to view summaries of the posterior.

Calling WinBUGS from R

---

Step 1) Create a file with the WinBUGS commands to specify the model. For example, create the file ``C:/Documents and Settings/bjreich/Desktop/BUGSmodel.txt'' with the following lines:

model{
   for(i in 1:n){
     y[i]~dnorm(mu[i],taue) #taue is a precision, not a variance!
     mu[i] <- alpha + beta*x[i]
   }
   taue~dgamma(0.01,0.01)
   alpha~dnorm(0,0.01)
   beta~dnorm(0,0.01)
}

Step 2) Load the data in R
x<-c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
y<-c(6, 7, 8, 4, 9, 11, 12, 14, 15, 19)
n<-10

Step 3) Call WinBUGS from R

First, load the package ``R2WinBUGS'' from CRAN. Then submit the following code in R.

library(``R2WinBUGS'')
modelfile<-``C:/Documents and Settings/bjreich/Desktop/BUGSmodel.txt''
data<-list(``n'',``x'',``y'')
vars2keep<-list(``taue'',``alpha'',``beta'')
inits<-function(){list(taue=1, alpha=0, beta=0)}

output<-bugs(
   model.file=file.path(modelfile),
   data=data,
   inits = inits,
   parameters.to.save = vars2keep,
   n.chains=3,
   n.iter=50000,
   n.burnin=5000,
   n.thin=10,
   debug=TRUE,
   DIC=TRUE
)

Step 4) Analyze the MCMC output in R

# Check convergence
   plot(output)

# Plot the posterior of beta
   samples<-output$sims.matrix #samples is a matrix with all the MCMC samples. Each variable is a column
   hist(samples[,3],main=``Posterior of beta'') #the third column is beta