Chapter 7:  Basic concepts of regression
FEV data from Rosner; Fundamentals of Biostatistics; 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 in years range from 10 to 19 

Exercise 7.11: Bivariate association

model
{
for(i in 1:n){
	FEV[i] ~ dnorm(mu[Smoke[i]+1], tau[Smoke[i]+1])
}
mu[1] ~ dnorm(0,0.001)
mu[2] ~ dnorm(0,0.001)
tau[1] ~ dgamma(0.001,0.001)
tau[2] ~ dgamma(0.001,0.001)
theta <- mu[2]-mu[1]   # smokers - nonsmokers
thetadiff <- step(theta)
junk1 <- Age[1]
}
DATA

list(n=345)
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=c(1,1), mu=c(0,0))

Main effects ACOVA model.  Diffuse priors for beta and tau.

model
{
for(i in 1:n){
	FEV[i] ~ dnorm(mu[i],tau)
	mu[i] <- beta[1] + beta[2]*Smoke[i] + beta[3]*Age[i]
}
beta[1] ~ dnorm(0,0.001)
beta[2] ~ dnorm(0,0.001)
beta[3] ~ dnorm(0,0.001)
tau ~ dgamma(0.001,0.001)
thetadiff <- 1-step(beta[2])
}
DATA
Use same data as above
INITS
list(tau=1, beta=c(0,0,0))
RESULTS

	 node	 mean	 sd	     MC error	2.5%	median	97.5%	   start	sample
	beta[1]	1.095	  0.2217	0.001093	0.6666    1.094	   1.53	      10000	50001
	beta[2]	-0.1648    0.09948   4.415E-4	-0.3601   -0.1651	0.03024	10000	50001
	beta[3]	0.1741	0.01845   8.984E-5	0.1378    0.1741	  0.21	      10000	50001
	tau	      2.201	  0.1681	7.812E-4	1.884	 2.196	   2.545	     10000	50001
	thetadiff   0.9511	0.2157	9.307E-4	0.0	    1.0	        1.0	         10000	50001
	
	Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
					Dbar		Dhat		pD		DIC	
		FEV	707.810	703.797	4.013	711.823	
		total	707.810	703.797	4.013	711.823	

Main effects and interaction.  Diffuse priors for beta and tau. 

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)
beta[2] ~ dnorm(0,0.001)
beta[3] ~ dnorm(0,0.001)
beta[4] ~ dnorm(0,0.001)
tau ~ dgamma(0.001,0.001)
int <- step(beta[4])

## 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 compute relative means and mean difference as well
RM <- meanFEVns[9]  / meanFEVs[9]    ## RM comparing 18 year old smoker to 18 year old non-smoker
MD1 <- meanFEVs[9]  - meanFEVns[4]   ## MD comparing 18 year old smoker to 13 year old nonsmoker
MD2 <- meanFEVs[2]  - meanFEVns[2]   ## MD comparing 11 year old smoker to 11 year old nonsmoker
MD3 <- meanFEVs[8]  - meanFEVns[8]   ## MD comparing 17 year old smoker to 17 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
Use same data as above

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

RESULTS

     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
	MD1	0.1897	0.188	8.951E-4	-0.1796	0.1898	0.5601	10000	50001
	MD2	0.1426	0.1341	5.715E-4	-0.1208	0.1434	0.4054	10000	50001
	MD3	-0.7167	0.1901	8.628E-4	-1.09	-0.7169	-0.3408	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
	int	4.2E-4	0.02049	9.609E-5	0.0	0.0	0.0	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