WinBUGS Examples

R‎ > ‎

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