R‎ > ‎

WinBUGS Examples

These are a series of WinBUGS models that I have run for a variety of models. Use at your own risk. See my other post on how to get WinBUGS going under WINE.

Hierarchical Normal Models:
1. The normal linear model with varying intercepts
2. The normal linear model with varying intercepts, some parameters varying slopes

Set up the data as follows:
Xvec=c(X1,X2...Xn) #create a vector for parameters that will have constant slopes over all levels
XVec2=c(XA1,XA1...XAn) #create a vector for parameters that will have slopes at the individual level
P=length(Xvec)
P2=length(Xvec2)
J=number of levels
station=c(1,1,1,2,2,2,3,3,3....) #indicator for what group the observation belongs to
N=length(Y)

Set up Winbugs to run under WINE
WINE="/opt/local/bin/wine"
WINEPATH="/opt/local/bin/winepath"
BUGS.DIR="/Users/bigbird/.wine/drive_c/Program Files/WinBUGS14/"

1. Varying intercepts only
for j in 1:J
for n in 1:N
Y=aj + B1x1 + B2x2 + Bnxn 
fishmodel <- function(){
    for (n in 1:N){
                     muj[n] <- alpha[station[n]] + inprod(beta[1:P],Xvec[n,1:P]) #the predictors
                      Y2[n] ~ dnorm(muj[n], phi) #normal distribution
                   }
    phi ~ dgamma(1.0E-6,1.0E-6) #prior on the variance of the D.Var
 
    for(j in 1:P){beta[j] ~ dnorm(0, 1.0E-6)} #priors on the Betas

    for (j in 1:J){alpha[j] ~ dnorm(mu, phi.mu) #individual intercepts depend on an overall intercept
                 }
    mu ~ dnorm (0.0, 1.0E-6) #priors on the overall intercepts
    phi.mu <- pow(sigma.mu, -2)
    sigma.mu ~ dunif (0, 1000)

    }
 inits = function() {list(mu=5,alpha=rnorm(J,5,sqrt(3)),phi=1/39.15,sigma.mu=sqrt(8.61),beta=rep(1,P))} #add initialization variables
 parameters.to.save = c("alpha","beta") #what parameters to save
data=list(N=N,P=P,Xvec=Xvec,Y2=Y2,J=J,station=station)

# write the model code out to a file
write.model(fishmodel, "nummodel.txt")
model.file = paste(getwd(),"nummodel.txt", sep="/")
## and let's take a look:
file.show("nummodel.txt")

sim = bugs(data, inits, parameters.to.save, model.file=model.file, n.chains=2, n.iter=5000, bugs.dir=BUGS.DIR, WINE=WINE, WINEPATH=WINEPATH, debug=T, DIC=F)
print(sim)
plot(sim)

2. Hybrid Model, varying intercepts, varying slopes, constant slopes
for j in 1:J
for n in number of cases
Y=aj + B1jx1 + B2jx2 +...+ Bnjxn + G1x1 + G2x2 +...+ Bnxn 
Xvec=cbind(X1,X2,X3) #can be as many or as few as you want - these vary at each level of the model
P=dim(Xvec)[2]
Xvec2=cbind(X4,X5,X6,X7,X8,X9,X10,X11,X12,X13,X14,X15,X16,X17,X18,X19,X20,X21,X22,X23) #can be as many as you want, these are constant across all levels
P2=dim(Xvec2)[2]
station=c(1,1,1,2,2,2,3,3,3,...) #indicator for what level case belongs to
fishmodel <- function(){
    for (n in 1:N){
                     muj[n] <- alpha[station[n]] + inprod(beta2[1:P2],Xvec2[n,1:P2]) + inprod(beta[station[n],1:P],Xvec[n,1:P])
                      Y2[n] ~ dnorm(muj[n], phi)
                   }
    phi ~ dgamma(1.0E-6,1.0E-6) #prior on the DVar slope

    for(p in 1:P){
    for(j in 1:J){
      beta[j,p] ~ dnorm(mub[p], phi.mub[p]) #each predictor is a function of an overall slope
     }
    mub[p] ~ dnorm (0.0, 1.0E-6) #these are the overall slopes for each predictor, normal prior
    phi.mub[p] <- pow(sigma.mub[p], -2) #prior on the overall variances
    sigma.mub[p] ~ dunif (0, 1000)
    }
    
    for(j in 1:P2){beta2[j] ~ dnorm(0, 1.0E-6)} #priors for predictors that don't vary across the predictors
    for (j in 1:J){alpha[j] ~ dnorm(mu, phi.mu)} #intercepts at each level depend on an overall intercept

    mu ~ dnorm (0.0, 1.0E-6) #prior on the overall intercept
    phi.mu <- pow(sigma.mu, -2)
    sigma.mu ~ dunif (0, 1000)

    }
inits = function() {list(mu=5,alpha=rnorm(J,5,sqrt(3)),phi=1/39.15,sigma.mu=sqrt(8.61),beta2=rep(1,P2),beta=matrix(rep(1,J*P),nrow=J,ncol=P))}
parameters.to.save = c("alpha","beta","beta2")
data=list(N=N,P=P,Xvec=Xvec,Y2=Y2,J=J,station=station,P2=P2,Xvec2=Xvec2)  # include only variables in the model
# write the model code out to a file
write.model(fishmodel, "nummodel.txt")
model.file = paste(getwd(),"nummodel.txt", sep="/")
## and let's take a look:
file.show("nummodel.txt")

sim = bugs(data, inits, parameters.to.save, model.file=model.file, n.chains=2, n.iter=5000, bugs.dir=BUGS.DIR, WINE=WINE, WINEPATH=WINEPATH, debug=T, DIC=F)
print(sim)
plot(sim)
Comments