Simple linear regression using R2OpenBUGS

First, let's load the package:
library(R2OpenBUGS)Next, write the model as though it was an R function. Note that the function does not have any arguments. Syntax used here is BUGS syntax. The function will be translated into BUGS code, so it's important to retain appropriate syntax for BUGS. If using constraints on parameters e.g. “I(0,)” then the syntax differs from the standard BUGS syntax - see the later example (Alternative 4). Note that in BUGS the ordering of likelihood and priors does not matter. The fact that we use alpha and beta before defining the prior does not matter. The convention is to specify the likelihood first and priors second.

linemodel <- function() {   
for (j in 1:N) {        Y[j ~ dnorm(mu[j, tau)  ## Response values Y are Normally distributed        
mu[j <- alpha + beta * (x[j - xbar)  ## linear model with x values centred   
}   

## Priors   
alpha ~ dnorm(0, 0.001)   
beta ~ dnorm(0, 0.001)   
tau ~ dgamma(0.001, 0.001)   
sigma <- 1/sqrt(tau)}

Next, let's prepare the data for analysis. BUGS accepts named list format data, including data structures with named columns (or mixtures of these). Note that BUGS doesn't react well if data passed to it is not used in the model… The best advice is to only pass what you need. Fortunately R is the perfect tool for preparing data for BUGS.

linedata <- list(Y = c(1, 3, 3, 3, 5), x = c(1, 2, 3, 4, 5), N = 5, xbar = 3)

The next step is to provide initial values for the MCMC. You should provide initial values for all random variables in the model. If you do not, BUGS will attempt to derive initial values for any random variables that are not initialised by using what you have provided. However sometimes it is unable to do so. The initial values are passed as a function with no arguments containing a named list with names corresponding to each random variable in the MCMC. R2OpenBUGS will call this function for every MCMC chain that is specified, generating initial values for each chain. Here we have integer values for alpha and beta. As we are running only 1 MCMC chain this is acceptable, but if we were running >1 chain, each chain would have the same starting values.

lineinits <- function() {    list(alpha = 1, beta = 1, tau = 1)}

Finally, let's run BUGS using the bugs(…) command and assign the output to an object called “lineout”. Note that we need to specify which nodes/parameters (random variables) in the model to save at each step of the MCMC. Basically, anything except data can be monitored by BUGS - model parameters, missing data values, functions of model parameters… This is a strength of BUGS since all of the parameters will have their posterior distribution samples saved for summarising, graphing etc. It can also place an overhead in data handling however: saving 1000 MCMC samples x k chains for elements such as fitted values e.g. mu[j] above can lead to a lot of data values. In this case it is trivial, but for other, larger datasets it can lead to problems.

Note also that here we only specify the number of MCMC iterations to be performed. By default R2OpenBUGS assumes that half of these iterations will be used as the burn-in. So in this instance the posterior distribution will be based on a sample of 500 MCMC iterations. In the case of this linear model that is probably sufficient - but also the number of iterations used as burn-in here is probably overkill.

## Uses default settings for n.burnin = n.iter/2; n.thin=1;
lineout <- bugs(data = linedata, inits = lineinits, parameters.to.save = c("alpha",     "beta", "sigma"), model.file = linemodel, n.chains = 1, n.iter = 10000)

Testing Bugs Installation with Examples

The purpose of this annotated R code is simply to check that your installation of OpenBUGS, JAGS, BRugs, R2OpenBUGS, Rjags & R2jags are all running and talking to each other in the way that they should. It is assumed that all of these programs have been installed (see, e.g., "Installing BUGS and the R to BUGS Interface"; the installing.bugs.jags.pdf is available on the Psych 548 website). The R-code for this document is shown as text at the end of this file, and also as an ascii file, test.bugs.install.txt. If you want to run the code as you are reading this document, I recomment that you load the text file, test.bugs.install.txt, into RStudio or your favorite programming editor. It will be easier to transfer the code to R from a programming editor than from this pdf file.

Two Way Anova using BUGS

Main Contents

• If you want to use the BRugs package and OpenBUGS for this analysis, ....
• Define the settings for either the OpenBUGS or JAGS runs
• Use BRugs to run OpenBUGS on the data and model file in ANOVAonewayBRugs.R
• Use rjags to run JAGS on the data and model file in ANOVAonewayBRugs.R
• Computing contrasts among group means
• Why Bayesian analysis does not introduce a correction for the number of contrasts that are tested?
• Why does the model file use the Student's t distribution to sample the standard deviation of the hyperprior of the group means? Why not use a gamma distribution?
  

Hierarchical Model using WinBUGS and R(Bayesian Data Analysis,Andrew Gelman)
http://www.stat.columbia.edu/~gelman/bugsR/software.pdf

library(R2OpenBUGS)
############################################################
###########             Load data               ############
############################################################
dat<-read.csv("F:\\ST733\\Data\\Obama2012.csv")
W<-as.matrix(read.csv("F:\\ST733\\Data\\NCADJ.csv",row.names=1))

Y<-dat[,2] #Percent that voted for Obama
X<-dat[,8] #Percent black


############################################################
###########           Set up for BUGS           ############
############################################################
n<-ncol(W)
num<-colSums(W)
adj<-NULL
for(j in 1:n){
  adj<-c(adj,which(W[j,]==1))
}
adj<-as.vector(adj)
num<-as.vector(num)
weights<-1+0*adj

# num[i] = number of regions that adjacent to region i
# The first num[1] elements of adj are region 1's neighbors,
# the next num[2] elements of adj are region 2's neighbors, etc.
# weights being all one specify that all neighbors are weighted equally


############################################################
###     CAR model with spatially-varying coefficients    ###
############################################################

sink("svc_car_model.txt")
cat("

model{               
   for(i in 1:n){  #Likelihood
       Y[i]  ~  dnorm(mu[i],tau[1])
       mu[i] <- a[i] + b[i]*X[i]
       a[i]  <- int[1] + R1[i]  #Add in an intercept to get the unrestricted CAR.
       b[i]  <- c*a[i] + (int[2] + R2[i])
    }

   #CAR model (in BUGS, the R[j,] are forced to sum to zero)
   R1[1:n] ~ car.normal(adj[], weights[], num[], tau[2])
   R2[1:n] ~ car.normal(adj[], weights[], num[], tau[3])

   #Priors
   int[1]~dflat()    # Must be flat to get the ICAR
   int[2]~dflat()   
   c~dnorm(0,0.1)   
   tau[1] ~ dgamma(0.1,0.1)
   tau[2] ~ dgamma(0.1,0.1)
   tau[3] ~ dgamma(0.1,0.1)
}

}", fill = TRUE)
sink()



############################################################
#########              Call BUGS                 ###########
############################################################


data<-list(X=X,Y=Y,n=n,num=num,adj=adj,weights=weights)
inits<-function(){
    list(int=c(.25,1),c=0,tau=rep(10,3))
}
keepers <- c("a","b","c")

out <- bugs(
        data=data,
        inits=inits,
        parameters.to.save=keepers,
        model.file="svc_car_model.txt",
        n.iter=10000,
        n.chains=3,
        n.burnin=1000,
        n.thin=1,   
        debug=FALSE,
#Change this to the location of the openBUGS executable on your computer
        OpenBUGS.pgm="F:\\OpenBUGS322\\OpenBUGS.exe",
        )



############################################################
#########       Summarize the results            ###########
############################################################

plot.nc<-function(data,main="",breaks=5,sig_digits=3){
   library(maps)
   d<-data[c(1:26,27,27,27:100)] #county 27 has three polygons!
   if(length(breaks)==1){
     breaks<-quantile(d,seq(0,1,length=breaks+1),na.rm=T)
   }
   nbreaks<-length(breaks)-1
   low<-breaks[-(nbreaks+1)]
   high<-breaks[-1]
   col<-rep(1,nbreaks)
   for(j in 2:nbreaks){
     col<-ifelse(d>low[j],j,col)
   }  
   shade<-1-col/nbreaks
   shades<-1-1:nbreaks/nbreaks

   low<-round(low,sig_digits)
   high<-round(high,sig_digits)
   map("county","north carolina",fill=T,col=gray(shade))
   title(main)
   legend("bottomleft",paste(low,"-",high),lwd=5,ncol=2,col=gray(shades))
}


plot.nc(X,main="Pct black")
plot.nc(Y,main="Pct Obama")
plot.nc(out$mean$a,main="Posterior mean a")
plot.nc(out$mean$b,main="Posterior mean b")

b<-out$sims.list$b
b_mean<-colMeans(b)
boxplot(b[,order(b_mean)],outline=F)

#The model:
model{
   for(i in 1:n){
       y[i]~dnorm(mu[i],tau)
       mu[i] <- alpha + beta*x[i]
   }
   tau<-1/sigma*sigma
   sigma~dunif(0,20)            
   alpha~dnorm(0,0.01)
   beta~dnorm(0,0.01)
}
===============================================================
#Load the data
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)

#Call WinBUGS:

library(R2OpenBUGS)
modelfile<-"F:\\ST740\\BUGSmodel.txt"
dat<-list("n","x","y")
vars2keep<-list("sigma","alpha","beta")
inits<-function(){list(sigma=1, alpha=rnorm(1))}

output<-bugs(
   model.file=file.path(modelfile),
   data=dat,
   inits = inits,
   parameters.to.save = vars2keep,
   n.chains=5,
   n.iter=1000,
   n.burnin=100,
   n.thin=1,
   debug=T,
   DIC=F,
   OpenBUGS.pgm="F:\\OpenBUGS322\\OpenBUGS.exe"
)

#Analyze the MCMC output

# Check convergence
  plot(output)   

# Plot the posterior of sigma
  samples<-output$sims.matrix #samples is a matrix with all the MCMC samples.
                              #each variable is a column
  hist(samples[,1],main="Posterior of sigma")    #the first column is sigma

#Print the posterior means and standard deviations
  print(output$mean)
  print(output$sd)

A Bayesian Nobel recipient

Chris Sims, who got the Nobel Prize in Economics, or more precisely the Sveriges Riksbank Prize in Economic Sciences, has [also] done work in Bayesian econometrics. See for instance his talk on Why Econometrics Should Always and Everywhere Be Bayesian. Or his analysis of a counterexample of Larry Wasserman’s. He even has a tech report on adaptive Metropolis-Hastings algorithms (that apparently did not get published). He has also been teaching Bayesian statistics and econometrics at Princeton for many years, so this is a cool day for Bayesian stats! (The picture of a switching regime estimation on his webpage is actually similar to a rendering of mine from the late 90′s)

我来标记一下