One and Two Population Inference:  WinBUGS code 
Chapter 5

Section 5.1: Proportions

Example 5.1.3. Comparing two proportions.  Medfly data
Group 1 = high producers
Group 2 = low producers

model
{
y1 ~ dbin(theta1,n1)
y2 ~ dbin(theta2,n2)
theta1 ~ dbeta(1,1)
theta2 ~ dbeta(1,1)
odds1 <- theta1/(1-theta1)
odds2 <- theta2/(1-theta2)
RD <- theta2-theta1
RR <- theta2/theta1
OR <- odds2/odds1
test <- step(RD)
}
list(n1=134,y1=54,n2=400,y2=224)
list(theta1=0.5,theta2=0.5)

Exercise 5.15. Case-control data for Reye's syndrome

model{
 for(i in 1:2){ y[i] ~ dbin(theta[i],n[i]) }
 theta[2] ~ dbeta(a,b)   
 delta ~ dnorm(mu, prec)   
 theta[1] <- exp(delta)*theta[2]/(1-theta[2]*(1-exp(delta)))   
 OR <- theta[1]/(1-theta[1])/(theta[2]/(1-theta[2])) 
}
list(y=c(7,8), n=c(7,16), mu=0,prec=0.5,a=0.5,b=0.5)
list(delta=0)

Section 5.2: Normal data

Exercise 5.23: Brass alloy zinc data

model{ 
 for(i in 1:n) { y[i] ~ dnorm(mu, tau) }
 mu ~ dflat()         
 tau ~ dgamma(c,d)
 sigma <- 1/sqrt(tau)
 gamma <- phi((4.4-mu)/sqrt(1/tau))
 prob <- step(4.4 - y[13])   
}
list(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, NA), c=0.001, d=0.001,n=13)
list(mu=0,tau =1)

Example 5.2.2. Comparing two log normal medians:  Diasorin data

model
 {
 for(i in 1:nL){
   low[i] ~ dlnorm(muL, tauL)
   }
 for(i in 1:nN){
   normal[i] ~ dlnorm(muN, tauN)
   }
 muL ~ dnorm(4.87, 347.12)
 muN ~ dnorm(5.39, 357.14)
 #muL ~ dnorm(4.87, 0.001)
 #muN ~ dnorm(5.39, 0.001)
 tauL ~ dgamma(1.0376,0.001)
 tauN ~ dgamma(1.04653,0.001)
 medianL <- exp(muL)
 medianN <- exp(muN)
 relmedian <- medianN / medianL
 test <- step(medianN-medianL)
 lowf ~ dlnorm(muL, tauL)
 normalf ~ dlnorm(muN, tauN)
 diffmu <- muN-muL
 ratiotau <- tauN/tauL
 }
 list(nL=19,nN=15,
 low = c(91, 46,  95,  60, 33, 410, 105, 43, 189,1097, 54,178, 114, 137, 233, 101, 
         25,70,357),
 normal = c(370, 267, 99,157, 75,1281, 48, 298, 268, 62,804,430,171,694,404))
 
list(muL=4.87, muN=5.39, tauL=1, tauN=1, lowf=50, normalf=50)
list(muL=0,      muN=0,      tauL=5,     tauN=15,    lowf=50, normalf=50)
list(muL=-5,     muN=-5,     tauL=10,  tauN=15,    lowf=50, normalf=50)
list(muL=10,    muN=10,    tauL=15,  tauN=10,    lowf=50, normalf=50)

Results

With prior  muL ~ dnorm(4.87, 0.001); muN ~ dnorm(5.39, 0.001)

	 node	 mean	 sd	     MC error	2.5%	  median	97.5%	  start	sample
	lowf	      178.4	275.1	0.2599	      16.37	  110.5	    744.4	  10000	1160004
	medianL	113.1	24.81	0.0228	      72.11	  110.5	    169.3	  10000	1160004
	medianN	249.9	63.89	0.06138	    147.6	  242.3	     397.4	 10000	1160004
	muL	       4.705	0.2157   1.983E-4	   4.278	  4.705	    5.132	  10000	1160004
	muN	       5.49	0.2492	2.39E-4	     4.994	  5.49	    5.985	   10000	1160004
	normalf	 408.0	958.7	 0.8724	      33.65	  242.3	    1747.0	10000	1160004
	tauL	       1.256	0.3963   3.795E-4	   0.6039	1.214	   2.144	   10000	1160004
	tauN	       1.227	0.4323   4.229E-4	   0.5318	1.176	   2.207	   10000	1160004

With prior muL ~ dnorm(4.87, 347.12);  muN ~ dnorm(5.39, 357.42)

	 node	 mean	 sd	       MC error	2.5%	median	97.5%	start	sample
	diffmu	0.5358	0.07324	6.915E-5	 0.3925	0.5359	 0.6794   10000	1160004
	lowf	    201.2	276.6	    0.2588	     20.46	  128.8	   813.0	10000	1160004
	medianL	129.1	6.716	  0.006184	116.4	  128.9	    142.8	10000	1160004
	medianN   220.6	11.38	    0.01083	  199.1	  220.3	    243.7	10000	1160004
	muL	     4.859	0.05198	4.786E-5	4.757	  4.859	    4.961	10000	1160004
	muN	     5.395	0.05157	4.905E-5	5.294	  5.395	    5.496	10000	1160004
	normalf	348.3	556.1	    0.5056	   34.19	  220.7	    1415.0    10000	1160004
	ratiotau	1.111	0.5532	  5.211E-4	0.3873	0.9985	  2.494	10000	1160004
	relmedian  1.713	0.1257	  1.185E-4	1.481	  1.709	    1.973	10000	1160004
	tauL	     1.279	0.3945	  3.638E-4	0.6269	1.238	    2.16	   10000	1160004
	tauN	     1.286	0.4397	  4.083E-4	0.5741	1.236	    2.279	10000	1160004
	test	      1.0	    0.0	        4.642E-14	1.0	    1.0	        1.0	     10000	1160004

Example 5.2.2.  Comparing two normal means using log transformed data:  Diasorin data

 model
 {
 for(i in 1:nL){
   loglow[i] <- log(low[i])
   loglow[i] ~ dnorm(muL, tauL)
   }
 for(i in 1:nN){
   lognormal[i] <- log(normal[i])
   lognormal[i] ~ dnorm(muN, tauN)
   }
 #muL ~ dnorm(4.87, 347.12)
 #muN ~ dnorm(5.39, 357.14)
 muL ~ dnorm(0,0.001)
 muN ~ dnorm(0,0.001)
 tauL ~ dgamma(1.0376,0.001)
 tauN ~ dgamma(1.04653,0.001)
 loglowf ~ dnorm(muL, tauL)
 lognormalf ~ dnorm(muN, tauN)
 diffmu <- muN-muL
 ratiotau <- tauN/tauL
 }
 list(nL=19,nN=15,
 low = c(91, 46,  95,  60, 33, 410, 105, 43, 189,1097, 54,178, 114, 137, 233, 101, 
         25,70,357),
 normal = c(370, 267, 99,157, 75,1281, 48, 298, 268, 62,804,430,171,694,404))
 list(muL=4.87, muN=5.39, tauL=1, tauN=1, loglowf=50, lognormalf=50)

Results

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	diffmu	0.5356	0.07314	3.407E-4	0.3922	0.5356	0.6792	10001	50000
	loglowf	4.852	0.9257	0.003869	3.022	4.85	6.696	10001	50000
	lognormalf	5.393	0.9333	0.004272	3.539	5.394	7.253	10001	50000
	muL	4.86	0.05213	2.346E-4	4.757	4.86	4.961	10001	50000
	muN	5.395	0.05141	2.383E-4	5.294	5.395	5.496	10001	50000
	ratiotau	1.114	0.5547	0.002562	0.3873	1.003	2.496	10001	50000
	tauL	1.275	0.3944	0.001796	0.6245	1.231	2.161	10001	50000
	tauN	1.285	0.4389	0.00195	0.5757	1.236	2.274	10001	50000

Sensitivity analysis:  Diasorin data

 model
 {
 for(i in 1:nL){
   low[i] ~ dlnorm(muL, tauL)
   }
 for(i in 1:nN){
   normal[i] ~ dlnorm(muN, tauN)
   }
 muL ~ dnorm(0, 0.01)
 muN ~ dnorm(0, 0.01)
 tauL ~ dgamma(0.001,0.001)
 tauN ~ dgamma(0.001,0.001)
 medianL <- exp(muL)
 medianN <- exp(muN)
 relmedian <- medianN / medianL
 test <- step(medianN-medianL)
 }
 list(nL=19,nN=15,
 low = c(91, 46,  95,  60, 33, 410, 105, 43, 189,1097, 54,178, 114, 137, 233, 101, 
        25,70,357),
 normal = c(370, 267, 99,157, 75,1281, 48, 298, 268, 62,804,430,171,694,404))
 list(muL=4.87, muN=5.39, tauL=1, tauN=1)

Results

Original analysis
     node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	medianL	129.1	6.675	0.02604	116.5	128.9	142.7	10000	50001
	medianN	220.6	11.39	0.05308	199.2	220.4	243.7	10000	50001
	muL	4.859	0.05167	2.027E-4	4.757	4.859	4.96	10000	50001
	muN	5.395	0.0516	2.401E-4	5.294	5.395	5.496	10000	50001
	relmedian	1.714	0.1252	5.693E-4	1.482	1.71	1.972	10000	50001
	tauL	1.279	0.3971	0.001763	0.6266	1.236	2.168	10000	50001
	tauN	1.285	0.4405	0.001885	0.5737	1.234	2.277	10000	50001
	test	1.0	0.0	4.472E-13	1.0	1.0	1.0	10000	50001

Sensitivity analysis
 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	medianL	113.0	26.16	0.1039	70.17	110.2	172.5	10000	50001
	medianN	250.8	70.1	0.3409	141.4	241.6	411.9	10000	50001
	muL	4.702	0.2274	9.254E-4	4.251	4.703	5.15	10000	50001
	muN	5.488	0.2703	0.001317	4.952	5.487	6.021	10000	50001
	relmedian	2.337	0.8609	0.004054	1.1	2.193	4.402	10000	50001
	tauL	1.128	0.3778	0.001702	0.5157	1.083	1.979	10000	50001
	tauN	1.064	0.4018	0.001971	0.4295	1.013	1.987	10000	50001
	test	0.9861	0.1169	5.05E-4	1.0	1.0	1.0	10000	50001

Exercise 5.29

model{ 
 for(i in 1:n[1]) { y[ i ] ~ dnorm(mu[1], tau[1]) }
 for(j in 1:n[2]) { x[ j ] ~ dnorm(mu[2], tau[2]) }
 for(r in 1:2){
   mu[ r ] ~ dnorm(a[ r ],b[ r ])
   tau[r] ~ dgamma(c[ r ],d[ r ])
  sigma[ r ] <- sqrt(1/tau[ r ]) 
 }
 meandiff <- mu[1] - mu[2]
 sdratio <- sigma[1]/sigma[2]
 prob[1] <- step(meandiff) # Gives Pr(meandiff >0|data)
 prob[2] <- step(sdratio -1)# Gives Pr(sdratio >1|data)
}
list(y = c(.,...,.), x = c(.,...,.), n = c(.,.),
a = c(.,.), b = c(.,.), c = c(.,.), d = c(.,.))
list(mu = c(0,0), tau = c(1,1))

Exercise 5.30. Diasorin data

model{ 
 for(i in 1:n[1]){ low[i] ~ dlnorm(mu[1], tau[1]) }
 for(i in 1:n[2]){ normal[i] ~ dlnorm(mu[2], tau[2]) }                                                                  # mu[1] ~ dnorm(4.87, 347.12)
 # mu[2] ~ dnorm(5.39, 357.14)
 mu[1] ~ dnorm(0,0.00001)
 mu[2] ~ dnorm(0,0.00001)
 # tau[1] ~ dgamma(1.0376,0.001)
 # tau[2] ~ dgamma(1.04653,0.001)
 tau[1] ~ dgamma(0.001,0.001)
 tau[2] ~ dgamma(0.001,0.001)
 med[1] <- exp(mu[1])
 med[2] <- exp(mu[2])
 rmed <- med[2] / med[1]
 test[1] <- step(med[2]-med[1])
 test[2] <- step(Nf - Lf)
 Lf ~ dlnorm(mu[1], tau[1])
 Nf ~ dlnorm(mu[2], tau[2])
 dmu <- mu[2]-mu[1]
 rtau <- tau[2]/tau[1] 
}
list(n=c(19,15),
     low = c(91,46,95,60,33,410,105,43,189,1097,54,178,114,137,233,101,25,70,357),
     normal = c(370,267,99,157,75,1281,48,298,268,62,804,430,171,694,404))
list(mu=c(0,0), tau=c(1,1), Lf=50, Nf=50)

Section 5.3: Count data

Example 5.3.1:  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)

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	theta	0.2803	0.08457	7.967E-4	0.1402	0.2717	0.469	1	10000

Exercise 5.31. Armadillo data

model{
 for(i in 1:38){ kills[i] ~ dpois(theta) }
 theta ~ dgamma(1.11,1.61)
 prob <- exp(-theta)
}
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))

Example 5.3.2.  Comparing two rates.  Breast cancer incidence data

model{
 y1 ~ dpois(lambda1)
 y2 ~ dpois(lambda2)
 lambda1 <- theta1*M1
 lambda2 <- theta2*M2
 theta1 ~ dgamma(0.5,0.001)
 theta2 ~ dgamma(0.5,0.001)
 r1 <- 1000*theta1
 r2 <- 1000*theta2
 RateRatio <- theta1/theta2
 test <-  step(RateRatio-1)
 }
 list(y1=123, M1=46524, y2=288, M2=145159)
 list(theta1=1, theta2=1)

Results

 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	RateRatio	1.34	0.1447	4.585E-4	1.076	1.333	1.642	10001	100000
	r1	2.654	0.2387	7.39E-4	2.208	2.647	3.144	10001	100000
	r2	1.987	0.1172	3.893E-4	1.765	1.984	2.225	10001	100000
	test	0.9958	0.06459	2.191E-4	1.0	1.0	1.0	10001	100000
	theta1	0.002654	2.387E-4	7.39E-7	0.002208	0.002647	0.003144	10001	100000
	theta2	0.001987	1.172E-4	3.893E-7	0.001765	0.001984	0.002225	10001	100000