Mixture Model using R2WinBUGS

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.

John Miyamoto.pdf (323.82 KB, 需要: 1 个论坛币)

2013-12-7 22:43:19 上传

Testing Bugs Installation with Examples
需要: 1 个论坛币  [购买]

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?
  

Two Way ANOVA using OpenBUGS and R John Miyamoto.pdf (402.05 KB, 需要: 1 个论坛币)

2013-12-7 23:16:40 上传

Two Way ANOVA using BUGS and R John Miyamoto
需要: 1 个论坛币  [购买]

hnd10-1.p548.spr13.zip (5.62 KB) 本附件包括：

• hnd10-1.p548.spr13.r

2013-12-7 23:26:31 上传

Twowayanova.r

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)

我来标记一下