# Bayesian Analysis in R/WinBUGS and SAS (Proc MCMC)

I have a Bayesian model here in R/WinBUGS. This data "pumps" has two columns, time and fail. Each observation is a pump and fail[1] tells how many failures there were in time[1]. I want to know how I can write this same model in SAS using Proc MCMC. Thank you!

pump <- function() {
for(i in 1:10){
fail[i] ~ dpois(lam[i])
theta[i] ~ dgamma(a,b)
lam[i] <- theta[i]*time[i]
}
a ~ dgamma(1.5,.25)
b ~ dgamma(1.5,.25)
}

filename <- file.path(tempdir(),'pump.bug')
write.model(pump,filename)
time <- pumps[,1]
fail <- pumps[,2]

data <- c('time','fail')
parameters <- c('fail','lam','theta','a','b')
pump.sim <- bugs(data, inits=NULL, parameters, model.file='pump.bug',
n.iter=20000, n.burnin=1000, n.chains=1, n.thin=1, debug=T)

• Could you do a little bit of explaining what the code does in order to make it clearer to those of us who might be good at SAS but not so much familiar with what the WinBUGS code is doing? On first pass, I'd suggest passing Proc MCMC and taking a look at Proc Lifereg with the MCMC option instead. Commented Mar 28, 2012 at 4:05
• Epi, it's all specified in the pump function: $a$ and $b$ are hyperparameters for $\theta$ (themselves with Gamma priors) and $\theta$ is a parameter for a Poisson model for the probability of failure (as a function of time).
– whuber
Commented Mar 28, 2012 at 15:33
• THis code does a Bayesian analysis of the survival time for a given pump to break. The data has 10 obs, so the first for loop is setting the model and priors for each obs. Then a and b are given hyper-priors. Time and Fail are the two variable in the data set. Hope this helps clarify.
– Lynn
Commented Apr 4, 2012 at 18:54

See if something like this works:

proc mcmc
data = /* your data here */
nmc = 20000, nbi = 1000; /* iterations and burn-in */

/* Specify parameters -- separate statements means they're proposed separately */
parms a;
parms b;
parms theta;

/* Specify priors */
hyper a ~ gamma(1.5, scale = .25);
hyper b ~ gamma(1.5, scale = .25);
prior theta ~ gamma(a, scale = b);

/* Sampling model */
lam = theta*time;
model fail ~ poisson(lam);
run;