Example of running BUGS for Bayesian inference from R

ROpenBUGS_example.R

####################################

# TITLE: ONEWAY.R

# PURPOSE: ILLUSTRATE ONE WAY ANALYSIS OF VARIANCE

# DATA: AGES AT WHICH INFANTS FIRST WALKED ALONE, FROM FISHER AND VAN BELLE (2nd edition),

# PAGE 359. FOUR GROUPS: ACTIVE GROUP,

#       PASSIVE GROUP, NO-EXCERSISE GROUP AND CONTROL GROUP.

# AUTHOR: Jeff Buzas

#######################################





library("arm")  #will load  package for running WINbugs from R

library("R2WinBUGS")

library("BRugs")

#setwd("C:/Rbayes")



#  data

age=c(9, 9.5, 9.75, 10, 13, 9.5, 11, 10, 10, 11.75, 10.5, 15,      
      11.5, 12, 9, 11.5, 13.25, 13, 13.25, 11.5, 12, 13.5, 11.5, 12.35)


group= sort(rep(1:4,6))

groupf=as.factor( group   )





######## FIT FIXED EFFECTS MODEL ###########



mean(age) # grand mean

g = lm(age ~ groupf - 1) # fit cell means model, i.e. no intercept

summary(g)

#  anova(lm(age ~ group - 1) )





######## FIT RANDOM EFFECTS MODEL ###########

gr = lmer(age ~ (1 | groupf) )

summary(gr)

coef(gr)









########## BUGS CODE ########################

model<-function()
  
{
  
  for( i in 1 : n ) {
    
    
    
    age[i] ~ dnorm (mu[i],tau)
    
    mu[i] <- alpha[group[i]]
    
    
    
    
    
  }
  
  
  
  for( j in 1 : NT ) {
    
    
    
    alpha[j] ~ dnorm(mu.alpha,tau.alpha)
    
    
    
  }
  
  
  
  tau ~ dgamma(.001,.001)
  
  sigma <- 1/sqrt(tau)
  
  
  
  mu.alpha ~ dnorm(0,.01)
  
  # mu.alpha ~ dgamma(1,.1)
  
  tau.alpha ~ dgamma(.001,.001)
  
  sigma.alpha <- 1/sqrt(tau.alpha)
  
  
  
}  # END FUNCTION 'MODEL'





write.model(model, "oneway.bug")



file.show("oneway.bug")







########## END BUGS CODE #####################







NT=4  # number of groups



n=length(age)





data <- list(  "age", "n", "NT", "group" )





parameters <- c("alpha", "mu.alpha", "sigma.alpha","sigma")





# NOTE: HERE RANDOM STARTING VALUES ARE USED FOR ALL VARIABLES.



inits = function () {
  
  list( alpha=rnorm(NT), tau=rgamma(1,1,1), mu.alpha=rnorm(1), tau.alpha=rgamma(1,1,1)  )
  
}



print(date())


################################################################
## WinBUGS
################################################################

## If using WinBUGS - uncomment following lines

#oneway.sim <- bugs(data, inits, parameters, "oneway.bug", n.chains=3, n.iter=2000,                  
#                   bugs.directory="c:/winbugs14/WinBUGS14/",debug=T,working.directory="M:/work/courses/bayes course/fall2013/R programs")

#print(date())


################################################################
## OpenBUGS
################################################################

# the following code calls openbugs from my computer

oneway.sim <- openbugs(data, inits, parameters, "oneway.bug", n.chains=3, n.iter=2000)

##  if OpenBUGS was not installed globally with admin priviledges, you will need to add 
##   'bugs.directory = "path/to/OpenBUGS"' to the opengubs function call

#oneway.sim <- openbugs(data, inits, parameters, "oneway.bug", n.chains=3, n.iter=2000, 
#                       bugs.directory="path/to/OpenBUGS") 

print(oneway.sim,digits=3)


arm::traceplot(oneway.sim)  # available from 'arm' package


plot(oneway.sim, display.parallel = F)


plot(oneway.sim)


attach.bugs(oneway.sim)



############# COMPARE MEANS FROM DIFFERENT GROUPS ###############


dim(alpha)  


### THE FOLLOWING FUNCTION COMPUTES PERCENTILES OF THE POSTERIOR DISTRIBUTION

### FOR THE DIFFERENCE OF MEANS BETWEEN TWO TRT GROUPS

pairwise = function(x,y) {
  
  return(quantile(x-y,probs=c(.025,.05,.25,.5,.75,.95,.975)))
  
}



pairwise(alpha[,1],alpha[,2])


hist(alpha[,1]-alpha[,2])

plot(density(alpha[,1]-alpha[,2]),xlab="difference",ylab="Density",main="group difference")