Chapter 9: Linear Regression
FEV data from Rosner; 5th edition;  p. 479
Use only the subset of data from individuals aged 10-19.  Age 10 was the youngest smoker
n=345
y=FEV= Forced expiratory volume = a measure of lung function; it ranges from about 2 to 5 in the data set;  higher is better
x=(Age, Smoke)
Smoke is binary (yes/no)
Ages range from 10 to 19 

Case 1.  Linear regression.  Need to compute C0inv and betastar in R

Model

model{
for(i in 1:n){
	FEV[i] ~ dnorm(mu[i],tau)
	mu[i] <- beta[1] + beta[2]*Age[i] + beta[3]*Smoke[i]  + beta[4]*Age[i]*Smoke[i]
    }
beta[1:r] ~ dmnorm(betastar[1:r],C0inv[1:r,1:r])   # BCJ prior
tau ~ dgamma(a,b)
#beta[1] ~ dnorm(0,0.001)  # Diffuse prior
#beta[2] ~ dnorm(0,0.001)
#beta[3] ~ dnorm(0,0.001)
#beta[4] ~ dnorm(0,0.001)
#tau ~ dgamma(0.001,0.001)

## Estimate mean FEV for smokers and nonsmokers who are 10, ..., 19 years old
for(i in 1:10){
   meanFEVs[i] <-  beta[1] + (beta[2]+beta[4])*(i+9) + beta[3]
   meanFEVns[i] <-  beta[1] + beta[2]*(i+9)
  }

## Easy to estimate relative means and mean differences as well
RM <- meanFEVns[9]  / meanFEVs[9]    ## RM comparing 18 year old smoker to 18 year old nonsmoker
MD <- meanFEVs[9]  - meanFEVns[4]   ## MD comparing 18 year old smoker to 13 year old nonsmoker

## Predict the FEV for a 20 year old smoker and nonsmoker
FEV20s ~ dnorm(mu20s,tau)
FEV20ns ~ dnorm(mu20ns,tau)
mu20s <-  beta[1] + (beta[2]+beta[4])*20 + beta[3]
mu20ns <-  beta[1] + beta[2]*20 
}

DATA  

list(n=345, r=4, betastar=c(0.16,0.24,2.06,-0.18),a = 1.73,b = 0.78)  # Load for BCJ prior
list(n=345)  # Load for diffuse prior

C0inv[,1] C0inv[,2] C0inv[,3]  C0inv[,4]    # Load for BCJ prior
86.11111  1200  36.11111  525
1200.00000 17250 525.00000 7825
36.11111   525  36.11111  525
525.00000  7825 525.00000 7825
END

Age[]	FEV[]	Smoke[]
19	5.102	0
19	3.519	1
19	3.345	1
18	2.906	0
18	3.082	0
18	4.22	0
18	4.086	1
18	4.404	1
18	2.853	0
17	4.429	0
17	3.082	1
17	4.724	0
17	3.406	1
17	3.5	0
17	5.633	0
17	3.96	0
17	5.638	0
16	3.688	1
16	3.387	0
16	2.903	1
16	3.674	0
16	2.608	1
16	3.645	0
16	4.07	1
16	4.299	0
16	2.981	0
16	4.504	0
16	4.872	1
16	4.27	1
16	2.795	1
15	4.284	0
15	4.506	1
15	4.279	0
15	4.5	0
15	2.635	0
15	2.679	1
15	2.198	1
15	3.004	1
15	5.793	0
15	3.985	0
15	3.731	0
15	3.122	1
15	3.33	1
15	3.799	1
15	2.887	0
15	2.264	1
15	2.278	1
15	3.727	1
15	3.211	0
14	3.381	0
14	3.741	0
14	4.842	0
14	2.236	1
14	3.436	0
14	3.428	1
14	4.683	0
14	3.957	1
14	3.074	1
14	4.273	0
14	3.395	0
14	3.169	0
14	4.309	1
14	4.111	0
14	2.538	0
14	3.68	0
14	3.585	0
14	2.891	0
14	3.78	0
14	4.271	0
14	2.934	0
14	2.276	1
14	2.997	0
14	3.806	0
14	4.763	1
13	4.336	0
13	3.816	0
13	3.549	0
13	3.147	0
13	4.305	0
13	3.906	0
13	2.819	0
13	3.994	0
13	3.208	1
13	3.193	0
13	4.789	1
13	3.745	0
13	3.887	0
13	2.646	0
13	2.976	0
13	2.599	1
13	3.06	0
13	4.756	1
13	2.704	0
13	2.434	0
13	3.086	1
13	4.225	0
13	3.089	0
13	4.877	0
13	2.216	1
13	3.078	1
13	3.297	1
13	4.448	0
13	3.984	0
13	2.677	1
13	3.255	0
13	3.297	1
13	3.331	0
13	5.083	0
13	3.056	0
13	3.785	1
13	2.449	0
13	4.232	0
13	2.531	0
13	3.577	0
13	3.152	1
13	3.141	0
13	4.045	1
12	3.058	0
12	4.55	0
12	2.841	0
12	2.889	0
12	4.393	0
12	2.384	1
12	2.971	0
12	3.835	1
12	2.751	0
12	2.556	0
12	2.868	0
12	3.255	0
12	1.916	0
12	2.241	0
12	5.224	0
12	4.08	0
12	4.411	0
12	3.791	0
12	3.343	1
12	2.913	0
12	3.279	0
12	2.347	0
12	3.751	1
12	3.079	0
12	3.403	0
12	3.501	0
12	3.186	1
12	4.073	0
12	2.752	0
12	2.304	1
12	3.692	0
12	3.082	0
12	4.72	0
12	2.417	0
12	2.759	1
12	3.231	0
12	3.529	0
12	2.866	0
12	2.605	0
12	2.569	0
12	3.924	0
12	2.752	0
12	3.53	0
12	2.332	0
12	3.035	0
12	4.831	0
12	2.714	0
12	3.519	0
12	3.341	0
12	2.822	0
12	2.935	0
12	2.499	0
12	3.001	0
12	3.222	0
12	4.203	0
12	3.082	0
12	2.579	0
11	2.884	0
11	2.17	0
11	3.47	0
11	2.524	0
11	3.654	0
11	2.665	0
11	3.222	0
11	3.49	0
11	2.957	0
11	4.007	0
11	2.386	0
11	2.762	0
11	3.011	0
11	3.583	0
11	3.236	0
11	3.058	0
11	3.108	0
11	1.694	1
11	2.346	0
11	3.515	0
11	2.754	0
11	2.463	0
11	2.633	0
11	3.111	0
11	3.977	0
11	3.411	0
11	3.171	0
11	3.587	0
11	3.845	0
11	2.417	0
11	2.735	0
11	4.065	0
11	2.318	0
11	3.596	0
11	2.542	0
11	2.354	0
11	2.491	0
11	2.993	0
11	3.774	0
11	2.988	0
11	2.498	0
11	2.887	0
11	3.515	0
11	3.425	0
11	4.593	0
11	3.223	0
11	4.073	0
11	2.606	0
11	3.169	1
11	2.465	0
11	2.827	0
11	2.578	0
11	2.081	0
11	2.974	0
11	3.102	1
11	3.023	0
11	3.681	0
11	3.258	0
11	3.247	0
11	4.324	0
11	2.362	0
11	2.563	0
11	3.206	0
11	2.953	1
11	3.078	0
11	3.369	0
11	3.022	0
11	2.866	0
11	2.501	0
11	3.32	0
11	2.123	0
11	3.847	0
11	2.928	0
11	2.689	0
11	2.894	0
11	4.637	1
11	2.812	0
11	2.387	0
11	4.13	0
11	3.28	0
11	2.659	0
11	2.822	0
11	2.14	0
11	3.12	1
11	2.562	0
11	3.339	1
11	2.458	0
11	3.104	1
11	3.069	1
11	2.785	0
10	2.328	0
10	1.811	0
10	2.642	0
10	3.166	0
10	2.561	0
10	2.481	0
10	3.203	0
10	3.111	0
10	2.52	0
10	2.292	0
10	2.246	0
10	1.937	0
10	2.646	0
10	3.251	0
10	3.007	0
10	3.489	0
10	2.864	0
10	2.25	0
10	2.352	0
10	2.592	0
10	2.72	0
10	3.048	0
10	2.094	0
10	3.183	0
10	3.354	0
10	2.387	1
10	2.504	0
10	2.891	0
10	1.823	0
10	2.175	0
10	2.673	0
10	2.608	0
10	1.458	0
10	3.795	0
10	2.545	0
10	3.305	0
10	2.855	0
10	2.287	0
10	2.365	0
10	2.696	0
10	2.813	0
10	3.413	1
10	1.858	0
10	2.975	1
10	3.35	0
10	2.901	0
10	3.2	0
10	2.358	0
10	2.581	0
10	2.691	0
10	1.873	0
10	2.758	0
10	3.05	0
10	2.201	0
10	1.858	0
10	1.665	0
10	2.25	0
10	2.862	0
10	2.356	0
10	4.591	0
10	2.216	0
10	3.498	1
10	2.364	0
10	2.341	0
10	3.127	0
10	2.132	0
10	3.456	0
10	3.073	0
10	2.688	0
10	3.329	0
10	3.11	0
10	2.435	0
10	2.838	0
10	3.086	0
10	2.77	0
10	3.09	0
10	3.038	1
10	2.568	0
10	3.132	0
10	2.391	0
10	2.1	0
END

INITS
list(tau=1, beta=c(0,0,0,0), FEV20s=2, FEV20ns=2)

RESULTS using BCJ prior

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	FEV20ns	4.864	0.683	0.003097	3.529	4.866	6.197	10000	50001
	FEV20s	3.611	0.6978	0.002921	2.229	3.611	4.977	10000	50001
	MD	0.1057	0.1612	7.758E-4	-0.2111	0.1059	0.4227	10000	50001
	RM	1.273	0.06682	3.06E-4	1.148	1.271	1.413	10000	50001
	beta[1]	0.6599	0.2327	0.001027	0.2042	0.6602	1.113	10000	50001
	beta[2]	0.2103	0.01923	8.594E-5	0.173	0.2103	0.248	10000	50001
	beta[3]	1.8	0.5077	0.002177	0.8072	1.803	2.787	10000	50001
	beta[4]	-0.1526	0.03757	1.63E-4	-0.2257	-0.1528	-0.07853	10000	50001
	meanFEVns[1]	2.763	0.05326	2.337E-4	2.659	2.763	2.867	10000	50001
	meanFEVns[2]	2.974	0.04212	1.893E-4	2.89	2.974	3.056	10000	50001
	meanFEVns[3]	3.184	0.03809	1.783E-4	3.109	3.184	3.259	10000	50001
	meanFEVns[4]	3.394	0.04321	2.063E-4	3.31	3.394	3.48	10000	50001
	meanFEVns[5]	3.605	0.05499	2.609E-4	3.497	3.605	3.713	10000	50001
	meanFEVns[6]	3.815	0.07014	3.291E-4	3.678	3.815	3.952	10000	50001
	meanFEVns[7]	4.025	0.08692	4.042E-4	3.856	4.025	4.196	10000	50001
	meanFEVns[8]	4.236	0.1045	4.829E-4	4.032	4.236	4.44	10000	50001
	meanFEVns[9]	4.446	0.1227	5.637E-4	4.207	4.446	4.687	10000	50001
	meanFEVns[10]	4.656	0.1411	6.459E-4	4.382	4.656	4.932	10000	50001
	meanFEVs[1]	3.038	0.1433	6.265E-4	2.755	3.038	3.318	10000	50001
	meanFEVs[2]	3.096	0.1169	5.05E-4	2.864	3.096	3.323	10000	50001
	meanFEVs[3]	3.153	0.09439	4.03E-4	2.968	3.154	3.338	10000	50001
	meanFEVs[4]	3.211	0.07893	3.383E-4	3.056	3.212	3.366	10000	50001
	meanFEVs[5]	3.269	0.07508	3.335E-4	3.122	3.269	3.416	10000	50001
	meanFEVs[6]	3.327	0.08444	3.909E-4	3.161	3.327	3.491	10000	50001
	meanFEVs[7]	3.384	0.1035	4.89E-4	3.181	3.384	3.587	10000	50001
	meanFEVs[8]	3.442	0.128	6.084E-4	3.191	3.442	3.693	10000	50001
	meanFEVs[9]	3.5	0.1553	7.389E-4	3.193	3.5	3.806	10000	50001
	meanFEVs[10]	3.558	0.1843	8.756E-4	3.192	3.558	3.922	10000	50001
	tau	2.269	0.1733	8.017E-4	1.94	2.265	2.62	10000	50001

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
FEV	697.164	692.817	4.347	701.510	
total	697.164	692.817	4.347	701.510	

RESULTS using diffuse prior

node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	FEV20ns	4.869	0.6876	0.003126	3.525	4.871	6.209	10000	50001
	FEV20s	3.72	0.7094	0.002973	2.316	3.721	5.103	10000	50001
	MD	0.1897	0.188	8.951E-4	-0.1796	0.1898	0.5601	10000	50001
	RM	1.243	0.07402	3.376E-4	1.106	1.24	1.398	10000	50001
	beta[1]	0.6727	0.2519	0.001117	0.1792	0.673	1.163	10000	50001
	beta[2]	0.2099	0.02106	9.471E-5	0.1691	0.2099	0.2511	10000	50001
	beta[3]	1.718	0.5655	0.002416	0.6095	1.721	2.817	10000	50001
	beta[4]	-0.1432	0.04233	1.832E-4	-0.2256	-0.1434	-0.05966	10000	50001
	meanFEVns[1]	2.772	0.05496	2.397E-4	2.664	2.772	2.879	10000	50001
	meanFEVns[2]	2.982	0.04308	1.924E-4	2.897	2.982	3.066	10000	50001
	meanFEVns[3]	3.192	0.03974	1.858E-4	3.113	3.192	3.27	10000	50001
	meanFEVns[4]	3.402	0.0468	2.235E-4	3.31	3.402	3.494	10000	50001
	meanFEVns[5]	3.612	0.06073	2.887E-4	3.493	3.611	3.731	10000	50001
	meanFEVns[6]	3.821	0.07793	3.67E-4	3.669	3.821	3.974	10000	50001
	meanFEVns[7]	4.031	0.09668	4.516E-4	3.843	4.031	4.221	10000	50001
	meanFEVns[8]	4.241	0.1162	5.396E-4	4.015	4.241	4.469	10000	50001
	meanFEVns[9]	4.451	0.1362	6.295E-4	4.186	4.451	4.718	10000	50001
	meanFEVns[10]	4.661	0.1565	7.207E-4	4.357	4.661	4.967	10000	50001
	meanFEVs[1]	3.058	0.1563	6.772E-4	2.749	3.058	3.365	10000	50001
	meanFEVs[2]	3.124	0.1268	5.432E-4	2.874	3.125	3.372	10000	50001
	meanFEVs[3]	3.191	0.1021	4.337E-4	2.99	3.191	3.391	10000	50001
	meanFEVs[4]	3.258	0.08637	3.713E-4	3.089	3.258	3.426	10000	50001
	meanFEVs[5]	3.324	0.08486	3.798E-4	3.159	3.325	3.491	10000	50001
	meanFEVs[6]	3.391	0.09821	4.553E-4	3.199	3.391	3.583	10000	50001
	meanFEVs[7]	3.458	0.1216	5.719E-4	3.218	3.457	3.698	10000	50001
	meanFEVs[8]	3.525	0.1505	7.096E-4	3.229	3.524	3.821	10000	50001
	meanFEVs[9]	3.591	0.1822	8.582E-4	3.231	3.591	3.952	10000	50001
	meanFEVs[10]	3.658	0.2155	0.001013	3.231	3.658	4.084	10000	50001
	tau	2.268	0.1744	8.091E-4	1.937	2.264	2.621	10000	50001

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
FEV	697.348	692.330	5.018	702.366	
total	697.348	692.330	5.018	702.366

Case 2.  Same model as case 1 above, slightly different coding.  Need to compute Xtildeinv in R

Model

model{
for(i in 1:n){
	FEV[i] ~ dnorm(mu[i],tau)
	mu[i] <- beta[1] + beta[2]*Age[i] + beta[3]*Smoke[i]  + beta[4]*Age[i]*Smoke[i]
    }
### BCJ prior
#mtilde[1] ~ dnorm(2.8, 25)
#mtilde[2] ~ dnorm(3.0,25)
#mtilde[3] ~ dnorm(4.0,25)
#mtilde[4] ~ dnorm(3.3,11.11)
#for(j in 1:r){
#    beta[j] <- inprod(Xtildeinv[j,1:r] , mtilde[1:r])
#}
#tau ~ dgamma(a,b)
beta[1] ~ dnorm(0,0.001)  # Diffuse prior
beta[2] ~ dnorm(0,0.001)
beta[3] ~ dnorm(0,0.001)
beta[4] ~ dnorm(0,0.001)
tau ~ dgamma(0.001,0.001)

## Estimate mean FEV for smokers and nonsmokers who are 10, ..., 19 years old
for(i in 1:10){
   meanFEVs[i] <-  beta[1] + (beta[2]+beta[4])*(i+9) + beta[3]
   meanFEVns[i] <-  beta[1] + beta[2]*(i+9)
  }

## Easy to estimate relative means and mean differences as well
RM <- meanFEVns[9]  / meanFEVs[9]    ## RM comparing 18 year old smoker to 18 year old non-smoker
MD <- meanFEVs[9]  - meanFEVns[4]   ## MD comparing 18 year old smoker to 13 year old nonsmoker

## Predict the FEV for a 20 year old smoker and for a 20 year old nonsmoker
FEV20s ~ dnorm(mu20s,tau)
FEV20ns ~ dnorm(mu20ns,tau)
mu20s <-  beta[1] + (beta[2]+beta[4])*20 + beta[3]
mu20ns <-  beta[1] + beta[2]*20 
}

DATA

list(n=345, r=4, a = 1.73,b = 0.78)   # Load for BCJ prior
Xtildeinv[,1] Xtildeinv[,2] Xtildeinv[,3]  Xtildeinv[,4]   # Load for BCJ prior
3.2     0   -2.2	  0
-0.2   0     0.2      0
-3.2  3.6  2.2   -2.6
0.2  -0.2  -0.2  0.2
END
list(n=345) # Load for diffuse prior

Age[]	FEV[]	Smoke[]
19	5.102	0
19	3.519	1
19	3.345	1
18	2.906	0
18	3.082	0
18	4.22	0
18	4.086	1
18	4.404	1
18	2.853	0
17	4.429	0
17	3.082	1
17	4.724	0
17	3.406	1
17	3.5	0
17	5.633	0
17	3.96	0
17	5.638	0
16	3.688	1
16	3.387	0
16	2.903	1
16	3.674	0
16	2.608	1
16	3.645	0
16	4.07	1
16	4.299	0
16	2.981	0
16	4.504	0
16	4.872	1
16	4.27	1
16	2.795	1
15	4.284	0
15	4.506	1
15	4.279	0
15	4.5	0
15	2.635	0
15	2.679	1
15	2.198	1
15	3.004	1
15	5.793	0
15	3.985	0
15	3.731	0
15	3.122	1
15	3.33	1
15	3.799	1
15	2.887	0
15	2.264	1
15	2.278	1
15	3.727	1
15	3.211	0
14	3.381	0
14	3.741	0
14	4.842	0
14	2.236	1
14	3.436	0
14	3.428	1
14	4.683	0
14	3.957	1
14	3.074	1
14	4.273	0
14	3.395	0
14	3.169	0
14	4.309	1
14	4.111	0
14	2.538	0
14	3.68	0
14	3.585	0
14	2.891	0
14	3.78	0
14	4.271	0
14	2.934	0
14	2.276	1
14	2.997	0
14	3.806	0
14	4.763	1
13	4.336	0
13	3.816	0
13	3.549	0
13	3.147	0
13	4.305	0
13	3.906	0
13	2.819	0
13	3.994	0
13	3.208	1
13	3.193	0
13	4.789	1
13	3.745	0
13	3.887	0
13	2.646	0
13	2.976	0
13	2.599	1
13	3.06	0
13	4.756	1
13	2.704	0
13	2.434	0
13	3.086	1
13	4.225	0
13	3.089	0
13	4.877	0
13	2.216	1
13	3.078	1
13	3.297	1
13	4.448	0
13	3.984	0
13	2.677	1
13	3.255	0
13	3.297	1
13	3.331	0
13	5.083	0
13	3.056	0
13	3.785	1
13	2.449	0
13	4.232	0
13	2.531	0
13	3.577	0
13	3.152	1
13	3.141	0
13	4.045	1
12	3.058	0
12	4.55	0
12	2.841	0
12	2.889	0
12	4.393	0
12	2.384	1
12	2.971	0
12	3.835	1
12	2.751	0
12	2.556	0
12	2.868	0
12	3.255	0
12	1.916	0
12	2.241	0
12	5.224	0
12	4.08	0
12	4.411	0
12	3.791	0
12	3.343	1
12	2.913	0
12	3.279	0
12	2.347	0
12	3.751	1
12	3.079	0
12	3.403	0
12	3.501	0
12	3.186	1
12	4.073	0
12	2.752	0
12	2.304	1
12	3.692	0
12	3.082	0
12	4.72	0
12	2.417	0
12	2.759	1
12	3.231	0
12	3.529	0
12	2.866	0
12	2.605	0
12	2.569	0
12	3.924	0
12	2.752	0
12	3.53	0
12	2.332	0
12	3.035	0
12	4.831	0
12	2.714	0
12	3.519	0
12	3.341	0
12	2.822	0
12	2.935	0
12	2.499	0
12	3.001	0
12	3.222	0
12	4.203	0
12	3.082	0
12	2.579	0
11	2.884	0
11	2.17	0
11	3.47	0
11	2.524	0
11	3.654	0
11	2.665	0
11	3.222	0
11	3.49	0
11	2.957	0
11	4.007	0
11	2.386	0
11	2.762	0
11	3.011	0
11	3.583	0
11	3.236	0
11	3.058	0
11	3.108	0
11	1.694	1
11	2.346	0
11	3.515	0
11	2.754	0
11	2.463	0
11	2.633	0
11	3.111	0
11	3.977	0
11	3.411	0
11	3.171	0
11	3.587	0
11	3.845	0
11	2.417	0
11	2.735	0
11	4.065	0
11	2.318	0
11	3.596	0
11	2.542	0
11	2.354	0
11	2.491	0
11	2.993	0
11	3.774	0
11	2.988	0
11	2.498	0
11	2.887	0
11	3.515	0
11	3.425	0
11	4.593	0
11	3.223	0
11	4.073	0
11	2.606	0
11	3.169	1
11	2.465	0
11	2.827	0
11	2.578	0
11	2.081	0
11	2.974	0
11	3.102	1
11	3.023	0
11	3.681	0
11	3.258	0
11	3.247	0
11	4.324	0
11	2.362	0
11	2.563	0
11	3.206	0
11	2.953	1
11	3.078	0
11	3.369	0
11	3.022	0
11	2.866	0
11	2.501	0
11	3.32	0
11	2.123	0
11	3.847	0
11	2.928	0
11	2.689	0
11	2.894	0
11	4.637	1
11	2.812	0
11	2.387	0
11	4.13	0
11	3.28	0
11	2.659	0
11	2.822	0
11	2.14	0
11	3.12	1
11	2.562	0
11	3.339	1
11	2.458	0
11	3.104	1
11	3.069	1
11	2.785	0
10	2.328	0
10	1.811	0
10	2.642	0
10	3.166	0
10	2.561	0
10	2.481	0
10	3.203	0
10	3.111	0
10	2.52	0
10	2.292	0
10	2.246	0
10	1.937	0
10	2.646	0
10	3.251	0
10	3.007	0
10	3.489	0
10	2.864	0
10	2.25	0
10	2.352	0
10	2.592	0
10	2.72	0
10	3.048	0
10	2.094	0
10	3.183	0
10	3.354	0
10	2.387	1
10	2.504	0
10	2.891	0
10	1.823	0
10	2.175	0
10	2.673	0
10	2.608	0
10	1.458	0
10	3.795	0
10	2.545	0
10	3.305	0
10	2.855	0
10	2.287	0
10	2.365	0
10	2.696	0
10	2.813	0
10	3.413	1
10	1.858	0
10	2.975	1
10	3.35	0
10	2.901	0
10	3.2	0
10	2.358	0
10	2.581	0
10	2.691	0
10	1.873	0
10	2.758	0
10	3.05	0
10	2.201	0
10	1.858	0
10	1.665	0
10	2.25	0
10	2.862	0
10	2.356	0
10	4.591	0
10	2.216	0
10	3.498	1
10	2.364	0
10	2.341	0
10	3.127	0
10	2.132	0
10	3.456	0
10	3.073	0
10	2.688	0
10	3.329	0
10	3.11	0
10	2.435	0
10	2.838	0
10	3.086	0
10	2.77	0
10	3.09	0
10	3.038	1
10	2.568	0
10	3.132	0
10	2.391	0
10	2.1	0
END

INITS
list(tau=1,beta=c(0,0,0,0), FEV20s=2, FEV20ns=2)  # Load for diffuse prior
list(tau=1,mtilde=c(0,0,0,0), FEV20s=2, FEV20ns=2)  # Load for BCJ prior

RESULTS using BCJ prior

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	FEV20ns	4.861	0.6816	0.003166	3.52	4.862	6.197	10000	50001
	FEV20s	3.61	0.6982	0.002906	2.222	3.611	4.969	10000	50001
	MD	0.1056	0.161	6.962E-4	-0.2082	0.105	0.4215	10000	50001
	RM	1.272	0.06691	2.844E-4	1.149	1.27	1.41	10000	50001
	beta[1]	0.6638	0.2327	0.001107	0.2063	0.6641	1.12	10000	50001
	beta[2]	0.21	0.01927	9.122E-5	0.1722	0.21	0.2478	10000	50001
	beta[3]	1.802	0.5074	0.002148	0.8067	1.803	2.793	10000	50001
	beta[4]	-0.1526	0.0376	1.579E-4	-0.2255	-0.1527	-0.07888	10000	50001
	meanFEVns[1]	2.764	0.05305	2.445E-4	2.66	2.764	2.868	10000	50001
	meanFEVns[2]	2.974	0.04204	1.862E-4	2.891	2.974	3.056	10000	50001
	meanFEVns[3]	3.184	0.03825	1.618E-4	3.109	3.184	3.259	10000	50001
	meanFEVns[4]	3.394	0.04361	1.853E-4	3.309	3.394	3.479	10000	50001
	meanFEVns[5]	3.604	0.05552	2.432E-4	3.495	3.604	3.712	10000	50001
	meanFEVns[6]	3.814	0.07076	3.172E-4	3.675	3.814	3.952	10000	50001
	meanFEVns[7]	4.024	0.0876	3.983E-4	3.853	4.024	4.195	10000	50001
	meanFEVns[8]	4.234	0.1053	4.831E-4	4.029	4.234	4.44	10000	50001
	meanFEVns[9]	4.444	0.1234	5.699E-4	4.204	4.443	4.687	10000	50001
	meanFEVns[10]	4.654	0.1419	6.579E-4	4.377	4.653	4.933	10000	50001
	meanFEVs[1]	3.04	0.1426	6.397E-4	2.761	3.041	3.32	10000	50001
	meanFEVs[2]	3.097	0.1164	5.346E-4	2.871	3.098	3.326	10000	50001
	meanFEVs[3]	3.155	0.09396	4.463E-4	2.971	3.155	3.34	10000	50001
	meanFEVs[4]	3.212	0.07861	3.866E-4	3.059	3.213	3.367	10000	50001
	meanFEVs[5]	3.27	0.07488	3.695E-4	3.123	3.27	3.418	10000	50001
	meanFEVs[6]	3.327	0.0843	4.006E-4	3.161	3.327	3.494	10000	50001
	meanFEVs[7]	3.384	0.1034	4.703E-4	3.182	3.384	3.589	10000	50001
	meanFEVs[8]	3.442	0.1278	5.646E-4	3.192	3.441	3.693	10000	50001
	meanFEVs[9]	3.499	0.1551	6.732E-4	3.197	3.499	3.804	10000	50001
	meanFEVs[10]	3.557	0.184	7.903E-4	3.2	3.556	3.919	10000	50001
	tau	2.269	0.1733	8.043E-4	1.94	2.265	2.62	10000	50001

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
FEV	697.155	692.809	4.346	701.501	
total	697.155	692.809	4.346	701.501	

RESULTS using diffuse prior

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	FEV20ns	4.869	0.6876	0.003126	3.525	4.871	6.209	10000	50001
	FEV20s	3.72	0.7094	0.002973	2.316	3.721	5.103	10000	50001
	MD	0.1897	0.188	8.951E-4	-0.1796	0.1898	0.5601	10000	50001
	RM	1.243	0.07402	3.376E-4	1.106	1.24	1.398	10000	50001
	beta[1]	0.6727	0.2519	0.001117	0.1792	0.673	1.163	10000	50001
	beta[2]	0.2099	0.02106	9.471E-5	0.1691	0.2099	0.2511	10000	50001
	beta[3]	1.718	0.5655	0.002416	0.6095	1.721	2.817	10000	50001
	beta[4]	-0.1432	0.04233	1.832E-4	-0.2256	-0.1434	-0.05966	10000	50001
	meanFEVns[1]	2.772	0.05496	2.397E-4	2.664	2.772	2.879	10000	50001
	meanFEVns[2]	2.982	0.04308	1.924E-4	2.897	2.982	3.066	10000	50001
	meanFEVns[3]	3.192	0.03974	1.858E-4	3.113	3.192	3.27	10000	50001
	meanFEVns[4]	3.402	0.0468	2.235E-4	3.31	3.402	3.494	10000	50001
	meanFEVns[5]	3.612	0.06073	2.887E-4	3.493	3.611	3.731	10000	50001
	meanFEVns[6]	3.821	0.07793	3.67E-4	3.669	3.821	3.974	10000	50001
	meanFEVns[7]	4.031	0.09668	4.516E-4	3.843	4.031	4.221	10000	50001
	meanFEVns[8]	4.241	0.1162	5.396E-4	4.015	4.241	4.469	10000	50001
	meanFEVns[9]	4.451	0.1362	6.295E-4	4.186	4.451	4.718	10000	50001
	meanFEVns[10]	4.661	0.1565	7.207E-4	4.357	4.661	4.967	10000	50001
	meanFEVs[1]	3.058	0.1563	6.772E-4	2.749	3.058	3.365	10000	50001
	meanFEVs[2]	3.124	0.1268	5.432E-4	2.874	3.125	3.372	10000	50001
	meanFEVs[3]	3.191	0.1021	4.337E-4	2.99	3.191	3.391	10000	50001
	meanFEVs[4]	3.258	0.08637	3.713E-4	3.089	3.258	3.426	10000	50001
	meanFEVs[5]	3.324	0.08486	3.798E-4	3.159	3.325	3.491	10000	50001
	meanFEVs[6]	3.391	0.09821	4.553E-4	3.199	3.391	3.583	10000	50001
	meanFEVs[7]	3.458	0.1216	5.719E-4	3.218	3.457	3.698	10000	50001
	meanFEVs[8]	3.525	0.1505	7.096E-4	3.229	3.524	3.821	10000	50001
	meanFEVs[9]	3.591	0.1822	8.582E-4	3.231	3.591	3.952	10000	50001
	meanFEVs[10]	3.658	0.2155	0.001013	3.231	3.658	4.084	10000	50001
	tau	2.268	0.1744	8.091E-4	1.937	2.264	2.621	10000	50001

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
FEV	697.348	692.330	5.018	702.366	
total	697.348	692.330	5.018	702.366	

Case 3.  Quadratic trend in age for smokers and nonsmokers

Model

model{
for(i in 1:n){
	FEV[i] ~ dnorm(mu[i],tau)
	mu[i] <- beta[1] + beta[2]*Age[i] + beta[3]*Age[i]*Age[i] + beta[4]*Smoke[i]  + beta[5]*Age[i]*Smoke[i] + beta[6]*Age[i]*Age[i]*Smoke[i] 
    }
for(j in 1:6){
   beta[j] ~ dnorm(0,0.001)  
   test[j] <- step(beta[j])
   }
tau ~ dgamma(0.001,0.001)
}
DATA

list(n=345) # Load for diffuse prior

Age[]	FEV[]	Smoke[]
19	5.102	0
19	3.519	1
19	3.345	1
18	2.906	0
18	3.082	0
18	4.22	0
18	4.086	1
18	4.404	1
18	2.853	0
17	4.429	0
17	3.082	1
17	4.724	0
17	3.406	1
17	3.5	0
17	5.633	0
17	3.96	0
17	5.638	0
16	3.688	1
16	3.387	0
16	2.903	1
16	3.674	0
16	2.608	1
16	3.645	0
16	4.07	1
16	4.299	0
16	2.981	0
16	4.504	0
16	4.872	1
16	4.27	1
16	2.795	1
15	4.284	0
15	4.506	1
15	4.279	0
15	4.5	0
15	2.635	0
15	2.679	1
15	2.198	1
15	3.004	1
15	5.793	0
15	3.985	0
15	3.731	0
15	3.122	1
15	3.33	1
15	3.799	1
15	2.887	0
15	2.264	1
15	2.278	1
15	3.727	1
15	3.211	0
14	3.381	0
14	3.741	0
14	4.842	0
14	2.236	1
14	3.436	0
14	3.428	1
14	4.683	0
14	3.957	1
14	3.074	1
14	4.273	0
14	3.395	0
14	3.169	0
14	4.309	1
14	4.111	0
14	2.538	0
14	3.68	0
14	3.585	0
14	2.891	0
14	3.78	0
14	4.271	0
14	2.934	0
14	2.276	1
14	2.997	0
14	3.806	0
14	4.763	1
13	4.336	0
13	3.816	0
13	3.549	0
13	3.147	0
13	4.305	0
13	3.906	0
13	2.819	0
13	3.994	0
13	3.208	1
13	3.193	0
13	4.789	1
13	3.745	0
13	3.887	0
13	2.646	0
13	2.976	0
13	2.599	1
13	3.06	0
13	4.756	1
13	2.704	0
13	2.434	0
13	3.086	1
13	4.225	0
13	3.089	0
13	4.877	0
13	2.216	1
13	3.078	1
13	3.297	1
13	4.448	0
13	3.984	0
13	2.677	1
13	3.255	0
13	3.297	1
13	3.331	0
13	5.083	0
13	3.056	0
13	3.785	1
13	2.449	0
13	4.232	0
13	2.531	0
13	3.577	0
13	3.152	1
13	3.141	0
13	4.045	1
12	3.058	0
12	4.55	0
12	2.841	0
12	2.889	0
12	4.393	0
12	2.384	1
12	2.971	0
12	3.835	1
12	2.751	0
12	2.556	0
12	2.868	0
12	3.255	0
12	1.916	0
12	2.241	0
12	5.224	0
12	4.08	0
12	4.411	0
12	3.791	0
12	3.343	1
12	2.913	0
12	3.279	0
12	2.347	0
12	3.751	1
12	3.079	0
12	3.403	0
12	3.501	0
12	3.186	1
12	4.073	0
12	2.752	0
12	2.304	1
12	3.692	0
12	3.082	0
12	4.72	0
12	2.417	0
12	2.759	1
12	3.231	0
12	3.529	0
12	2.866	0
12	2.605	0
12	2.569	0
12	3.924	0
12	2.752	0
12	3.53	0
12	2.332	0
12	3.035	0
12	4.831	0
12	2.714	0
12	3.519	0
12	3.341	0
12	2.822	0
12	2.935	0
12	2.499	0
12	3.001	0
12	3.222	0
12	4.203	0
12	3.082	0
12	2.579	0
11	2.884	0
11	2.17	0
11	3.47	0
11	2.524	0
11	3.654	0
11	2.665	0
11	3.222	0
11	3.49	0
11	2.957	0
11	4.007	0
11	2.386	0
11	2.762	0
11	3.011	0
11	3.583	0
11	3.236	0
11	3.058	0
11	3.108	0
11	1.694	1
11	2.346	0
11	3.515	0
11	2.754	0
11	2.463	0
11	2.633	0
11	3.111	0
11	3.977	0
11	3.411	0
11	3.171	0
11	3.587	0
11	3.845	0
11	2.417	0
11	2.735	0
11	4.065	0
11	2.318	0
11	3.596	0
11	2.542	0
11	2.354	0
11	2.491	0
11	2.993	0
11	3.774	0
11	2.988	0
11	2.498	0
11	2.887	0
11	3.515	0
11	3.425	0
11	4.593	0
11	3.223	0
11	4.073	0
11	2.606	0
11	3.169	1
11	2.465	0
11	2.827	0
11	2.578	0
11	2.081	0
11	2.974	0
11	3.102	1
11	3.023	0
11	3.681	0
11	3.258	0
11	3.247	0
11	4.324	0
11	2.362	0
11	2.563	0
11	3.206	0
11	2.953	1
11	3.078	0
11	3.369	0
11	3.022	0
11	2.866	0
11	2.501	0
11	3.32	0
11	2.123	0
11	3.847	0
11	2.928	0
11	2.689	0
11	2.894	0
11	4.637	1
11	2.812	0
11	2.387	0
11	4.13	0
11	3.28	0
11	2.659	0
11	2.822	0
11	2.14	0
11	3.12	1
11	2.562	0
11	3.339	1
11	2.458	0
11	3.104	1
11	3.069	1
11	2.785	0
10	2.328	0
10	1.811	0
10	2.642	0
10	3.166	0
10	2.561	0
10	2.481	0
10	3.203	0
10	3.111	0
10	2.52	0
10	2.292	0
10	2.246	0
10	1.937	0
10	2.646	0
10	3.251	0
10	3.007	0
10	3.489	0
10	2.864	0
10	2.25	0
10	2.352	0
10	2.592	0
10	2.72	0
10	3.048	0
10	2.094	0
10	3.183	0
10	3.354	0
10	2.387	1
10	2.504	0
10	2.891	0
10	1.823	0
10	2.175	0
10	2.673	0
10	2.608	0
10	1.458	0
10	3.795	0
10	2.545	0
10	3.305	0
10	2.855	0
10	2.287	0
10	2.365	0
10	2.696	0
10	2.813	0
10	3.413	1
10	1.858	0
10	2.975	1
10	3.35	0
10	2.901	0
10	3.2	0
10	2.358	0
10	2.581	0
10	2.691	0
10	1.873	0
10	2.758	0
10	3.05	0
10	2.201	0
10	1.858	0
10	1.665	0
10	2.25	0
10	2.862	0
10	2.356	0
10	4.591	0
10	2.216	0
10	3.498	1
10	2.364	0
10	2.341	0
10	3.127	0
10	2.132	0
10	3.456	0
10	3.073	0
10	2.688	0
10	3.329	0
10	3.11	0
10	2.435	0
10	2.838	0
10	3.086	0
10	2.77	0
10	3.09	0
10	3.038	1
10	2.568	0
10	3.132	0
10	2.391	0
10	2.1	0
END

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

RESULTS 

node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	beta[1]	-2.664	1.396	0.004109	-5.399	-2.66	0.08227	10000	110001
	beta[2]	0.7343	0.2168	6.425E-4	0.3093	0.734	1.159	10000	110001
	beta[3]	-0.01997	0.008213	2.461E-5	-0.03604	-0.01997	-0.003905	10000	110001
	beta[4]	5.428	2.973	0.008598	-0.3879	5.428	11.24	10000	110001
	beta[5]	-0.722	0.4372	0.001264	-1.578	-0.7215	0.1341	10000	110001
	beta[6]	0.0219	0.01579	4.579E-5	-0.009024	0.0219	0.05293	10000	110001
	tau	2.295	0.1764	5.6E-4	1.962	2.291	2.652	10000	110001
	test[1]	0.02857	0.1666	5.407E-4	0.0	0.0	1.0	10000	110001
	test[2]	0.9997	0.01834	5.958E-5	1.0	1.0	1.0	10000	110001
	test[3]	0.0077	0.08741	2.696E-4	0.0	0.0	0.0	10000	110001
	test[4]	0.9664	0.1802	5.515E-4	0.0	1.0	1.0	10000	110001
	test[5]	0.04955	0.217	6.984E-4	0.0	0.0	1.0	10000	110001
	test[6]	0.9176	0.275	7.979E-4	0.0	1.0	1.0	10000	110001

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
	Dbar	Dhat	pD	DIC	
FEV	693.329	686.333	6.996	700.326	
total	693.329	686.333	6.996	700.326	

Case 4.  Convergence and mixing diagnostics.

model{
for(i in 1:n){
	FEV[i] ~ dnorm(mu[i],tau)
	mu[i] <- beta[1] + beta[2]*Age[i] + beta[3]*Smoke[i]  + beta[4]*Age[i]*Smoke[i]
    }
beta[1] ~ dnorm(0,0.001)  # Diffuse prior
beta[2] ~ dnorm(0,0.001)
beta[3] ~ dnorm(0,0.001)
beta[4] ~ dnorm(0,0.001)
tau ~ dgamma(0.001,0.001)
}
DATA  

list(n=345)  # Load for diffuse prior

Age[]	FEV[]	Smoke[]
19	5.102	0
19	3.519	1
19	3.345	1
18	2.906	0
18	3.082	0
18	4.22	0
18	4.086	1
18	4.404	1
18	2.853	0
17	4.429	0
17	3.082	1
17	4.724	0
17	3.406	1
17	3.5	0
17	5.633	0
17	3.96	0
17	5.638	0
16	3.688	1
16	3.387	0
16	2.903	1
16	3.674	0
16	2.608	1
16	3.645	0
16	4.07	1
16	4.299	0
16	2.981	0
16	4.504	0
16	4.872	1
16	4.27	1
16	2.795	1
15	4.284	0
15	4.506	1
15	4.279	0
15	4.5	0
15	2.635	0
15	2.679	1
15	2.198	1
15	3.004	1
15	5.793	0
15	3.985	0
15	3.731	0
15	3.122	1
15	3.33	1
15	3.799	1
15	2.887	0
15	2.264	1
15	2.278	1
15	3.727	1
15	3.211	0
14	3.381	0
14	3.741	0
14	4.842	0
14	2.236	1
14	3.436	0
14	3.428	1
14	4.683	0
14	3.957	1
14	3.074	1
14	4.273	0
14	3.395	0
14	3.169	0
14	4.309	1
14	4.111	0
14	2.538	0
14	3.68	0
14	3.585	0
14	2.891	0
14	3.78	0
14	4.271	0
14	2.934	0
14	2.276	1
14	2.997	0
14	3.806	0
14	4.763	1
13	4.336	0
13	3.816	0
13	3.549	0
13	3.147	0
13	4.305	0
13	3.906	0
13	2.819	0
13	3.994	0
13	3.208	1
13	3.193	0
13	4.789	1
13	3.745	0
13	3.887	0
13	2.646	0
13	2.976	0
13	2.599	1
13	3.06	0
13	4.756	1
13	2.704	0
13	2.434	0
13	3.086	1
13	4.225	0
13	3.089	0
13	4.877	0
13	2.216	1
13	3.078	1
13	3.297	1
13	4.448	0
13	3.984	0
13	2.677	1
13	3.255	0
13	3.297	1
13	3.331	0
13	5.083	0
13	3.056	0
13	3.785	1
13	2.449	0
13	4.232	0
13	2.531	0
13	3.577	0
13	3.152	1
13	3.141	0
13	4.045	1
12	3.058	0
12	4.55	0
12	2.841	0
12	2.889	0
12	4.393	0
12	2.384	1
12	2.971	0
12	3.835	1
12	2.751	0
12	2.556	0
12	2.868	0
12	3.255	0
12	1.916	0
12	2.241	0
12	5.224	0
12	4.08	0
12	4.411	0
12	3.791	0
12	3.343	1
12	2.913	0
12	3.279	0
12	2.347	0
12	3.751	1
12	3.079	0
12	3.403	0
12	3.501	0
12	3.186	1
12	4.073	0
12	2.752	0
12	2.304	1
12	3.692	0
12	3.082	0
12	4.72	0
12	2.417	0
12	2.759	1
12	3.231	0
12	3.529	0
12	2.866	0
12	2.605	0
12	2.569	0
12	3.924	0
12	2.752	0
12	3.53	0
12	2.332	0
12	3.035	0
12	4.831	0
12	2.714	0
12	3.519	0
12	3.341	0
12	2.822	0
12	2.935	0
12	2.499	0
12	3.001	0
12	3.222	0
12	4.203	0
12	3.082	0
12	2.579	0
11	2.884	0
11	2.17	0
11	3.47	0
11	2.524	0
11	3.654	0
11	2.665	0
11	3.222	0
11	3.49	0
11	2.957	0
11	4.007	0
11	2.386	0
11	2.762	0
11	3.011	0
11	3.583	0
11	3.236	0
11	3.058	0
11	3.108	0
11	1.694	1
11	2.346	0
11	3.515	0
11	2.754	0
11	2.463	0
11	2.633	0
11	3.111	0
11	3.977	0
11	3.411	0
11	3.171	0
11	3.587	0
11	3.845	0
11	2.417	0
11	2.735	0
11	4.065	0
11	2.318	0
11	3.596	0
11	2.542	0
11	2.354	0
11	2.491	0
11	2.993	0
11	3.774	0
11	2.988	0
11	2.498	0
11	2.887	0
11	3.515	0
11	3.425	0
11	4.593	0
11	3.223	0
11	4.073	0
11	2.606	0
11	3.169	1
11	2.465	0
11	2.827	0
11	2.578	0
11	2.081	0
11	2.974	0
11	3.102	1
11	3.023	0
11	3.681	0
11	3.258	0
11	3.247	0
11	4.324	0
11	2.362	0
11	2.563	0
11	3.206	0
11	2.953	1
11	3.078	0
11	3.369	0
11	3.022	0
11	2.866	0
11	2.501	0
11	3.32	0
11	2.123	0
11	3.847	0
11	2.928	0
11	2.689	0
11	2.894	0
11	4.637	1
11	2.812	0
11	2.387	0
11	4.13	0
11	3.28	0
11	2.659	0
11	2.822	0
11	2.14	0
11	3.12	1
11	2.562	0
11	3.339	1
11	2.458	0
11	3.104	1
11	3.069	1
11	2.785	0
10	2.328	0
10	1.811	0
10	2.642	0
10	3.166	0
10	2.561	0
10	2.481	0
10	3.203	0
10	3.111	0
10	2.52	0
10	2.292	0
10	2.246	0
10	1.937	0
10	2.646	0
10	3.251	0
10	3.007	0
10	3.489	0
10	2.864	0
10	2.25	0
10	2.352	0
10	2.592	0
10	2.72	0
10	3.048	0
10	2.094	0
10	3.183	0
10	3.354	0
10	2.387	1
10	2.504	0
10	2.891	0
10	1.823	0
10	2.175	0
10	2.673	0
10	2.608	0
10	1.458	0
10	3.795	0
10	2.545	0
10	3.305	0
10	2.855	0
10	2.287	0
10	2.365	0
10	2.696	0
10	2.813	0
10	3.413	1
10	1.858	0
10	2.975	1
10	3.35	0
10	2.901	0
10	3.2	0
10	2.358	0
10	2.581	0
10	2.691	0
10	1.873	0
10	2.758	0
10	3.05	0
10	2.201	0
10	1.858	0
10	1.665	0
10	2.25	0
10	2.862	0
10	2.356	0
10	4.591	0
10	2.216	0
10	3.498	1
10	2.364	0
10	2.341	0
10	3.127	0
10	2.132	0
10	3.456	0
10	3.073	0
10	2.688	0
10	3.329	0
10	3.11	0
10	2.435	0
10	2.838	0
10	3.086	0
10	2.77	0
10	3.09	0
10	3.038	1
10	2.568	0
10	3.132	0
10	2.391	0
10	2.1	0
END

INITS
list(tau=1, beta=c(0,0,0,0))
list(tau=.1, beta=c(10,10,10,10))
list(tau=5, beta=c(-10,-10,-10,-10))
list(tau=10, beta=c(10,-10,10,-10))
list(tau=15, beta=c(-10,10,-10,10))
list(tau=20, beta=c(10,-10,0,0))
list(tau=.5, beta=c(0,0,5,-10))
list(tau=1,  beta=c(5,5,-10,0))

Exercise 9.19.  Remove Case 1, the person with the highest Cook's distance

Quadratic trend in age for nonsmokers; linear trend in age for smokers
model{
for(i in 1:n){
	FEV[i] ~ dnorm(mu[i],tau)
	mu[i] <- beta[1] + beta[2]*Age[i] + beta[3]*Age[i]*Age[i]*(1-Smoke[i]) + beta[4]*Smoke[i] + beta[5]*Age[i]*Smoke[i] 
}
for(j in 1:5){
   beta[j] ~ dnorm(0,0.001)  
   test[j] <- step(beta[j])
   }
tau ~ dgamma(0.001,0.001)
}

DATA

list(n=344) # Load for diffuse prior

Age[]	FEV[]	Smoke[]
#19	5.102	0
19	3.519	1
19	3.345	1
18	2.906	0
18	3.082	0
18	4.22	0
18	4.086	1
18	4.404	1
18	2.853	0
17	4.429	0
17	3.082	1
17	4.724	0
17	3.406	1
17	3.5	0
17	5.633	0
17	3.96	0
17	5.638	0
16	3.688	1
16	3.387	0
16	2.903	1
16	3.674	0
16	2.608	1
16	3.645	0
16	4.07	1
16	4.299	0
16	2.981	0
16	4.504	0
16	4.872	1
16	4.27	1
16	2.795	1
15	4.284	0
15	4.506	1
15	4.279	0
15	4.5	0
15	2.635	0
15	2.679	1
15	2.198	1
15	3.004	1
15	5.793	0
15	3.985	0
15	3.731	0
15	3.122	1
15	3.33	1
15	3.799	1
15	2.887	0
15	2.264	1
15	2.278	1
15	3.727	1
15	3.211	0
14	3.381	0
14	3.741	0
14	4.842	0
14	2.236	1
14	3.436	0
14	3.428	1
14	4.683	0
14	3.957	1
14	3.074	1
14	4.273	0
14	3.395	0
14	3.169	0
14	4.309	1
14	4.111	0
14	2.538	0
14	3.68	0
14	3.585	0
14	2.891	0
14	3.78	0
14	4.271	0
14	2.934	0
14	2.276	1
14	2.997	0
14	3.806	0
14	4.763	1
13	4.336	0
13	3.816	0
13	3.549	0
13	3.147	0
13	4.305	0
13	3.906	0
13	2.819	0
13	3.994	0
13	3.208	1
13	3.193	0
13	4.789	1
13	3.745	0
13	3.887	0
13	2.646	0
13	2.976	0
13	2.599	1
13	3.06	0
13	4.756	1
13	2.704	0
13	2.434	0
13	3.086	1
13	4.225	0
13	3.089	0
13	4.877	0
13	2.216	1
13	3.078	1
13	3.297	1
13	4.448	0
13	3.984	0
13	2.677	1
13	3.255	0
13	3.297	1
13	3.331	0
13	5.083	0
13	3.056	0
13	3.785	1
13	2.449	0
13	4.232	0
13	2.531	0
13	3.577	0
13	3.152	1
13	3.141	0
13	4.045	1
12	3.058	0
12	4.55	0
12	2.841	0
12	2.889	0
12	4.393	0
12	2.384	1
12	2.971	0
12	3.835	1
12	2.751	0
12	2.556	0
12	2.868	0
12	3.255	0
12	1.916	0
12	2.241	0
12	5.224	0
12	4.08	0
12	4.411	0
12	3.791	0
12	3.343	1
12	2.913	0
12	3.279	0
12	2.347	0
12	3.751	1
12	3.079	0
12	3.403	0
12	3.501	0
12	3.186	1
12	4.073	0
12	2.752	0
12	2.304	1
12	3.692	0
12	3.082	0
12	4.72	0
12	2.417	0
12	2.759	1
12	3.231	0
12	3.529	0
12	2.866	0
12	2.605	0
12	2.569	0
12	3.924	0
12	2.752	0
12	3.53	0
12	2.332	0
12	3.035	0
12	4.831	0
12	2.714	0
12	3.519	0
12	3.341	0
12	2.822	0
12	2.935	0
12	2.499	0
12	3.001	0
12	3.222	0
12	4.203	0
12	3.082	0
12	2.579	0
11	2.884	0
11	2.17	0
11	3.47	0
11	2.524	0
11	3.654	0
11	2.665	0
11	3.222	0
11	3.49	0
11	2.957	0
11	4.007	0
11	2.386	0
11	2.762	0
11	3.011	0
11	3.583	0
11	3.236	0
11	3.058	0
11	3.108	0
11	1.694	1
11	2.346	0
11	3.515	0
11	2.754	0
11	2.463	0
11	2.633	0
11	3.111	0
11	3.977	0
11	3.411	0
11	3.171	0
11	3.587	0
11	3.845	0
11	2.417	0
11	2.735	0
11	4.065	0
11	2.318	0
11	3.596	0
11	2.542	0
11	2.354	0
11	2.491	0
11	2.993	0
11	3.774	0
11	2.988	0
11	2.498	0
11	2.887	0
11	3.515	0
11	3.425	0
11	4.593	0
11	3.223	0
11	4.073	0
11	2.606	0
11	3.169	1
11	2.465	0
11	2.827	0
11	2.578	0
11	2.081	0
11	2.974	0
11	3.102	1
11	3.023	0
11	3.681	0
11	3.258	0
11	3.247	0
11	4.324	0
11	2.362	0
11	2.563	0
11	3.206	0
11	2.953	1
11	3.078	0
11	3.369	0
11	3.022	0
11	2.866	0
11	2.501	0
11	3.32	0
11	2.123	0
11	3.847	0
11	2.928	0
11	2.689	0
11	2.894	0
11	4.637	1
11	2.812	0
11	2.387	0
11	4.13	0
11	3.28	0
11	2.659	0
11	2.822	0
11	2.14	0
11	3.12	1
11	2.562	0
11	3.339	1
11	2.458	0
11	3.104	1
11	3.069	1
11	2.785	0
10	2.328	0
10	1.811	0
10	2.642	0
10	3.166	0
10	2.561	0
10	2.481	0
10	3.203	0
10	3.111	0
10	2.52	0
10	2.292	0
10	2.246	0
10	1.937	0
10	2.646	0
10	3.251	0
10	3.007	0
10	3.489	0
10	2.864	0
10	2.25	0
10	2.352	0
10	2.592	0
10	2.72	0
10	3.048	0
10	2.094	0
10	3.183	0
10	3.354	0
10	2.387	1
10	2.504	0
10	2.891	0
10	1.823	0
10	2.175	0
10	2.673	0
10	2.608	0
10	1.458	0
10	3.795	0
10	2.545	0
10	3.305	0
10	2.855	0
10	2.287	0
10	2.365	0
10	2.696	0
10	2.813	0
10	3.413	1
10	1.858	0
10	2.975	1
10	3.35	0
10	2.901	0
10	3.2	0
10	2.358	0
10	2.581	0
10	2.691	0
10	1.873	0
10	2.758	0
10	3.05	0
10	2.201	0
10	1.858	0
10	1.665	0
10	2.25	0
10	2.862	0
10	2.356	0
10	4.591	0
10	2.216	0
10	3.498	1
10	2.364	0
10	2.341	0
10	3.127	0
10	2.132	0
10	3.456	0
10	3.073	0
10	2.688	0
10	3.329	0
10	3.11	0
10	2.435	0
10	2.838	0
10	3.086	0
10	2.77	0
10	3.09	0
10	3.038	1
10	2.568	0
10	3.132	0
10	2.391	0
10	2.1	0
END

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

Exercise 9.19.  Remove Case 132, the person with the highest Studentized residual

Quadratic trend in age for nonsmokers; linear trend in age for smokers
model{
for(i in 1:n){
	FEV[i] ~ dnorm(mu[i],tau)
	mu[i] <- beta[1] + beta[2]*Age[i] + beta[3]*Age[i]*Age[i]*(1-Smoke[i]) + beta[4]*Smoke[i] + beta[5]*Age[i]*Smoke[i] 
}
for(j in 1:5){
   beta[j] ~ dnorm(0,0.001)  
   test[j] <- step(beta[j])
   }
tau ~ dgamma(0.001,0.001)
}

DATA

list(n=344) # Load for diffuse prior

Age[]	FEV[]	Smoke[]
19	5.102	0
19	3.519	1
19	3.345	1
18	2.906	0
18	3.082	0
18	4.22	0
18	4.086	1
18	4.404	1
18	2.853	0
17	4.429	0
17	3.082	1
17	4.724	0
17	3.406	1
17	3.5	0
17	5.633	0
17	3.96	0
17	5.638	0
16	3.688	1
16	3.387	0
16	2.903	1
16	3.674	0
16	2.608	1
16	3.645	0
16	4.07	1
16	4.299	0
16	2.981	0
16	4.504	0
16	4.872	1
16	4.27	1
16	2.795	1
15	4.284	0
15	4.506	1
15	4.279	0
15	4.5	0
15	2.635	0
15	2.679	1
15	2.198	1
15	3.004	1
15	5.793	0
15	3.985	0
15	3.731	0
15	3.122	1
15	3.33	1
15	3.799	1
15	2.887	0
15	2.264	1
15	2.278	1
15	3.727	1
15	3.211	0
14	3.381	0
14	3.741	0
14	4.842	0
14	2.236	1
14	3.436	0
14	3.428	1
14	4.683	0
14	3.957	1
14	3.074	1
14	4.273	0
14	3.395	0
14	3.169	0
14	4.309	1
14	4.111	0
14	2.538	0
14	3.68	0
14	3.585	0
14	2.891	0
14	3.78	0
14	4.271	0
14	2.934	0
14	2.276	1
14	2.997	0
14	3.806	0
14	4.763	1
13	4.336	0
13	3.816	0
13	3.549	0
13	3.147	0
13	4.305	0
13	3.906	0
13	2.819	0
13	3.994	0
13	3.208	1
13	3.193	0
13	4.789	1
13	3.745	0
13	3.887	0
13	2.646	0
13	2.976	0
13	2.599	1
13	3.06	0
13	4.756	1
13	2.704	0
13	2.434	0
13	3.086	1
13	4.225	0
13	3.089	0
13	4.877	0
13	2.216	1
13	3.078	1
13	3.297	1
13	4.448	0
13	3.984	0
13	2.677	1
13	3.255	0
13	3.297	1
13	3.331	0
13	5.083	0
13	3.056	0
13	3.785	1
13	2.449	0
13	4.232	0
13	2.531	0
13	3.577	0
13	3.152	1
13	3.141	0
13	4.045	1
12	3.058	0
12	4.55	0
12	2.841	0
12	2.889	0
12	4.393	0
12	2.384	1
12	2.971	0
12	3.835	1
12	2.751	0
12	2.556	0
12	2.868	0
12	3.255	0
12	1.916	0
12	2.241	0
#12	5.224	0
12	4.08	0
12	4.411	0
12	3.791	0
12	3.343	1
12	2.913	0
12	3.279	0
12	2.347	0
12	3.751	1
12	3.079	0
12	3.403	0
12	3.501	0
12	3.186	1
12	4.073	0
12	2.752	0
12	2.304	1
12	3.692	0
12	3.082	0
12	4.72	0
12	2.417	0
12	2.759	1
12	3.231	0
12	3.529	0
12	2.866	0
12	2.605	0
12	2.569	0
12	3.924	0
12	2.752	0
12	3.53	0
12	2.332	0
12	3.035	0
12	4.831	0
12	2.714	0
12	3.519	0
12	3.341	0
12	2.822	0
12	2.935	0
12	2.499	0
12	3.001	0
12	3.222	0
12	4.203	0
12	3.082	0
12	2.579	0
11	2.884	0
11	2.17	0
11	3.47	0
11	2.524	0
11	3.654	0
11	2.665	0
11	3.222	0
11	3.49	0
11	2.957	0
11	4.007	0
11	2.386	0
11	2.762	0
11	3.011	0
11	3.583	0
11	3.236	0
11	3.058	0
11	3.108	0
11	1.694	1
11	2.346	0
11	3.515	0
11	2.754	0
11	2.463	0
11	2.633	0
11	3.111	0
11	3.977	0
11	3.411	0
11	3.171	0
11	3.587	0
11	3.845	0
11	2.417	0
11	2.735	0
11	4.065	0
11	2.318	0
11	3.596	0
11	2.542	0
11	2.354	0
11	2.491	0
11	2.993	0
11	3.774	0
11	2.988	0
11	2.498	0
11	2.887	0
11	3.515	0
11	3.425	0
11	4.593	0
11	3.223	0
11	4.073	0
11	2.606	0
11	3.169	1
11	2.465	0
11	2.827	0
11	2.578	0
11	2.081	0
11	2.974	0
11	3.102	1
11	3.023	0
11	3.681	0
11	3.258	0
11	3.247	0
11	4.324	0
11	2.362	0
11	2.563	0
11	3.206	0
11	2.953	1
11	3.078	0
11	3.369	0
11	3.022	0
11	2.866	0
11	2.501	0
11	3.32	0
11	2.123	0
11	3.847	0
11	2.928	0
11	2.689	0
11	2.894	0
11	4.637	1
11	2.812	0
11	2.387	0
11	4.13	0
11	3.28	0
11	2.659	0
11	2.822	0
11	2.14	0
11	3.12	1
11	2.562	0
11	3.339	1
11	2.458	0
11	3.104	1
11	3.069	1
11	2.785	0
10	2.328	0
10	1.811	0
10	2.642	0
10	3.166	0
10	2.561	0
10	2.481	0
10	3.203	0
10	3.111	0
10	2.52	0
10	2.292	0
10	2.246	0
10	1.937	0
10	2.646	0
10	3.251	0
10	3.007	0
10	3.489	0
10	2.864	0
10	2.25	0
10	2.352	0
10	2.592	0
10	2.72	0
10	3.048	0
10	2.094	0
10	3.183	0
10	3.354	0
10	2.387	1
10	2.504	0
10	2.891	0
10	1.823	0
10	2.175	0
10	2.673	0
10	2.608	0
10	1.458	0
10	3.795	0
10	2.545	0
10	3.305	0
10	2.855	0
10	2.287	0
10	2.365	0
10	2.696	0
10	2.813	0
10	3.413	1
10	1.858	0
10	2.975	1
10	3.35	0
10	2.901	0
10	3.2	0
10	2.358	0
10	2.581	0
10	2.691	0
10	1.873	0
10	2.758	0
10	3.05	0
10	2.201	0
10	1.858	0
10	1.665	0
10	2.25	0
10	2.862	0
10	2.356	0
10	4.591	0
10	2.216	0
10	3.498	1
10	2.364	0
10	2.341	0
10	3.127	0
10	2.132	0
10	3.456	0
10	3.073	0
10	2.688	0
10	3.329	0
10	3.11	0
10	2.435	0
10	2.838	0
10	3.086	0
10	2.77	0
10	3.09	0
10	3.038	1
10	2.568	0
10	3.132	0
10	2.391	0
10	2.1	0
END

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