WinBUGS Code for Chapter 11:  Count Data

Fabric fault data

1a.  Bissell (1972) data.   MODEL 1
      y=number of faults on a piece of fabric
      M=length of the piece of fabric
      n=32
model
{
for(i in 1:32){
   y[i] ~ dpois(lambda[i])
   lambda[i] <- theta * M[i]
   }
theta ~ dgamma(0.001,0.001)
junk <- ID[1]
}

DATA

ID[ ] M[ ]  y[ ] 
 1 551  6 
 2 651  4 
 3 832 17 
 4 375  9 
 5 715 14
 6 868  8 
 7 271  5 
 8 630  7 
 9 491  7 
10 372  7 
11 645  6 
12 441  8 
13 895 28 
14 458  4 
15 642 10 
16 492  4 
17 543  8 
18 842  9 
19 905 23 
20 542  9 
21 522  6 
22 122  1 
23 657  9 
24 170  4 
25 738  9 
26 371 14 
27 735 17 
28 749 10 
29 495  7 
30 716  3  
31 952  9  
32 417  2 
END

INITS
list(theta=1)

RESULTS

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
	         Dbar	     Dhat	     pD	     DIC	
y	    188.836	187.836	1.000	189.835	
total	188.836	187.836	1.000	189.835	

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	theta	0.0151	8.965E-4	2.863E-6	0.01339	0.01508	0.01692	1	100000
	lambda[1]	8.319	0.494	0.001577	7.377	8.309	9.324	1	100000
	lambda[2]	9.829	0.5836	0.001864	8.716	9.817	11.02	1	100000
	lambda[3]	12.56	0.7459	0.002382	11.14	12.55	14.08	1	100000
	lambda[4]	5.662	0.3362	0.001074	5.021	5.655	6.346	1	100000
	lambda[5]	10.8	0.641	0.002047	9.573	10.78	12.1	1	100000
	lambda[6]	13.11	0.7781	0.002485	11.62	13.09	14.69	1	100000
	lambda[7]	4.092	0.2429	7.758E-4	3.628	4.086	4.586	1	100000
	lambda[8]	9.512	0.5648	0.001804	8.435	9.5	10.66	1	100000
	lambda[9]	7.413	0.4402	0.001406	6.574	7.404	8.308	1	100000
	lambda[10]	5.617	0.3335	0.001065	4.981	5.609	6.295	1	100000
	lambda[11]	9.738	0.5782	0.001847	8.636	9.726	10.91	1	100000
	lambda[12]	6.658	0.3953	0.001263	5.904	6.65	7.462	1	100000
	lambda[13]	13.51	0.8023	0.002562	11.98	13.5	15.14	1	100000
	lambda[14]	6.915	0.4106	0.001311	6.132	6.906	7.75	1	100000
	lambda[15]	9.693	0.5755	0.001838	8.596	9.681	10.86	1	100000
	lambda[16]	7.428	0.4411	0.001409	6.587	7.419	8.325	1	100000
	lambda[17]	8.198	0.4868	0.001555	7.27	8.188	9.188	1	100000
	lambda[18]	12.71	0.7548	0.002411	11.27	12.7	14.25	1	100000
	lambda[19]	13.66	0.8113	0.002591	12.12	13.65	15.31	1	100000
	lambda[20]	8.183	0.4859	0.001552	7.257	8.173	9.171	1	100000
	lambda[21]	7.881	0.468	0.001494	6.989	7.871	8.833	1	100000
	lambda[22]	1.842	0.1094	3.493E-4	1.633	1.84	2.064	1	100000
	lambda[23]	9.92	0.589	0.001881	8.796	9.907	11.12	1	100000
	lambda[24]	2.567	0.1524	4.867E-4	2.276	2.563	2.877	1	100000
	lambda[25]	11.14	0.6616	0.002113	9.881	11.13	12.49	1	100000
	lambda[26]	5.601	0.3326	0.001062	4.967	5.594	6.278	1	100000
	lambda[27]	11.1	0.6589	0.002104	9.841	11.08	12.44	1	100000
	lambda[28]	11.31	0.6715	0.002144	10.03	11.29	12.67	1	100000
	lambda[29]	7.474	0.4438	0.001417	6.627	7.464	8.376	1	100000
	lambda[30]	10.81	0.6419	0.00205	9.586	10.8	12.12	1	100000
	lambda[31]	14.37	0.8534	0.002725	12.75	14.36	16.11	1	100000
	lambda[32]	6.296	0.3738	0.001194	5.583	6.288	7.056	1	100000

1b.  Bissell (1972) data.    MODEL 2
      y=number of faults on a piece of fabric
      M=length of the piece of fabric
      n=32
model
{
for(i in 1:32){
   y[i] ~ dpois(lambda[i])
   log(lambda[i]) <- beta1 + log(M[i])
   }
beta1 ~ dnorm(0,0.001)
theta <- exp(beta1)
junk <- ID[1]
}
DATA

ID[ ] M[ ]  y[ ] 
 1 551  6 
 2 651  4 
 3 832 17 
 4 375  9 
 5 715 14
 6 868  8 
 7 271  5 
 8 630  7 
 9 491  7 
10 372  7 
11 645  6 
12 441  8 
13 895 28 
14 458  4 
15 642 10 
16 492  4 
17 543  8 
18 842  9 
19 905 23 
20 542  9 
21 522  6 
22 122  1 
23 657  9 
24 170  4 
25 738  9 
26 371 14 
27 735 17 
28 749 10 
29 495  7 
30 716  3  
31 952  9  
32 417  2 
END

INITS
list(beta1=0)

RESULTS

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
y		188.832	187.837	0.996	189.828	
total	188.832	187.837	0.996	189.828	

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	beta1	-4.195	0.05927	1.834E-4	-4.312	-4.194	-4.081	10001	100000
	lambda[1]	8.321	0.4928	0.001528	7.387	8.312	9.311	10001	100000
	lambda[2]	9.831	0.5823	0.001806	8.727	9.82	11.0	10001	100000
	lambda[3]	12.56	0.7441	0.002308	11.15	12.55	14.06	10001	100000
	lambda[4]	5.663	0.3354	0.00104	5.027	5.657	6.337	10001	100000
	lambda[5]	10.8	0.6395	0.001983	9.585	10.79	12.08	10001	100000
	lambda[6]	13.11	0.7763	0.002407	11.64	13.09	14.67	10001	100000
	lambda[7]	4.092	0.2424	7.516E-4	3.633	4.088	4.579	10001	100000
	lambda[8]	9.514	0.5635	0.001747	8.446	9.503	10.65	10001	100000
	lambda[9]	7.415	0.4392	0.001362	6.582	7.407	8.297	10001	100000
	lambda[10]	5.618	0.3327	0.001032	4.987	5.611	6.286	10001	100000
	lambda[11]	9.74	0.5769	0.001789	8.647	9.73	10.9	10001	100000
	lambda[12]	6.66	0.3944	0.001223	5.912	6.652	7.452	10001	100000
	lambda[13]	13.52	0.8005	0.002482	12.0	13.5	15.12	10001	100000
	lambda[14]	6.916	0.4096	0.00127	6.14	6.909	7.739	10001	100000
	lambda[15]	9.695	0.5742	0.001781	8.607	9.684	10.85	10001	100000
	lambda[16]	7.43	0.44	0.001365	6.596	7.422	8.314	10001	100000
	lambda[17]	8.2	0.4857	0.001506	7.28	8.191	9.176	10001	100000
	lambda[18]	12.72	0.7531	0.002335	11.29	12.7	14.23	10001	100000
	lambda[19]	13.67	0.8094	0.00251	12.13	13.65	15.29	10001	100000
	lambda[20]	8.185	0.4848	0.001503	7.266	8.176	9.159	10001	100000
	lambda[21]	7.883	0.4669	0.001448	6.998	7.874	8.821	10001	100000
	lambda[22]	1.842	0.1091	3.384E-4	1.636	1.84	2.062	10001	100000
	lambda[23]	9.921	0.5876	0.001822	8.808	9.911	11.1	10001	100000
	lambda[24]	2.567	0.152	4.715E-4	2.279	2.564	2.873	10001	100000
	lambda[25]	11.14	0.6601	0.002047	9.894	11.13	12.47	10001	100000
	lambda[26]	5.602	0.3318	0.001029	4.974	5.596	6.269	10001	100000
	lambda[27]	11.1	0.6574	0.002039	9.853	11.09	12.42	10001	100000
	lambda[28]	11.31	0.6699	0.002077	10.04	11.3	12.66	10001	100000
	lambda[29]	7.475	0.4427	0.001373	6.636	7.467	8.365	10001	100000
	lambda[30]	10.81	0.6404	0.001986	9.599	10.8	12.1	10001	100000
	lambda[31]	14.38	0.8515	0.00264	12.76	14.36	16.09	10001	100000
	lambda[32]	6.297	0.373	0.001157	5.59	6.29	7.046	10001	100000
	theta	0.0151	8.944E-4	2.773E-6	0.01341	0.01508	0.0169	10001	100000

1c.  Bissell (1972) data.  MODEL 3
model
{
for(i in 1:32){
   y[i] ~ dpois(lambda[i])
   log(lambda[i]) <- beta1 + beta2*log(M[i])
   }
beta1 ~ dnorm(0,0.001)
beta2 ~ dnorm(0,0.001)
testslope1 <- step(beta2)
testslope2 <- step(beta2-1)
junk <- ID[1]
}

DATA

ID[ ] M[ ]  y[ ] 
 1 551  6 
 2 651  4 
 3 832 17 
 4 375  9 
 5 715 14
 6 868  8 
 7 271  5 
 8 630  7 
 9 491  7 
10 372  7 
11 645  6 
12 441  8 
13 895 28 
14 458  4 
15 642 10 
16 492  4 
17 543  8 
18 842  9 
19 905 23 
20 542  9 
21 522  6 
22 122  1 
23 657  9 
24 170  4 
25 738  9 
26 371 14 
27 735 17 
28 749 10 
29 495  7 
30 716  3  
31 952  9  
32 417  2 
END

INITS
list(beta1=0, beta2=0)

RESULTS

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
y		189.808	187.839	1.969	191.777	
total	189.808	187.839	1.969	191.777	

 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	beta1	-4.17	1.119	0.05432	-6.547	-4.115	-2.039	10001	100000
	beta2	0.9959	0.1734	0.008416	0.6659	0.9877	1.364	10001	100000
	testslope1	1.0	0.0	3.162E-13	1.0	1.0	1.0	10001	100000
	testslope2	0.47	0.4991	0.02065	0.0	0.0	1.0	10001	100000

	lambda[1]	8.314	0.5305	0.01006	7.304	8.303	9.387	10001	100000
	lambda[2]	9.814	0.5863	0.002731	8.699	9.802	11.0	10001	100000
	lambda[3]	12.54	0.961	0.02858	10.75	12.52	14.52	10001	100000
	lambda[4]	5.689	0.609	0.02491	4.521	5.679	6.921	10001	100000
	lambda[5]	10.78	0.6829	0.01085	9.484	10.76	12.16	10001	100000
	lambda[6]	13.09	1.067	0.03452	11.11	13.05	15.29	10001	100000
	lambda[7]	4.144	0.6469	0.02915	2.935	4.121	5.503	10001	100000
	lambda[8]	9.499	0.5649	0.001982	8.426	9.486	10.64	10001	100000
	lambda[9]	7.417	0.544	0.01603	6.375	7.41	8.514	10001	100000
	lambda[10]	5.644	0.6106	0.02509	4.474	5.634	6.88	10001	100000
	lambda[11]	9.724	0.5796	0.002229	8.622	9.712	10.89	10001	100000
	lambda[12]	6.671	0.5699	0.02032	5.576	6.665	7.819	10001	100000
	lambda[13]	13.5	1.151	0.03912	11.38	13.46	15.89	10001	100000
	lambda[14]	6.925	0.5602	0.01894	5.848	6.917	8.052	10001	100000
	lambda[15]	9.679	0.5764	0.002044	8.583	9.666	10.84	10001	100000
	lambda[16]	7.432	0.5436	0.01593	6.391	7.425	8.529	10001	100000
	lambda[17]	8.195	0.5307	0.0109	7.182	8.183	9.266	10001	100000
	lambda[18]	12.7	0.9895	0.03021	10.85	12.67	14.73	10001	100000
	lambda[19]	13.65	1.183	0.04085	11.48	13.61	16.11	10001	100000
	lambda[20]	8.18	0.5308	0.01101	7.166	8.168	9.251	10001	100000
	lambda[21]	7.881	0.5338	0.01305	6.859	7.87	8.958	10001	100000
	lambda[22]	1.927	0.5535	0.02606	0.9972	1.875	3.205	10001	100000
	lambda[23]	9.904	0.5933	0.003344	8.777	9.892	11.1	10001	100000
	lambda[24]	2.643	0.6124	0.02861	1.565	2.601	4.007	10001	100000
	lambda[25]	11.12	0.7283	0.01412	9.747	11.11	12.6	10001	100000
	lambda[26]	5.629	0.6112	0.02514	4.458	5.619	6.867	10001	100000
	lambda[27]	11.08	0.7221	0.01369	9.713	11.06	12.54	10001	100000
	lambda[28]	11.29	0.7518	0.01573	9.87	11.27	12.81	10001	100000
	lambda[29]	7.477	0.5424	0.01565	6.437	7.469	8.571	10001	100000
	lambda[30]	10.79	0.6848	0.01099	9.495	10.78	12.18	10001	100000
	lambda[31]	14.37	1.341	0.0492	11.92	14.3	17.19	10001	100000
	lambda[32]	6.314	0.5842	0.02214	5.192	6.307	7.49	10001	100000

	

1d.  Example 11.2.2  
      Bissell (1972) data; with overdispersion.  
      y=number of faults on a piece of fabric
      M=length of the piece of fabric
      n=32
model
{
for(i in 1:32){
   y[i] ~ dpois(lambda[i])
   log(lambda[i]) <- beta1 + log(M[i]) + delta[i]
   delta[i] ~ dnorm(0,tau)
   }
beta1 ~ dnorm(0,0.001)
tau ~ dgamma(0.001, 0.001)
junk <- ID[1]
}

DATA

ID[ ] M[ ]  y[ ] 
 1 551  6 
 2 651  4 
 3 832 17 
 4 375  9 
 5 715 14
 6 868  8 
 7 271  5 
 8 630  7 
 9 491  7 
10 372  7 
11 645  6 
12 441  8 
13 895 28 
14 458  4 
15 642 10 
16 492  4 
17 543  8 
18 842  9 
19 905 23 
20 542  9 
21 522  6 
22 122  1 
23 657  9 
24 170  4 
25 738  9 
26 371 14 
27 735 17 
28 749 10 
29 495  7 
30 716  3  
31 952  9  
32 417  2 
END

INITS
list(beta1=0, tau=1)

RESULTS

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
		Dbar			Dhat		pD		DIC	
y		155.342	138.959	16.383	171.725	
total	155.342	138.959	16.383	171.725	

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	beta1	-4.252	0.09162	6.165E-4	-4.441	-4.25	-4.08	10001	100000
	lambda[1]	7.141	1.843	0.00662	3.97	6.987	11.2	10001	100000
	lambda[2]	6.932	1.894	0.009749	3.632	6.798	11.0	10001	100000
	lambda[3]	15.01	3.149	0.01162	9.73	14.7	22.0	10001	100000
	lambda[4]	6.963	1.856	0.007227	4.026	6.716	11.3	10001	100000
	lambda[5]	12.42	2.741	0.009095	7.806	12.14	18.55	10001	100000
	lambda[6]	10.11	2.421	0.01008	5.832	9.943	15.26	10001	100000
	lambda[7]	4.343	1.247	0.004056	2.348	4.186	7.238	10001	100000
	lambda[8]	8.165	2.033	0.006821	4.659	8.013	12.61	10001	100000
	lambda[9]	7.135	1.821	0.005736	4.059	6.96	11.2	10001	100000
	lambda[10]	6.093	1.624	0.005321	3.446	5.904	9.813	10001	100000
	lambda[11]	7.784	1.99	0.007999	4.331	7.636	12.1	10001	100000
	lambda[12]	7.145	1.834	0.005974	4.133	6.939	11.34	10001	100000
	lambda[13]	22.71	4.463	0.02371	15.06	22.36	32.43	10001	100000
	lambda[14]	5.626	1.545	0.006454	2.972	5.51	9.018	10001	100000
	lambda[15]	9.723	2.298	0.007391	5.838	9.504	14.87	10001	100000
	lambda[16]	5.882	1.621	0.007211	3.093	5.755	9.407	10001	100000
	lambda[17]	7.997	1.976	0.006478	4.665	7.823	12.42	10001	100000
	lambda[18]	10.51	2.461	0.00916	6.197	10.34	15.82	10001	100000
	lambda[19]	19.37	3.866	0.01661	12.86	19.02	27.93	10001	100000
	lambda[20]	8.463	2.071	0.006559	4.996	8.251	13.11	10001	100000
	lambda[21]	6.927	1.787	0.006438	3.863	6.78	10.85	10001	100000
	lambda[22]	1.687	0.5565	0.001847	0.7917	1.623	2.978	10001	100000
	lambda[23]	9.319	2.227	0.006704	5.513	9.128	14.25	10001	100000
	lambda[24]	2.892	0.9132	0.002954	1.501	2.751	5.086	10001	100000
	lambda[25]	9.857	2.322	0.008437	5.823	9.68	14.93	10001	100000
	lambda[26]	9.284	2.516	0.01476	5.423	8.933	15.16	10001	100000
	lambda[27]	14.33	3.103	0.01228	9.192	14.01	21.35	10001	100000
	lambda[28]	10.45	2.402	0.008028	6.309	10.26	15.73	10001	100000
	lambda[29]	7.167	1.828	0.005601	4.092	6.997	11.26	10001	100000
	lambda[30]	6.883	1.973	0.01155	3.431	6.744	11.1	10001	100000
	lambda[31]	11.11	2.594	0.01089	6.539	10.95	16.64	10001	100000
	lambda[32]	4.597	1.357	0.006398	2.247	4.496	7.542	10001	100000
	tau	11.36	24.98	0.3999	3.444	8.671	32.06	10001	100000

Smoking and CHD 

2a.  Poisson Regression for Rates
       Death rates for smokers and non smokers.  British male doctor study.
model{
  for(i in 1:10){
    y[i] ~ dpois(lambda[i])
    lambda[i] <- M[i]*theta[i]
    log(theta[i]) <- beta[1] + beta[2]*A[i] + beta[3]*S[i] 
    log(c[i]) <- logfact(y[i]) - y[i] *(log(M[i]) + log(theta[i])) + M[i]*theta[i]
  }
  for(i in 1:3){beta[i] ~ dnorm(0,0.01)}
  for(i in 1:5){RateRatio[i] <- theta[i] / theta[i+5]}
}
INITS
list(beta=c(0,0,0))

DATA
list(A=c(1,2,3,4,5,1,2,3,4,5), 
      S=c(1,1,1,1,1,0,0,0,0,0),
       y=c(32,104,206,186,102, 2,12,28,28,31),          
       M=c(52407,43248,28612,12663,5317,18790,10673,5710,2585,1462))

RESULTS
With log(theta[i]) <- beta[1] + beta[2]*A[i] + beta[3]*S[i] 
 	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	RateRatio[1]	1.515	0.1638	2.956E-4	1.225	1.504	1.866	10001	500000
	RateRatio[2]	1.515	0.1638	2.956E-4	1.225	1.504	1.866	10001	500000
	RateRatio[3]	1.515	0.1638	2.956E-4	1.225	1.504	1.866	10001	500000
	RateRatio[4]	1.515	0.1638	2.956E-4	1.225	1.504	1.866	10001	500000
	RateRatio[5]	1.515	0.1638	2.956E-4	1.225	1.504	1.866	10001	500000
	beta[1]	-8.123	0.1394	2.979E-4	-8.399	-8.121	-7.853	10001	500000
	beta[2]	0.8357	0.02905	5.416E-5	0.7788	0.8356	0.8928	10001	500000
	beta[3]	0.4094	0.1074	1.938E-4	0.2028	0.4078	0.6239	10001	500000
	c[1]	35390.0	928200.0	1303.0	156.0	2788.0	170400.0	10001	500000
	c[2]	31.73	13.45	0.01956	25.59	27.7	61.63	10001	500000
	c[3]	931600.0	7.77E+7	109700.0	1237.0	38420.0	3.813E+6	10001	500000
	c[4]	753.9	4779.0	7.015	50.61	246.0	4064.0	10001	500000
	c[5]	9.095E+9	2.645E+12	3.732E+9	3817.0	8.14E+5	1.295E+9	10001	500000
	c[6]	14680.0	66020.0	105.5	481.3	4726.0	85730.0	10001	500000
	c[7]	25.01	21.7	0.03952	9.647	19.15	75.64	10001	500000
	c[8]	64.69	101.9	0.1676	15.48	40.51	258.6	10001	500000
	c[9]	45.45	61.07	0.09464	13.99	30.3	167.7	10001	500000
	c[10]	20.88	16.92	0.0246	14.0	16.29	54.76	10001	500000
	theta[1]	0.001034	8.338E-5	1.448E-7	8.783E-4	0.001031	0.001205	10001	500000
	theta[2]	0.002381	1.362E-4	2.215E-7	0.002122	0.002379	0.002656	10001	500000
	theta[3]	0.005488	2.285E-4	3.266E-7	0.005048	0.005485	0.005947	10001	500000
	theta[4]	0.01266	5.49E-4	8.284E-7	0.01161	0.01265	0.01376	10001	500000
	theta[5]	0.02923	0.001781	3.049E-6	0.02584	0.02918	0.03283	10001	500000
	theta[6]	6.895E-4	8.334E-5	1.757E-7	5.381E-4	6.854E-4	8.643E-4	10001	500000
	theta[7]	0.001588	1.699E-4	3.457E-7	0.001274	0.001581	0.00194	10001	500000
	theta[8]	0.00366	3.659E-4	6.918E-7	0.002979	0.003648	0.004411	10001	500000
	theta[9]	0.008442	8.542E-4	1.46E-6	0.006852	0.008413	0.0102	10001	500000
	theta[10]	0.01949	0.002151	3.374E-6	0.01552	0.0194	0.02396	10001	500000
Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
y		127.252	124.254	2.998	130.251	
total	127.252	124.254	2.998	130.251	
2b.  Poisson Regression for Rates
       Death rates for smokers and non smokers.  British male doctor study.
model{
  for(i in 1:10){
    y[i] ~ dpois(lambda[i])
    lambda[i] <- M[i]*theta[i]
    log(theta[i]) <- beta[1] + beta[2]*A[i] + beta[3]*S[i] + beta[4]*A[i]*S[i]
    log(c[i]) <- logfact(y[i]) - y[i] *(log(M[i]) + log(theta[i])) + M[i]*theta[i]
  }
  for(i in 1:4){beta[i] ~ dnorm(0,0.01)}
  for(i in 1:5){RateRatio[i] <- theta[i] / theta[i+5]}
}

INITS
list(beta=c(0,0,0,0))

DATA
list(A=c(1,2,3,4,5,1,2,3,4,5), 
      S=c(1,1,1,1,1,0,0,0,0,0),
       y=c(32,104,206,186,102, 2,12,28,28,31),          
       M=c(52407,43248,28612,12663,5317,18790,10673,5710,2585,1462))

RESULTS
With  log(theta[i]) <- beta[1] + beta[2]*A[i] + beta[3]*S[i] + beta[4]*A[i]*S[i]
	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	RateRatio[1]	2.924	0.754	0.004828	1.767	2.819	4.694	10001	500000
	RateRatio[2]	2.245	0.4073	0.002432	1.579	2.201	3.169	10001	500000
	RateRatio[3]	1.736	0.2149	9.351E-4	1.364	1.719	2.203	10001	500000
	RateRatio[4]	1.352	0.1521	2.339E-4	1.083	1.342	1.679	10001	500000
	RateRatio[5]	1.06	0.1662	6.393E-4	0.777	1.046	1.428	10001	500000
	beta[1]	-8.876	0.3064	0.001999	-9.498	-8.868	-8.297	10001	500000
	beta[2]	1.047	0.0776	5.035E-4	0.8972	1.046	1.202	10001	500000
	beta[3]	1.29	0.3267	0.002127	0.6688	1.284	1.949	10001	500000
	beta[4]	-0.2488	0.08376	5.394E-4	-0.4148	-0.2479	-0.08667	10001	500000
	c[1]	2.882E+6	4.004E+8	563600.0	585.0	25090.0	4.737E+6	10001	500000
	c[2]	36.97	33.43	0.04916	25.59	29.2	93.53	10001	500000
	c[3]	321300.0	9.304E+6	12990.0	768.2	19290.0	1.533E+6	10001	500000
	c[4]	2082.0	19670.0	27.7	66.59	465.8	11900.0	10001	500000
	c[5]	1.124E+8	1.51E+10	2.122E+7	448.0	42690.0	3.712E+7	10001	500000
	c[6]	303.5	4007.0	9.282	9.596	65.52	1648.0	10001	500000
	c[7]	10.75	4.126	0.01516	8.745	9.445	20.51	10001	500000
	c[8]	306.3	1415.0	3.596	22.49	109.7	1663.0	10001	500000
	c[9]	24.81	21.64	0.03211	13.33	18.83	72.22	10001	500000
	c[10]	869.4	302400.0	427.4	14.04	30.27	941.5	10001	500000
	theta[1]	0.001131	9.529E-5	1.792E-7	9.54E-4	0.001128	0.001327	10001	500000
	theta[2]	0.002509	1.465E-4	2.629E-7	0.00223	0.002506	0.002804	10001	500000
	theta[3]	0.005568	2.308E-4	3.535E-7	0.005125	0.005565	0.006029	10001	500000
	theta[4]	0.01237	5.546E-4	8.045E-7	0.01131	0.01236	0.01348	10001	500000
	theta[5]	0.02751	0.001802	2.963E-6	0.0241	0.02747	0.03116	10001	500000
	theta[6]	4.088E-4	9.528E-5	5.982E-7	2.467E-4	4.005E-4	6.184E-4	10001	500000
	theta[7]	0.001149	1.911E-4	1.124E-6	8.049E-4	0.001139	0.001551	10001	500000
	theta[8]	0.00325	3.723E-4	1.621E-6	0.002558	0.003238	0.004014	10001	500000
	theta[9]	0.009247	9.425E-4	1.417E-6	0.00749	0.009216	0.01118	10001	500000
	theta[10]	0.02646	0.003701	1.44E-5	0.01974	0.02628	0.03425	10001	500000

	Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
y		118.971	114.971	4.000	122.971	
total	118.971	114.971	4.000	122.971	
2c.  Poisson Regression for Rates
       Death rates for smokers and non smokers.  British male doctor study.

model{
  for(i in 1:10){
    y[i] ~ dpois(lambda[i])
    lambda[i] <- M[i]*theta[i]
    log(theta[i]) <- beta[1] + beta[2]*A[i] + beta[3]*A[i]*A[i] + beta[4]*S[i] + beta[5]*A[i]*S[i] + beta[6]*A[i]*A[i]*S[i]
    log(c[i]) <- logfact(y[i]) - y[i] *(log(M[i]) + log(theta[i])) + M[i]*theta[i]
    Rate[i] <- 10000*theta[i]
  }
  for(i in 1:6){beta[i] ~ dnorm(0,0.01)}
  for(i in 1:5){
    RateRatio[i] <- theta[i] / theta[i+5]
    RateDifference[i] <- 10000*(theta[i] - theta[i+5])
  }
}

INITS
list(beta=c(0,0,0,0,0,0))

DATA
list(A=c(1,2,3,4,5,1,2,3,4,5), 
  S=c(1,1,1,1,1,0,0,0,0,0),
  y=c(32,104,206,186,102, 2,12,28,28,31),          
  M=c(52407,43248,28612,12663,5317,18790,10673,5710,2585,1462))

RESULTS

With  log(theta[i]) <- beta[1] + beta[2]*A[i] + beta[3]*A[i]*A[i] + beta[4]*S[i] + beta[5]*A[i]*S[i] + beta[6]*A[i]*A[i]*S[i]
	node	 mean	 sd	 MC error	             2.5%	median	97.5%	start	sample
	RateRatio[1]	4.623	2.54	0.06538	      1.702	4.014	      11.03	50001	500000
	RateRatio[2]	2.479	0.5433	0.007858	1.633	2.405	       3.747	50001	500000
	RateRatio[3]	1.621	0.2533	0.003387	1.192	1.598	       2.182	50001	500000
	RateRatio[4]	1.195	0.1629	0.002654	0.9119   1.183	       1.549	50001	500000
	RateRatio[5]	0.9987	0.2032	0.002285   0.6692   0.9755	     1.461	50001	500000
	beta[1]	-11.35	0.8968	0.02503	-13.18	-11.31	-9.709	50001	500000
	beta[2]	2.729	0.5375	0.01546	1.721	2.713	3.804	50001	500000
	beta[3]	-0.2483	0.07655	0.002171	-0.4002	-0.2465	-0.103	50001	500000
	beta[4]	2.049	0.9415	0.02638	0.3049	2.013	3.966	50001	500000
	beta[5]	-0.6937	0.5683	0.01636	-1.834	-0.6792	0.3821	50001	500000
	beta[6]	0.05595	0.08155	0.002305	-0.1002	0.05425	0.2181	50001	500000
	c[1]	27.2	72.03	0.206	14.22	17.0	95.5	50001	500000
	c[2]	36.34	31.33	0.05093	25.59	28.79	91.39	50001	500000
	c[3]	55.76	83.2	0.1468	36.0	41.12	156.9	50001	500000
	c[4]	50.24	54.52	0.1042	34.21	38.81	132.3	50001	500000
	c[5]	64.93	1055.0	1.507	25.35	30.85	227.9	50001	500000
	c[6]	6.958	38.5	0.2156	3.696	4.339	23.21	50001	500000
	c[7]	13.03	14.69	0.1406	8.746	9.918	35.23	50001	500000
	c[8]	24.03	66.1	0.3018	13.31	15.83	78.96	50001	500000
	c[9]	27.26	76.15	0.5181	13.31	16.64	97.61	50001	500000
	c[10]	59.11	3169.0	4.673	14.0	17.75	186.9	50001	500000
	theta[1]	5.805E-4	8.518E-5	6.022E-7	4.268E-4	5.758E-4	7.602E-4	50001	500000
	theta[2]	0.002472	1.59E-4	4.776E-7	0.00217	0.002469	0.002794	50001	500000
	theta[3]	0.007225	3.845E-4	1.805E-6	0.006493	0.007218	0.007999	50001	500000
	theta[4]	0.01438	7.165E-4	2.687E-6	0.01301	0.01436	0.01582	50001	500000
	theta[5]	0.01953	0.00177	7.314E-6	0.01623	0.01948	0.02317	50001	500000
	theta[6]	1.549E-4	6.969E-5	1.708E-6	5.407E-5	1.434E-4	3.22E-4	50001	500000
	theta[7]	0.001039	2.069E-4	2.885E-6	6.697E-4	0.001027	0.001478	50001	500000
	theta[8]	0.004552	6.564E-4	8.74E-6	0.003366	0.004518	0.005934	50001	500000
	theta[9]	0.01222	0.001536	2.595E-5	0.009435	0.01215	0.01545	50001	500000
	theta[10]	0.02018	0.003553	4.031E-5	0.01384	0.01997	0.02773	50001	500000

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
		Dbar		Dhat		pD		DIC	
y		62.271	56.346	5.925	68.196	
total	62.271	56.346	5.925	68.196	

 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	Rate[1]	5.8	0.8546	0.004743	4.259	5.752	7.604	100001	900000
	Rate[2]	24.72	1.592	0.003608	21.69	24.69	27.94	100001	900000
	Rate[3]	72.26	3.856	0.01436	64.92	72.19	80.04	100001	900000
	Rate[4]	143.8	7.168	0.02097	130.1	143.7	158.2	100001	900000
	Rate[5]	195.3	17.7	0.05866	162.2	194.7	231.6	100001	900000
	Rate[6]	1.576	0.7128	0.0144	0.558	1.452	3.294	100001	900000
	Rate[7]	10.43	2.07	0.02308	6.749	10.31	14.84	100001	900000
	Rate[8]	45.4	6.558	0.07124	33.54	45.07	59.2	100001	900000
	Rate[9]	121.9	15.34	0.2109	93.86	121.2	153.9	100001	900000
	Rate[10]	202.4	35.63	0.3272	138.9	200.2	278.2	100001	900000
	RateDifference[1]	4.223	1.109	0.0174	1.976	4.24	6.362	100001	900000
	RateDifference[2]	14.29	2.61	0.02505	8.967	14.36	19.21	100001	900000
	RateDifference[3]	26.86	7.605	0.07986	11.25	27.1	41.11	100001	900000
	RateDifference[4]	21.94	16.92	0.2227	-12.81	22.49	53.49	100001	900000
	RateDifference[5]	-7.112	39.79	0.3611	-89.78	-5.494	66.23	100001	900000
	RateRatio[1]	4.526	2.405	0.04897	1.659	3.961	10.66	100001	900000
	RateRatio[2]	2.467	0.5376	0.006114	1.626	2.395	3.718	100001	900000
	RateRatio[3]	1.625	0.2544	0.002796	1.194	1.602	2.189	100001	900000
	RateRatio[4]	1.199	0.1641	0.002197	0.9153	1.186	1.557	100001	900000
	RateRatio[5]	0.9956	0.2024	0.001822	0.6671	0.9725	1.457	100001	900000
	beta[1]	-11.32	0.8932	0.02026	-13.12	-11.29	-9.652	100001	900000
	beta[2]	2.706	0.5362	0.01255	1.682	2.7	3.765	100001	900000
	beta[3]	-0.2451	0.07643	0.001764	-0.3946	-0.2448	-0.09736	100001	900000
	beta[4]	2.009	0.9401	0.02134	0.2419	1.989	3.893	100001	900000
	beta[5]	-0.6693	0.5682	0.01327	-1.789	-0.6662	0.4218	100001	900000
	beta[6]	0.05253	0.08158	0.001871	-0.106	0.05244	0.2121	100001	900000

3. Example 11.2.1.   Poisson regression for longitudinal data:  Armadillo hunting

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