MSE 598-DM Spring 2021

Bayesian exercise:

R code for the interested

First install JAGS (Windows, Mac OS X, Linux): http://mcmc-jags.sourceforge.net

We will use R package rjags to access JAGS within R:

install.packages("rjags")
library(rjags)

rattumor.txt has two column, y and N. y is the number of tumors and N is the number of control units.

d <- read.table("rattumor.txt", header=TRUE)
head(d)
summary(d)
with(d, hist(y/N, freq=FALSE, xlim=c(0,1)))

JAGS model specification in file rattumor1.bug:

model {
  for (j in 1:length(y)) {
    y[j] ~ dbin(theta[j], N[j])
    theta[j] ~ dbeta(alpha, beta)
}
  alpha ~ dexp(0.001)
  beta ~ dexp(0.001)
}

m <- jags.model("rattumor1.bug", d)
update(m, 2500) # burn-in
x <- coda.samples(m, c("alpha","beta"), n.iter=10000)
head(x)
head(as.matrix(x))
summary(x)
require(lattice)
densityplot(x) #posterior density of alpha and beta parameters