Bayes Ideas and Data Analysis
Chapter 1 

Section 1.1:  Probability of a defective:  Binomial data

y=219 defectives
n=2430 Bernoulli trials

model{
 y ~ dbin(theta, 2430)
 theta ~ dbeta(12.05,116.06)
}
list(y=219)
list(theta=0.5)

Section 1.2:  Brass alloy zinc content:  Normal data

model{
 for(i in 1:n){ y[i] ~ dnorm(mu,tau) }
 mu ~ dnorm(4.75, prec)
 tau ~ dgamma(0.001, 0.001)
 y13 ~ dnorm(mu,tau)
 prob <- phi((4.4-mu)/sigma)
 prec <- 1/0.0163
 sigma <- 1/sqrt(tau)
}
list(n=12, y=c(4.20, 4.36, 4.11, 3.96, 5.63, 4.50, 5.64, 4.38,4.45, 3.67, 5.26, 4.66))
list(mu=0, tau=1, y13=1)

R code to get a histogram and normal plot of y
y=c(4.20, 4.36, 4.11, 3.96, 5.63, 4.50, 5.64, 4.38,4.45, 3.67, 5.26, 4.66)
par(mfrow=c(2,1))
hist(y)
qqnorm(y)
qqline(y)

Section 1.3:  Armadillo hunting:  Poisson data

model{
    for(i in 1:38){ kills[i] ~ dpois(theta) }
    theta ~ dgamma(1.11,1.61)
}
list(kills=c(2,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0))
list(theta=1)

Section 1.4:  Abortion in dairy cattle:  Survival data

Group 1 is infected cows
Group 2 is uninfected cows
model{
 for(i in 1:n1){ y1[i] ~ dlnorm(mu1,tau1) }
 for(i in 1:n2){ y2[i] ~ dlnorm(mu2,tau2) }
 mu1 ~ dnorm(0, 0.001)
 mu2 ~ dnorm(0, 0.001)
 tau1 ~ dgamma(0.001,0.001)
 tau2 ~ dgamma(0.001, 0.001)
 med1 <- exp(mu1)
 med2 <- exp(mu2)
 Delta <- med1 - med2
}
list(n1=19, n2=26, 
y1=c(50,130,100,130, 50,140,129,76,138,69,70,144,70,130,70,150,251,10,120),
y2=c(60,74,37,45,75,40,50,50,146,70,50,84,60,149,50,90,259,40,90,101,70,90,                  
      254, 130, 80, 40))
list(mu1=0, mu2=0, tau1=1, tau2=1)

Section 1.5:  Armadillo hunting with age trends:  Poisson regression for longitudinal count data

model{
 for(i in 1:38){
	kills[i]~dpois(rate.offset[i])
	rate.offset[i]<-lambda[i]*trials[i]
	log(lambda[i])<-a[1]+a[2]*(age[i]-50)+a[3]*pow(age[i]-50,2)+b[id[i]]
	cpo[i]<-1/(exp(-rate.offset[i]+kills[i]*log(rate.offset[i])-logfact(kills[i])))
	b[i]~dnorm(0,tau)
	e[i]<-kills[i]-rate.offset[i] }
 for(i in 1:3){ a[i] ~ dnorm(0,0.00001) }
 for(i in 20:70){ 
  b.new[i]~dnorm(0,tau)
  rate.new[i]<-a[1]+a[2]*(i-50)+a[3]*pow(i-50,2)+b.new[i] }	
tau <- 1/pow(sigma,2)
sigma ~ dunif(0,100)}

list(a=c(-1,0,0),sigma=1)

id[]  age[]  kills[]  trials[]
1    65.8      17    64
2    62.5      30   105
3    60.8       3    24
4    58.8      29    48
5    57.0       3     9
6    55.5       30   97
7    54.7       45  107
8    53.8       34   37
9    50.5       16   97
10   48.5       33   64
11   46.7       37   48
12   41.4       10   17
13   39.7       33   72
14   39.5        6   12
15   39.0        2   11
16   36.7       10   28
17   34.5       24   46
18   35.0        3   33
19   32.6        0    8
20   32.7       34   70
21   32.0       12   45
22   32.0        2   12
23   31.0        2   11
24   29.8        5   18
25   30.0        9   22
26   28.0        2    4
27   26.7        1   13
28   24.2        0   10
29   22.0        1   16
30   21.5        3   31
31   21.0        1    7
32   19.9        0   32
33   18.0        0    8
34   17.0        0    3
35   16.5        1   12
36   17.0        0    3
37   56.0       19   47
38   20.0        2   11
END

Section 1.6:  Lung cancer:  Log normal regression

model{
 for(i in 1:121){
  t[i] ~ dlnorm(mu[i], tau)I(left[i],)
  mu[i] <- beta[1]+beta[2]*a[i]+beta[3]*g[i]
 }
 tau ~ dgamma(0.001,0.001) 
 for(i in 1:3){ beta[i] ~ dnorm(0,0.001) }
 a.mean <- mean(a[1:121])
 mean[1] <- beta[1]+beta[2]*a.mean
 mean[2] <- beta[1]+beta[2]*(a.mean+10)
 mean[3] <- beta[1]+beta[2]*a.mean+beta[3]
 mean[4] <- beta[1]+beta[2]*(a.mean+10)+beta[3]
 for( j in 1:4){ S[j,1] <- 1 }
   for(i in 2:100){ 
    grid[i] <- (i-1)*20 
    for(j in 1:4){ S[j,i] <- 1-phi((log(grid[i])-mean[j])*sqrt(tau)) }
   }
}
list(beta=c(10,0,0), tau=1) # Initial values

Section 1.7:  Armadillo hunting:  Survival with random effects

model{
	for(i in 1:324){
		SEARCH[i] ~ dlnorm(mean[i],tau)I(LOWER[i],)
		mean[i] <- mu+delta[ID[i]]
		dummy[i] <- FIRSTLAST[i]+AGE[i]  # not used in simple analysis in Ch 1
     }
     for(i in 1:15){ delta[i] ~ dnorm(0,tau.d) }  # delta[15] is predictive
	mu ~ dnorm(4.5,10)
	tau <- 1/(sigma*sigma)
	tau.d <- 1/(sigma.d*sigma.d)
	sigma ~ dgamma(2.29,2.92)
	sigma.d ~ dgamma(8.1,24.5)
	mean.pred <- mu+delta[15]
	t.pred ~ dlnorm(mean.pred,tau)
	for(i in 2:100){ grid[i] <- i*2 }
	grid[1] <- 0.01
	for(i in 1:100){ d[i] <- exp(-0.5*pow(log(grid[i])-mean.pred,2)*tau)/(sqrt(2*3.14159/tau)*grid[i]) }
}
list(mu = 4, sigma=1, sigma.d=1)