Chapter 10: Correlated Data

############################
## IL-1beta longitudinal data    ##
############################

1a.  MVN specification with independence covariance matrix.
        Biomarker (y=IL-1beta) measured from saliva at 6 time points for 30 healthy subjects.  

model{ 
for( j in 1:30){
     for(i in 1:6){
           logy[ j, i ] <-  log(y[ j, i ])
     }
  }
  for( j in 1:30){
    logy[ j, 1:6 ] ~ dmnorm(m[1:6], precision[1:6,1:6])
  }
  for(j in 1:6){
	for(k in 1:6){
	      covariance[j,k] <-  pow(sigma,2)*equals(j , k) 
	}
  }
 for(i in 1:6){   m[ i ] <- mu  }
 precision[1:6,1:6] <- inverse(covariance[1:6,1:6])
 sigma ~dunif(0,100)
 tau <- 1/(sigma*sigma)
 mu ~ dnorm(0,0.001)
}

Inits
list(mu =0,       sigma = 1)
list(mu = -50,  sigma = 0.1)
list(mu = -30,  sigma = 5)
list(mu = -10,  sigma = 10)
list(mu = 100, sigma = 15)

IL-1beta Data

y[,1]	y[,2]    y[,3]     y[,4]   y[,5]   y[,6]  
22.8	17.6	11.6	7	8.78	6.6
3.14	14	7.13	3.66	1.992	7.46
12.7	5.63	6.4	4.14	5.3	0.784
48.2	20.1	6.21	6.26	8.12	4.06
3.09	4.21	3.14	3.04	0.656	3.04
26.7	14.8	11.8	23.8	30.6	12.14
1.31	1.31	1.8	10.86	3.38	5.44
30.2	29.3	12.5	17.36	9.06	11.88
1.87	37.5	3.86	7.94	0.402	4.48
2.72	1.87	1.87	3.66	4.06	19.48
36.9	81.3	25.1	38.8	59.4	16.8
24	13.5	24.2	20.6	14.86	6.66
3.14	1.87	0.772	2.67	1.87	0.289
11.78	8.32	17.68	4.47	4.47	5.88
2.62	1.834	0.1644	0.489	1.87	4.88
1.092	1.91	1.782	1.71	4.01	4.31
37.2	19.54	11.26	8.05	7.18	29
0.214	7.58	10.5	5.38	25.9	26.9
1.87	2.22	0.64	0.887	1.87	1.87
2.88	3.7	0.64	1.15	1.96	9.43
9.06	17.8	1.91	5.77	13.3	7.85
29.6	10.54	12.6	9.77	0.108	6.86
10.1	3.42	14.58	3.18	1.83	3.3
6.5	1.962	5.42	3.37	5.04	2.21
2.56	3.24	1.284	1.87	1.87	6.69
0.581	2.14	3.77	7.3	4.82	1.28
1.51	1.51	0.799	1.83	0.48	7.3
11.2	3.56	29	3.42	15.9	4.99
8.86	3.09	26.1	1.67	3.83	2.67
1.51	2.29	1.45	1.67	3.96	1.67
END

Results

 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	mu	1.564	0.09029	3.765E-4	1.387	1.565	1.741	4001	50000
	sigma	1.211	0.06423	2.95E-4	1.092	1.208	1.344	4001	50000
	tau	0.688	0.07255	3.265E-4	0.5533	0.6854	0.8386	4001	50000

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
		Dbar			Dhat		pD		DIC	
logy	578.199	576.220	1.978	580.177	
total	578.199	576.220	1.978	580.177	

1b.  One-way random effects model.  Biomarker (y=IL-1beta) measured from saliva at 6 time points for 30 healthy subjects. 

model{
  for(i in 1:180){
     logy[i] <- log(y[i])
     logy[i] ~ dnorm(m[i],tau)
     m[i] <- mu + b[ID[i]]
   }
  for(i in 1:30){
     b[i] ~ dnorm(0,taub)
   }
 sigma ~dunif(0,100)
 tau <- 1/(sigma*sigma)
 sigmab ~ dunif(0,100) 
 taub <- 1/(sigmab*sigmab)
 mu ~ dnorm(0,0.001)
 rho <- tau/(taub+tau)
}

Inits
list(mu=0,rho=0.5,sigmab=1)
list(mu=-50,sigma=.1,sigmab=25)
list(mu=-30,sigma=5,sigmab=20)
list(mu=-10,sigma=10,sigmab=15)
list(mu=100,sigma=15,sigmab=10)
list(mu=30,sigma=20,sigmab=5)
list(mu=50,sigma=25,sigmab=.1)

IL-1beta Data

ID[]	y[]
1	22.8
1	17.6
1	11.6
1	7
1	8.78
1	6.6
2	3.14
2	14
2	7.13
2	3.66
2	1.992
2	7.46
3	12.7
3	5.63
3	6.4
3	4.14
3	5.3
3	0.784
4	48.2
4	20.1
4	6.21
4	6.26
4	8.12
4	4.06
5	3.09
5	4.21
5	3.14
5	3.04
5	0.656
5	3.04
6	26.7
6	14.8
6	11.8
6	23.8
6	30.6
6	12.14
7	1.31
7	1.31
7	1.8
7	10.86
7	3.38
7	5.44
8	30.2
8	29.3
8	12.5
8	17.36
8	9.06
8	11.88
9	1.87
9	37.5
9	3.86
9	7.94
9	0.402
9	4.48
10	2.72
10	1.87
10	1.87
10	3.66
10	4.06
10	19.48
11	36.9
11	81.3
11	25.1
11	38.8
11	59.4
11	16.8
12	24
12	13.5
12	24.2
12	20.6
12	14.86
12	6.66
13	3.14
13	1.87
13	0.772
13	2.67
13	1.87
13	0.289
14	11.78
14	8.32
14	17.68
14	4.47
14	4.47
14	5.88
15	2.62
15	1.834
15	0.1644
15	0.489
15	1.87
15	4.88
16	1.092
16	1.91
16	1.782
16	1.71
16	4.01
16	4.31
17	37.2
17	19.54
17	11.26
17	8.05
17	7.18
17	29
18	0.214
18	7.58
18	10.5
18	5.38
18	25.9
18	26.9
19	1.87
19	2.22
19	0.64
19	0.887
19	1.87
19	1.87
20	2.88
20	3.7
20	0.64
20	1.15
20	1.96
20	9.43
21	9.06
21	17.8
21	1.91
21	5.77
21	13.3
21	7.85
22	29.6
22	10.54
22	12.6
22	9.77
22	0.108
22	6.86
23	10.1
23	3.42
23	14.58
23	3.18
23	1.83
23	3.3
24	6.5
24	1.962
24	5.42
24	3.37
24	5.04
24	2.21
25	2.56
25	3.24
25	1.284
25	1.87
25	1.87
25	6.69
26	0.581
26	2.14
26	3.77
26	7.3
26	4.82
26	1.28
27	1.51
27	1.51
27	0.799
27	1.83
27	0.48
27	7.3
28	11.2
28	3.56
28	29
28	3.42
28	15.9
28	4.99
29	8.86
29	3.09
29	26.1
29	1.67
29	3.83
29	2.67
30	1.51
30	2.29
30	1.45
30	1.67
30	3.96
30	1.67
END

Results

 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	mu	1.567	0.1677	0.001639	1.236	1.566	1.898	10001	100000
	rho	0.4486	0.08651	4.217E-4	0.2846	0.4469	0.6217	10001	100000
	sigma	0.9174	0.05396	2.048E-4	0.8195	0.9145	1.031	10001	100000
	sigmab	0.8348	0.1406	6.963E-4	0.597	0.8216	1.148	10001	100000
	tau	1.201	0.14	5.23E-4	0.9409	1.196	1.489	10001	100000
	taub	1.558	0.5271	0.002501	0.7591	1.481	2.806	10001	100000

 

1c.  Same model as 1b, but with a MVN specification with compound symmetric covariance.
Biomarker (y=IL-1beta) measured from saliva at 6 time points for 30 healthy subjects.  

model{ 
for( j in 1:30){
     for(i in 1:6){
           logy[ j, i ] <-  log(y[ j, i ])
     }
  }
  for( j in 1:30){
    logy[ j, 1:6 ] ~ dmnorm(m[1:6], precision[1:6,1:6])
  }
  for(j in 1:6){
	for(k in 1:6){
	      covariance[j,k] <-  pow(sigmab,2) + pow(sigma,2)*equals(j , k) 
	}
  }
 for(i in 1:6){   m[ i ] <- mu  }
 precision[1:6,1:6] <- inverse(covariance[1:6,1:6])
 sigma ~dunif(0,100)
 tau <- 1/(sigma*sigma)
 sigmab ~ dunif(0,100) 
 taub <- 1/(sigmab*sigmab)
 mu ~ dnorm(0,0.001)
 rho <- tau/(taub+tau) 
 ICC <- taub/(taub+tau)
}

Inits
list(mu=0,sigma=1,sigmab=1)
list(mu=-50,sigma=.1,sigmab=25)
list(mu=-30,sigma=5,sigmab=20)
list(mu=-10,sigma=10,sigmab=15)
list(mu=100,sigma=15,sigmab=10)

IL-1beta Data

y[,1]	y[,2]    y[,3]     y[,4]   y[,5]   y[,6]  
22.8	17.6	11.6	7	8.78	6.6
3.14	14	7.13	3.66	1.992	7.46
12.7	5.63	6.4	4.14	5.3	0.784
48.2	20.1	6.21	6.26	8.12	4.06
3.09	4.21	3.14	3.04	0.656	3.04
26.7	14.8	11.8	23.8	30.6	12.14
1.31	1.31	1.8	10.86	3.38	5.44
30.2	29.3	12.5	17.36	9.06	11.88
1.87	37.5	3.86	7.94	0.402	4.48
2.72	1.87	1.87	3.66	4.06	19.48
36.9	81.3	25.1	38.8	59.4	16.8
24	13.5	24.2	20.6	14.86	6.66
3.14	1.87	0.772	2.67	1.87	0.289
11.78	8.32	17.68	4.47	4.47	5.88
2.62	1.834	0.1644	0.489	1.87	4.88
1.092	1.91	1.782	1.71	4.01	4.31
37.2	19.54	11.26	8.05	7.18	29
0.214	7.58	10.5	5.38	25.9	26.9
1.87	2.22	0.64	0.887	1.87	1.87
2.88	3.7	0.64	1.15	1.96	9.43
9.06	17.8	1.91	5.77	13.3	7.85
29.6	10.54	12.6	9.77	0.108	6.86
10.1	3.42	14.58	3.18	1.83	3.3
6.5	1.962	5.42	3.37	5.04	2.21
2.56	3.24	1.284	1.87	1.87	6.69
0.581	2.14	3.77	7.3	4.82	1.28
1.51	1.51	0.799	1.83	0.48	7.3
11.2	3.56	29	3.42	15.9	4.99
8.86	3.09	26.1	1.67	3.83	2.67
1.51	2.29	1.45	1.67	3.96	1.67
END

Results

		 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	1/ICC	0.4487	0.08661	4.373E-4	0.2846	0.4468	0.6226	4001	50000
	mu	1.564	0.1692	7.131E-4	1.231	1.564	1.898	4001	50000
	sigma	0.9172	0.05404	2.544E-4	0.8194	0.9146	1.03	4001	50000
	sigmab	0.8349	0.141	7.109E-4	0.5991	0.8213	1.15	4001	50000
	tau	1.201	0.1405	6.381E-4	0.9424	1.195	1.489	4001	50000
	taub	1.558	0.5262	0.002561	0.7563	1.483	2.786	4001	50000
	
	Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
		Dbar		Dhat			pD		DIC	
logy	530.462	527.482	2.980	533.441	
total	530.462	527.482	2.980	533.441

1d.  Compound symmetric covariance with a uniform(-1/(k-1), 1) prior on rho
     Biomarker (y=IL-1beta) measured from saliva at 6 time points for 30 healthy 
     subjects.  
        This is the model presented in the book in Ch 10.

model{ 
for(i in 1:30){
  for(j in 1:6){
    logy[i,j] <- log(y[i,j]) 
  }
}
for( i in 1:30){  logy[i,1:6] ~ dmnorm(m[1:6],precision[1:6,1:6])  }
for(j in 1:6){ 
  for(k in 1:6){
    covariance[j,k] <- sigma2*pow(rho, step(abs(j-k)-0.5))
  }
}
for(i in 1:6){ m[i] <- mu }
precision[1:6,1:6] <- inverse(covariance[1:6,1:6])
sigma ~ dunif(0,100) 
mu ~ dnorm(0,0.001) 
L <- -1/(6-1)
rho ~ dunif(L, 1)
sigma2 <- sigma*sigma
tau <- 1/sigma2
}

Inits
list(mu=0,rho=0.5,sigma=1)

IL-1beta Data

y[,1]	y[,2]    y[,3]     y[,4]   y[,5]   y[,6]  
22.8	17.6	11.6	7	8.78	6.6
3.14	14	7.13	3.66	1.992	7.46
12.7	5.63	6.4	4.14	5.3	0.784
48.2	20.1	6.21	6.26	8.12	4.06
3.09	4.21	3.14	3.04	0.656	3.04
26.7	14.8	11.8	23.8	30.6	12.14
1.31	1.31	1.8	10.86	3.38	5.44
30.2	29.3	12.5	17.36	9.06	11.88
1.87	37.5	3.86	7.94	0.402	4.48
2.72	1.87	1.87	3.66	4.06	19.48
36.9	81.3	25.1	38.8	59.4	16.8
24	13.5	24.2	20.6	14.86	6.66
3.14	1.87	0.772	2.67	1.87	0.289
11.78	8.32	17.68	4.47	4.47	5.88
2.62	1.834	0.1644	0.489	1.87	4.88
1.092	1.91	1.782	1.71	4.01	4.31
37.2	19.54	11.26	8.05	7.18	29
0.214	7.58	10.5	5.38	25.9	26.9
1.87	2.22	0.64	0.887	1.87	1.87
2.88	3.7	0.64	1.15	1.96	9.43
9.06	17.8	1.91	5.77	13.3	7.85
29.6	10.54	12.6	9.77	0.108	6.86
10.1	3.42	14.58	3.18	1.83	3.3
6.5	1.962	5.42	3.37	5.04	2.21
2.56	3.24	1.284	1.87	1.87	6.69
0.581	2.14	3.77	7.3	4.82	1.28
1.51	1.51	0.799	1.83	0.48	7.3
11.2	3.56	29	3.42	15.9	4.99
8.86	3.09	26.1	1.67	3.83	2.67
1.51	2.29	1.45	1.67	3.96	1.67
END

Results

 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	mu	1.564	0.1685	0.001305	1.234	1.563	1.902	1001	20000
	rho	0.4473	0.08428	0.001032	0.2876	0.4458	0.6141	1001	20000
	sigma	1.24	0.09798	0.00124	1.077	1.23	1.458	1001	20000
	tau	0.6617	0.1003	0.001231	0.4701	0.6605	0.8615	1001	20000

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
	Dbar	Dhat	pD	DIC	
logy	530.372	527.487	2.885	533.257	
total	530.372	527.487	2.885	533.257	

1e.  Compound symmetric precision,  with a uniform(-1/(k-1), 1) prior on rho
     Biomarker (y=IL-1beta) measured from saliva at 6 time points for 30 healthy   
     subjects.

model{ 
for(i in 1:30){
  for(j in 1:6){
    logy[i,j] <- log(y[i,j]) 
  }
}
for( i in 1:30){
  logy[i,1:6] ~ dmnorm(m[1:6],precision[1:6,1:6]) 
}
for(s in 1:6){ 
  for(t in 1:6){
    precision[s, t] <- 1/(sigma2*(1-rho))*equals(s,t) - rho/(sigma2*(1-rho)*(1+(k-1)*rho)) 
  }
}
for(i in 1:6){ m[i] <- mu }
sigma ~ dunif(0,100) 
mu ~ dnorm(0,0.001) 
L <- -1/(6-1)
rho ~ dunif(L, 1)
sigma2 <- sigma*sigma
tau <- 1/sigma2
k <- 6
}

Inits
list(mu=0,rho=0.5,sigma=1)

IL-1beta Data

y[,1]	y[,2]    y[,3]     y[,4]   y[,5]   y[,6]  
22.8	17.6	11.6	7	8.78	6.6
3.14	14	7.13	3.66	1.992	7.46
12.7	5.63	6.4	4.14	5.3	0.784
48.2	20.1	6.21	6.26	8.12	4.06
3.09	4.21	3.14	3.04	0.656	3.04
26.7	14.8	11.8	23.8	30.6	12.14
1.31	1.31	1.8	10.86	3.38	5.44
30.2	29.3	12.5	17.36	9.06	11.88
1.87	37.5	3.86	7.94	0.402	4.48
2.72	1.87	1.87	3.66	4.06	19.48
36.9	81.3	25.1	38.8	59.4	16.8
24	13.5	24.2	20.6	14.86	6.66
3.14	1.87	0.772	2.67	1.87	0.289
11.78	8.32	17.68	4.47	4.47	5.88
2.62	1.834	0.1644	0.489	1.87	4.88
1.092	1.91	1.782	1.71	4.01	4.31
37.2	19.54	11.26	8.05	7.18	29
0.214	7.58	10.5	5.38	25.9	26.9
1.87	2.22	0.64	0.887	1.87	1.87
2.88	3.7	0.64	1.15	1.96	9.43
9.06	17.8	1.91	5.77	13.3	7.85
29.6	10.54	12.6	9.77	0.108	6.86
10.1	3.42	14.58	3.18	1.83	3.3
6.5	1.962	5.42	3.37	5.04	2.21
2.56	3.24	1.284	1.87	1.87	6.69
0.581	2.14	3.77	7.3	4.82	1.28
1.51	1.51	0.799	1.83	0.48	7.3
11.2	3.56	29	3.42	15.9	4.99
8.86	3.09	26.1	1.67	3.83	2.67
1.51	2.29	1.45	1.67	3.96	1.67
END

Results

node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	mu	1.564	0.1685	0.001305	1.234	1.563	1.902	1001	20000
	rho	0.4473	0.08428	0.001032	0.2876	0.4458	0.6141	1001	20000
	sigma	1.24	0.09798	0.00124	1.077	1.23	1.458	1001	20000
	tau	0.6617	0.1003	0.001231	0.4701	0.6605	0.8615	1001	20000

  

1f.  MVN analysis with AR(1) correlation.
    Biomarker (y=IL-1beta) measured from saliva at 6 time points for 30 healthy 
    subjects.  

model{ 
for( j in 1:30){
     for(i in 1:6){
           logy[ j, i ] <-  log(y[ j, i ])
     }
  }
  for( j in 1:30){
    logy[ j, 1:6 ] ~ dmnorm(m[1:6], precision[1:6,1:6])
  }
  for(j in 1:6){
	for(k in 1:6){
	      covariance[j,k] <-  pow(sigma,2)*pow(rho, abs(j-k)) 
	}
  }
 for(i in 1:6){   m[ i ] <- mu  }
 precision[1:6,1:6] <- inverse(covariance[1:6,1:6])
 sigma ~dunif(0,100)
 tau <- 1/(sigma*sigma)
 rho ~ dunif(-1,1) 
 mu ~ dnorm(0,0.001)
}

Inits
list(mu=0,sigma=1,rho=0)
list(mu=-50,sigma=.1,rho=.25)
list(mu=-30,sigma=5,rho=.50)
list(mu=-10,sigma=10,rho=.75)
list(mu=100,sigma=15,rho=-.25)

IL-1beta Data

y[,1]	y[,2]    y[,3]     y[,4]   y[,5]   y[,6]  
22.8	17.6	11.6	7	8.78	6.6
3.14	14	7.13	3.66	1.992	7.46
12.7	5.63	6.4	4.14	5.3	0.784
48.2	20.1	6.21	6.26	8.12	4.06
3.09	4.21	3.14	3.04	0.656	3.04
26.7	14.8	11.8	23.8	30.6	12.14
1.31	1.31	1.8	10.86	3.38	5.44
30.2	29.3	12.5	17.36	9.06	11.88
1.87	37.5	3.86	7.94	0.402	4.48
2.72	1.87	1.87	3.66	4.06	19.48
36.9	81.3	25.1	38.8	59.4	16.8
24	13.5	24.2	20.6	14.86	6.66
3.14	1.87	0.772	2.67	1.87	0.289
11.78	8.32	17.68	4.47	4.47	5.88
2.62	1.834	0.1644	0.489	1.87	4.88
1.092	1.91	1.782	1.71	4.01	4.31
37.2	19.54	11.26	8.05	7.18	29
0.214	7.58	10.5	5.38	25.9	26.9
1.87	2.22	0.64	0.887	1.87	1.87
2.88	3.7	0.64	1.15	1.96	9.43
9.06	17.8	1.91	5.77	13.3	7.85
29.6	10.54	12.6	9.77	0.108	6.86
10.1	3.42	14.58	3.18	1.83	3.3
6.5	1.962	5.42	3.37	5.04	2.21
2.56	3.24	1.284	1.87	1.87	6.69
0.581	2.14	3.77	7.3	4.82	1.28
1.51	1.51	0.799	1.83	0.48	7.3
11.2	3.56	29	3.42	15.9	4.99
8.86	3.09	26.1	1.67	3.83	2.67
1.51	2.29	1.45	1.67	3.96	1.67
END

Results

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	mu	1.589	0.1399	6.58E-4	1.312	1.589	1.865	4001	50000
	rho	0.5088	0.06671	4.502E-4	0.3748	0.5101	0.6365	4001	50000
	sigma	1.223	0.08179	5.701E-4	1.08	1.217	1.4	4001	50000
	tau	0.6774	0.08858	6.003E-4	0.5099	0.6757	0.8574	4001	50000
	
	Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
logy	535.216	532.305	2.911	538.127	
total	535.216	532.305	2.911	538.127	

1g.  MVN analysis with banded correlation.
     Biomarker (y=IL-1beta) measured from saliva at 6 time points for 30 healthy 
     subjects.  

model{ 
for( j in 1:30){
     for(i in 1:6){
           logy[ j, i ] <-  log(y[ j, i ])
     }
  }
  for( j in 1:30){
    logy[ j, 1:6 ] ~ dmnorm(m[1:6], precision[1:6,1:6])
  }
 for(j in 1:6){covariance[ j , j ] <-  pow(sigma,2) }
 for(k in 2:6){k1[k] <- k-1; covariance[1,k] <-  pow(sigma,2) * rho[k1[k]]  ; covariance[k,1] <- covariance[1,k]}
 for(k in 3:6){k2[k] <- k-2; covariance[2,k] <-  pow(sigma,2) * rho[k2[k]]  ; covariance[k,2] <- covariance[2,k]}
 for(k in 4:6){k3[k] <- k-3; covariance[3,k] <-  pow(sigma,2) * rho[k3[k]]  ; covariance[k,3] <- covariance[3,k]}
 for(k in 5:6){k4[k] <- k-4; covariance[4,k] <-  pow(sigma,2) * rho[k4[k]]  ; covariance[k,4] <- covariance[3,4]}
 covariance[5,6] <-  pow(sigma,2) * rho[1] ; covariance[6,5] <-  pow(sigma,2) * rho[1]
 for(i in 1:6){   m[ i ] <- mu  }
 precision[1:6,1:6] <- inverse(covariance[1:6,1:6])
 sigma ~dunif(0,100)
 tau <- 1/(sigma*sigma)
 mu ~ dnorm(0,0.001)
  for(i in 1:5){   rho[ i ] ~ dunif(-1,1)  }
}

Inits
list(mu=0,sigma=1,rho=c(0,0,0,0,0))

IL-1beta Data

y[,1]	y[,2]    y[,3]     y[,4]   y[,5]   y[,6]  
22.8	17.6	11.6	7	8.78	6.6
3.14	14	7.13	3.66	1.992	7.46
12.7	5.63	6.4	4.14	5.3	0.784
48.2	20.1	6.21	6.26	8.12	4.06
3.09	4.21	3.14	3.04	0.656	3.04
26.7	14.8	11.8	23.8	30.6	12.14
1.31	1.31	1.8	10.86	3.38	5.44
30.2	29.3	12.5	17.36	9.06	11.88
1.87	37.5	3.86	7.94	0.402	4.48
2.72	1.87	1.87	3.66	4.06	19.48
36.9	81.3	25.1	38.8	59.4	16.8
24	13.5	24.2	20.6	14.86	6.66
3.14	1.87	0.772	2.67	1.87	0.289
11.78	8.32	17.68	4.47	4.47	5.88
2.62	1.834	0.1644	0.489	1.87	4.88
1.092	1.91	1.782	1.71	4.01	4.31
37.2	19.54	11.26	8.05	7.18	29
0.214	7.58	10.5	5.38	25.9	26.9
1.87	2.22	0.64	0.887	1.87	1.87
2.88	3.7	0.64	1.15	1.96	9.43
9.06	17.8	1.91	5.77	13.3	7.85
29.6	10.54	12.6	9.77	0.108	6.86
10.1	3.42	14.58	3.18	1.83	3.3
6.5	1.962	5.42	3.37	5.04	2.21
2.56	3.24	1.284	1.87	1.87	6.69
0.581	2.14	3.77	7.3	4.82	1.28
1.51	1.51	0.799	1.83	0.48	7.3
11.2	3.56	29	3.42	15.9	4.99
8.86	3.09	26.1	1.67	3.83	2.67
1.51	2.29	1.45	1.67	3.96	1.67
END

Results

	node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	mu	1.591	0.1556	0.00127	1.284	1.591	1.899	4000	19962
	rho[1]	0.4621	0.08507	0.001501	0.2861	0.4654	0.619	4000	19962
	rho[2]	0.4729	0.08145	0.001476	0.3011	0.4761	0.6232	4000	19962
	rho[3]	0.2961	0.1038	0.002047	0.08718	0.2996	0.4936	4000	19962
	rho[4]	0.3035	0.1146	0.00211	0.06713	0.3065	0.5129	4000	19962
	rho[5]	0.1609	0.1538	0.002566	-0.147	0.1644	0.448	4000	19962
	sigma	1.221	0.08675	0.00137	1.069	1.215	1.411	4000	19962
	tau	0.6808	0.09447	0.001437	0.502	0.6778	0.8743	4000	19962

	
Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
		Dbar			Dhat		pD		DIC	
logy	525.183	518.996	6.187	531.370	
total	525.183	518.996	6.187	531.370	

1h.    MVN specification with MA(1) covariance matrix.
       Biomarker (y=IL-1beta) measured from saliva at 6 time points for 30 healthy 
       subjects. 

model{ 
for( j in 1:30){
     for(i in 1:6){
           logy[ j, i ] <-  log(y[ j, i ])
     }
  }
  for( j in 1:30){
    logy[ j, 1:6 ] ~ dmnorm(m[1:6], precision[1:6,1:6])
  }
  for(j in 1:6){
	for(k in 1:6){
	      covariance[j,k] <-  pow(sigma,2)*equals(j , k)  + pow(sigma,2)*rho*equals(abs(j-k) , 1)
	}
  }
 for(i in 1:6){   m[ i ] <- mu  }
 precision[1:6,1:6] <- inverse(covariance[1:6,1:6])
 sigma ~dunif(0,100)
 tau <- 1/(sigma*sigma)
 mu ~ dnorm(0,0.001)
 rho ~ dunif(-0.5, 0.5)
}
Inits   
list(mu =0,       sigma = 1   ,     rho=0)

IL-1beta Data

y[,1]	y[,2]    y[,3]     y[,4]   y[,5]   y[,6]  
22.8	17.6	11.6	7	8.78	6.6
3.14	14	7.13	3.66	1.992	7.46
12.7	5.63	6.4	4.14	5.3	0.784
48.2	20.1	6.21	6.26	8.12	4.06
3.09	4.21	3.14	3.04	0.656	3.04
26.7	14.8	11.8	23.8	30.6	12.14
1.31	1.31	1.8	10.86	3.38	5.44
30.2	29.3	12.5	17.36	9.06	11.88
1.87	37.5	3.86	7.94	0.402	4.48
2.72	1.87	1.87	3.66	4.06	19.48
36.9	81.3	25.1	38.8	59.4	16.8
24	13.5	24.2	20.6	14.86	6.66
3.14	1.87	0.772	2.67	1.87	0.289
11.78	8.32	17.68	4.47	4.47	5.88
2.62	1.834	0.1644	0.489	1.87	4.88
1.092	1.91	1.782	1.71	4.01	4.31
37.2	19.54	11.26	8.05	7.18	29
0.214	7.58	10.5	5.38	25.9	26.9
1.87	2.22	0.64	0.887	1.87	1.87
2.88	3.7	0.64	1.15	1.96	9.43
9.06	17.8	1.91	5.77	13.3	7.85
29.6	10.54	12.6	9.77	0.108	6.86
10.1	3.42	14.58	3.18	1.83	3.3
6.5	1.962	5.42	3.37	5.04	2.21
2.56	3.24	1.284	1.87	1.87	6.69
0.581	2.14	3.77	7.3	4.82	1.28
1.51	1.51	0.799	1.83	0.48	7.3
11.2	3.56	29	3.42	15.9	4.99
8.86	3.09	26.1	1.67	3.83	2.67
1.51	2.29	1.45	1.67	3.96	1.67
END

Results

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	rho	0.3026	0.04594	2.218E-4	0.2043	0.3056	0.3836	4001	50000
	mu	1.573	0.1067	5.283E-4	1.364	1.573	1.784	4001	50000
	sigma	1.18	0.06546	3.195E-4	1.061	1.177	1.318	4001	50000
	tau	0.7251	0.07985	3.842E-4	0.576	0.7224	0.8887	4001	50000
	
		Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
logy	552.319	549.422	2.896	555.215	
total	552.319	549.422	2.896	555.215	

 

1i. MVN specification with ARMA(1,1) covariance matrix.
   Biomarker (y=IL-1beta) measured from saliva at 6 time points for 30 healthy 
   subjects. 

model{ 
for( j in 1:30){
     for(i in 1:6){
           logy[ j, i ] <-  log(y[ j, i ])
     }
  }
  for( j in 1:30){
    logy[ j, 1:6 ] ~ dmnorm(m[1:6], precision[1:6,1:6])
  }
  for(j in 1:6){
	for(k in 1:6){
	  covariance[j,k] <-  sigma2*equals(j , k) + sigma2*gamma*equals(abs(j-k) , 1)  
                       +sigma2*gamma*pow(phi,abs(j-k)-1)*step(abs(j-k)-2)
	}
  }
 for(i in 1:6){   m[ i ] <- mu  }
 precision[1:6,1:6] <- inverse(covariance[1:6,1:6])
 sigma ~dunif(0,100)
 sigma2 <- sigma*sigma
 tau <- 1/sigma2
 mu ~ dnorm(0,0.001)
 theta ~ dunif(-1,1)
 phi ~ dunif(-1,1)
 gamma <- (phi-theta)*(1-phi*theta)  /  (  (1-phi*phi) + (phi-theta)*(phi-theta)  )
 logyf ~ dnorm(mu, tau)
 yf <- exp(logyf)
}

Inits   
list(mu =0, sigma = 1, theta=0, phi=0, logyf=0)

IL-1beta Data

y[,1]	y[,2]    y[,3]     y[,4]   y[,5]   y[,6]  
22.8	17.6	11.6	7	8.78	6.6
3.14	14	7.13	3.66	1.992	7.46
12.7	5.63	6.4	4.14	5.3	0.784
48.2	20.1	6.21	6.26	8.12	4.06
3.09	4.21	3.14	3.04	0.656	3.04
26.7	14.8	11.8	23.8	30.6	12.14
1.31	1.31	1.8	10.86	3.38	5.44
30.2	29.3	12.5	17.36	9.06	11.88
1.87	37.5	3.86	7.94	0.402	4.48
2.72	1.87	1.87	3.66	4.06	19.48
36.9	81.3	25.1	38.8	59.4	16.8
24	13.5	24.2	20.6	14.86	6.66
3.14	1.87	0.772	2.67	1.87	0.289
11.78	8.32	17.68	4.47	4.47	5.88
2.62	1.834	0.1644	0.489	1.87	4.88
1.092	1.91	1.782	1.71	4.01	4.31
37.2	19.54	11.26	8.05	7.18	29
0.214	7.58	10.5	5.38	25.9	26.9
1.87	2.22	0.64	0.887	1.87	1.87
2.88	3.7	0.64	1.15	1.96	9.43
9.06	17.8	1.91	5.77	13.3	7.85
29.6	10.54	12.6	9.77	0.108	6.86
10.1	3.42	14.58	3.18	1.83	3.3
6.5	1.962	5.42	3.37	5.04	2.21
2.56	3.24	1.284	1.87	1.87	6.69
0.581	2.14	3.77	7.3	4.82	1.28
1.51	1.51	0.799	1.83	0.48	7.3
11.2	3.56	29	3.42	15.9	4.99
8.86	3.09	26.1	1.67	3.83	2.67
1.51	2.29	1.45	1.67	3.96	1.67
END

Results

node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	gamma	0.5064	0.07811	5.984E-4	0.3516	0.5076	0.6569	4001	50000
	mu	1.582	0.1593	7.082E-4	1.269	1.582	1.896	4001	50000
	phi	0.8418	0.08065	4.634E-4	0.6541	0.8527	0.9673	4001	50000
	sigma	1.219	0.09223	7.655E-4	1.067	1.209	1.424	4001	50000
	tau	0.6842	0.09848	7.666E-4	0.4931	0.6839	0.8791	4001	50000
	
	Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
logy	525.090	521.374	3.716	528.807	
total	525.090	521.374	3.716	528.807	

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	gamma	0.507	0.07623	9.276E-4	0.3552	0.5079	0.6529	3999	20000
	logyf	1.59	1.234	0.008016	-0.8055	1.584	4.011	3999	20000
	mu	1.582	0.1586	0.001188	1.269	1.582	1.892	3999	20000
	phi	0.8412	0.08025	7.17E-4	0.6563	0.8522	0.9657	3999	20000
	sigma	1.218	0.09007	0.001106	1.065	1.21	1.417	3999	20000
	tau	0.6846	0.09778	0.001162	0.4983	0.6836	0.8819	3999	20000
	yf	10.72	22.18	0.1588	0.4469	4.875	55.18	3999	20000

 # this set of results includes inference for yf

1j.  A MVN specification with unstructured covariance matrix.
    Biomarker (y=IL-1beta) measured from saliva at 6 time points for 30 healthy 
    subjects.

model{ 
for( j in 1:30){
     for(i in 1:6){
           logy[ j, i ] <-  log(y[ j, i ])
     }
  }
  for( j in 1:30){
    logy[ j, 1:6 ] ~ dmnorm(m[1:6], precision[1:6,1:6])
  }
  for(i in 1:6){   m[ i ] <- mu  }
 precision[1:6,1:6] ~ dwish(R[1:6 , 1:6 ], 6)
 mu ~ dnorm(0,0.001)
 #covariance[1:6,1:6] <- inverse(precision[1:6,1:6])
}

Inits   
list(mu =0,  precision=structure(.Data=c(.001,0,0,0,0,0,
                                        0,.001,0,0,0,0,
                                        0,0,.001,0,0,0,
                                        0,0,0,.001,0,0,
                                        0,0,0,0,.001,0,
                                        0,0,0,0,0,.001),.Dim=c(6,6))))
IL-1beta Data

list(R=structure(.Data=c(.01,0,0,0,0,0,
                         0,.01,0,0,0,0,
                         0,0,.01,0,0,0,
                         0,0,0,.01,0,0,
                         0,0,0,0,.01,0,
                         0,0,0,0,0,.01),.Dim=c(6,6)))  )
y[,1]	y[,2]    y[,3]     y[,4]   y[,5]   y[,6]  
22.8	17.6	11.6	7	8.78	6.6
3.14	14	7.13	3.66	1.992	7.46
12.7	5.63	6.4	4.14	5.3	0.784
48.2	20.1	6.21	6.26	8.12	4.06
3.09	4.21	3.14	3.04	0.656	3.04
26.7	14.8	11.8	23.8	30.6	12.14
1.31	1.31	1.8	10.86	3.38	5.44
30.2	29.3	12.5	17.36	9.06	11.88
1.87	37.5	3.86	7.94	0.402	4.48
2.72	1.87	1.87	3.66	4.06	19.48
36.9	81.3	25.1	38.8	59.4	16.8
24	13.5	24.2	20.6	14.86	6.66
3.14	1.87	0.772	2.67	1.87	0.289
11.78	8.32	17.68	4.47	4.47	5.88
2.62	1.834	0.1644	0.489	1.87	4.88
1.092	1.91	1.782	1.71	4.01	4.31
37.2	19.54	11.26	8.05	7.18	29
0.214	7.58	10.5	5.38	25.9	26.9
1.87	2.22	0.64	0.887	1.87	1.87
2.88	3.7	0.64	1.15	1.96	9.43
9.06	17.8	1.91	5.77	13.3	7.85
29.6	10.54	12.6	9.77	0.108	6.86
10.1	3.42	14.58	3.18	1.83	3.3
6.5	1.962	5.42	3.37	5.04	2.21
2.56	3.24	1.284	1.87	1.87	6.69
0.581	2.14	3.77	7.3	4.82	1.28
1.51	1.51	0.799	1.83	0.48	7.3
11.2	3.56	29	3.42	15.9	4.99
8.86	3.09	26.1	1.67	3.83	2.67
1.51	2.29	1.45	1.67	3.96	1.67
END

Results

	 node			 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	covariance[1,1]	1.982	0.5408	0.002312	1.185	1.897	3.285	4001	50000
	covariance[1,2]	1.018	0.3703	0.001634	0.4544	0.9648	1.899	4001	50000
	covariance[1,3]	1.079	0.4166	0.001812	0.4324	1.02	2.066	4001	50000
	covariance[1,4]	0.7095	0.3033	0.001307	0.2294	0.6699	1.415	4001	50000
	covariance[1,5]	0.4737	0.3752	0.001726	-0.1972	0.4449	1.304	4001	50000
	covariance[1,6]	0.2926	0.2812	0.001262	-0.2162	0.2739	0.908	4001	50000
	covariance[2,1]	1.018	0.3703	0.001634	0.4544	0.9648	1.899	4001	50000
	covariance[2,2]	1.327	0.3679	0.001584	0.7906	1.265	2.209	4001	50000
	covariance[2,3]	0.872	0.341	0.001484	0.3422	0.8229	1.678	4001	50000
	covariance[2,4]	0.7927	0.2714	0.001132	0.3853	0.7504	1.435	4001	50000
	covariance[2,5]	0.4959	0.3143	0.001356	-0.04068	0.4666	1.205	4001	50000
	covariance[2,6]	0.5577	0.2531	0.001065	0.1555	0.5251	1.147	4001	50000
	covariance[3,1]	1.079	0.4166	0.001812	0.4324	1.02	2.066	4001	50000
	covariance[3,2]	0.872	0.341	0.001484	0.3422	0.8229	1.678	4001	50000
	covariance[3,3]	1.836	0.4975	0.00202	1.103	1.757	3.039	4001	50000
	covariance[3,4]	0.9326	0.3185	0.001366	0.4537	0.885	1.687	4001	50000
	covariance[3,5]	0.7921	0.3828	0.001695	0.1681	0.7458	1.679	4001	50000
	covariance[3,6]	0.3925	0.2734	0.001277	-0.08634	0.3686	1.003	4001	50000
	covariance[4,1]	0.7095	0.3033	0.001307	0.2294	0.6699	1.415	4001	50000
	covariance[4,2]	0.7927	0.2714	0.001132	0.3853	0.7504	1.435	4001	50000
	covariance[4,3]	0.9326	0.3185	0.001366	0.4537	0.885	1.687	4001	50000
	covariance[4,4]	1.027	0.2824	0.001192	0.6119	0.9813	1.705	4001	50000
	covariance[4,5]	0.6318	0.2931	0.001307	0.1631	0.5949	1.313	4001	50000
	covariance[4,6]	0.3864	0.2113	9.07E-4	0.0347	0.364	0.8669	4001	50000
	covariance[5,1]	0.4737	0.3752	0.001726	-0.1972	0.4449	1.304	4001	50000
	covariance[5,2]	0.4959	0.3143	0.001356	-0.04068	0.4666	1.205	4001	50000
	covariance[5,3]	0.7921	0.3828	0.001695	0.1681	0.7458	1.679	4001	50000
	covariance[5,4]	0.6318	0.2931	0.001307	0.1631	0.5949	1.313	4001	50000
	covariance[5,5]	1.866	0.5143	0.002291	1.11	1.781	3.095	4001	50000
	covariance[5,6]	0.4096	0.2784	0.001306	-0.06993	0.3815	1.038	4001	50000
	covariance[6,1]	0.2926	0.2812	0.001262	-0.2162	0.2739	0.908	4001	50000
	covariance[6,2]	0.5577	0.2531	0.001065	0.1555	0.5251	1.147	4001	50000
	covariance[6,3]	0.3925	0.2734	0.001277	-0.08634	0.3686	1.003	4001	50000
	covariance[6,4]	0.3864	0.2113	9.07E-4	0.0347	0.364	0.8669	4001	50000
	covariance[6,5]	0.4096	0.2784	0.001306	-0.06993	0.3815	1.038	4001	50000
	covariance[6,6]	1.039	0.2841	0.001242	0.6235	0.992	1.725	4001	50000
	mu	1.54	0.1471	7.517E-4	1.251	1.541	1.83	4001	50000

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
			Dbar		Dhat		pD		DIC	
logy	519.158	499.550	19.608	538.766	
total	519.158	499.550	19.608	538.766

	

Modeling a Wishart

model{
 for(i in 1:a){  
   b[i,1:2] ~ dmnorm(m[1:2], taub[1:2,1:2])  
 }
 taub[1:2,1:2] ~ dwish(B[1:2,1:2], 3)
}
list(a=1000, m=c(0,0))
B[,1] B[,2]    
0.001  0       
0      0.001   
END

#################################
## Hierarchical models: Meta-analysis ##
## Example 10.2.3 and Exercise 10.10 ##
##################################

2.  Model, Data, and Inits

model{
 for(i in 1:m){
   LOR[i] ~ dnorm(eta[i] , tau[i])
   eta[i] ~ dnorm(mu, taub)
   tau[i] <- pow(1/Se[i],2)
 }
 mu ~ dnorm(0, 0.001)
 taub <- 1/(sigmab*sigmab)
 sigmab ~ dunif(0,100)
 # Median of OR distribution
 medianOR <- exp(mu)
 prob <- step(medianOR - 1)
 # 10th and 90th percentiles of distribution of ORs
 ninetypct <- exp(mu + 1.28*sigmab)
 tenpct <- exp(mu - 1.28*sigmab)
# Predictive density of ORs
 muf ~ dnorm(mu, taub)
 ORf <- exp(muf)
 for(i in 1:m){ OR[i] <- exp(eta[i]) }
}
list(m=14,
LOR=c(-0.79493,-0.72392,-0.63227,-0.45604,
      -0.40547,-0.3071,-0.30157,-0.28041,
      -0.15861,0.028988,0.068993,0.100155,
      0.145954,0.200933),
Se=c(0.436627,0.384353,0.332176,0.431258,
     0.629153,0.246629,0.302944,0.29574,
     0.444878,0.448152,0.440058,0.277911,
     0.241851,0.240026))
list(mu=0,sigmab=1)

############################
##  The dental data             ##
############################

3a.  Mixed model with random intercepts and random slopes
        Response:  Distance at ages 8, 10, 12, and 14
        Covariate:  Sex, Age
         Example 10.2.4

Model and Inits

model{ 
 for( i in 1:27){
   for(j in 1:4){
    Distance[ i , j ] ~ dnorm(m[ i , j ], tau)
    m[ i , j  ]  <- (beta[1] + b[ i , 1 ]) + (beta[2] + b[ i , 2 ])*Age[ j ]  + beta[3]*Male[ i ]  
    } 
}
 for( i in 1:27){
    b[ i , 1:2 ] ~ dmnorm(mu0[1:2], prec[1:2,1:2])
 }
 prec[1:2,1:2] ~ dwish(R[ 1:2 , 1:2 ],  2)
 for(i in 1:3){   beta[ i ] ~ dnorm(0,0.0001)  }
 tau ~dgamma(0.0001,0.0001)
 sigma <- sqrt(1/(tau*tau))
}
Inits
list(beta =c(0,0,0), tau=1, prec=structure(.Data=c(.001,0,0,.001),.Dim=c(2,2))))

Dental Data

list(R=structure(.Data=c(.001,0,0,.001),.Dim=c(2,2)), mu0=c(0,0)  ,
Male=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
               0,0,0,0,0,0,0,0,0,0,0),Age=c(8,10,12,14),
Distance=structure(.Data=c(
26,	25,	29,	31,
21.5,	22.5,	23,	26.5,
23,	22.5,	24,	27.5,
25.5,	27.5,	26.5,	27,
20,	23.5,	22.5,	26,
24.5,	25.5,	27,	28.5,
22,	22,	24.5,	26.5,
24,	21.5,	24.5,	25.5,
23,	20.5,	31,	26,
27.5,	28,	31,	31.5,
23,	23,	23.5,	25,
21.5,	23.5,	24,	28,
17,	24.5,	26,	29.5,
22.5,	25.5,	25.5,	26,
23,	24.5,	26,	30,
22,	21.5,	23.5,	25,
21,	20,	21.5,	23,
21,	21.5,	24,	25.5,
20.5,	24,	24.5,	26,
23.5,	24.5,	25,	26.5,
21.5,	23,	22.5,	23.5,
20, 21,	21,	22.5,
21.5,	22.5,	23,	25,
23,	23,	23.5,	24,
20,	21,	22,	21.5,
16.5,	19,	19,	19.5,
24.5,	25,	28,	28),.Dim=c(27,4)) )

Results

 DIC = 412.8
 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	beta[1]	15.92	0.8536	0.002906	14.24	15.92	17.6	501001	509000
	beta[2]	0.6601	0.07046	2.03E-4	0.5225	0.66	0.7984	501001	509000
	beta[3]	1.335	0.7644	0.004283	-0.2156	1.352	2.793	501001	509000
	m[1,1]	24.62	0.5645	0.00361	23.52	24.61	25.74	501001	509000
	m[1,2]	26.43	0.6305	0.001734	25.18	26.43	27.66	501001	509000
	m[1,3]	28.23	0.7407	0.001271	26.77	28.24	29.68	501001	509000
	m[1,4]	30.04	0.8785	0.002976	28.31	30.04	31.77	501001	509000
	m[2,1]	21.86	0.5513	0.001593	20.78	21.86	22.94	501001	509000
	m[2,2]	23.03	0.6213	0.001191	21.8	23.03	24.25	501001	509000
	m[2,3]	24.19	0.7278	0.001082	22.76	24.19	25.62	501001	509000
	m[2,4]	25.35	0.8575	0.001341	23.67	25.35	27.04	501001	509000
	m[3,1]	22.4	0.5473	9.424E-4	21.33	22.4	23.47	501001	509000
	m[3,2]	23.69	0.6176	9.209E-4	22.47	23.69	24.9	501001	509000
	m[3,3]	24.97	0.724	0.001016	23.56	24.97	26.39	501001	509000
	m[3,4]	26.26	0.8531	0.001201	24.59	26.26	27.93	501001	509000
	m[4,1]	23.71	0.5554	0.003092	22.63	23.71	24.81	501001	509000
	m[4,2]	25.31	0.6228	0.002102	24.09	25.31	26.53	501001	509000
	m[4,3]	26.9	0.7293	0.001342	25.47	26.91	28.33	501001	509000
	m[4,4]	28.5	0.8605	0.00131	26.81	28.5	30.19	501001	509000
	m[5,1]	21.65	0.5537	0.001951	20.56	21.65	22.74	501001	509000
	m[5,2]	22.76	0.6226	0.001333	21.54	22.76	23.99	501001	509000
	m[5,3]	23.87	0.7285	0.00107	22.45	23.87	25.31	501001	509000
	m[5,4]	24.99	0.8579	0.001381	23.31	24.99	26.68	501001	509000
	m[6,1]	23.69	0.5525	0.002442	22.61	23.69	24.79	501001	509000
	m[6,2]	25.28	0.6218	0.001483	24.06	25.28	26.5	501001	509000
	m[6,3]	26.88	0.73	0.001085	25.44	26.88	28.31	501001	509000
	m[6,4]	28.47	0.8625	0.001703	26.77	28.47	30.16	501001	509000
	m[7,1]	22.1	0.5498	0.001224	21.02	22.1	23.18	501001	509000
	m[7,2]	23.32	0.6203	0.001023	22.1	23.32	24.54	501001	509000
	m[7,3]	24.54	0.727	0.001019	23.11	24.54	25.97	501001	509000
	m[7,4]	25.76	0.8565	0.001213	24.08	25.76	27.44	501001	509000
	m[8,1]	22.06	0.5493	9.709E-4	20.99	22.06	23.14	501001	509000
	m[8,2]	23.27	0.6202	9.044E-4	22.06	23.27	24.49	501001	509000
	m[8,3]	24.48	0.7273	0.001136	23.06	24.48	25.92	501001	509000
	m[8,4]	25.69	0.8572	0.001536	24.01	25.69	27.38	501001	509000
	m[9,1]	23.0	0.547	0.001053	21.93	23.0	24.07	501001	509000
	m[9,2]	24.43	0.6182	8.785E-4	23.21	24.43	25.64	501001	509000
	m[9,3]	25.86	0.726	0.001104	24.43	25.86	27.28	501001	509000
	m[9,4]	27.29	0.8567	0.001566	25.6	27.29	28.97	501001	509000
	m[10,1]	25.64	0.5914	0.005765	24.5	25.64	26.83	501001	509000
	m[10,2]	27.69	0.6485	0.002929	26.41	27.7	28.96	501001	509000
	m[10,3]	29.74	0.7586	0.001352	28.24	29.75	31.22	501001	509000
	m[10,4]	31.79	0.9024	0.003602	30.01	31.8	33.55	501001	509000
	m[11,1]	21.9	0.5506	0.001091	20.81	21.89	22.98	501001	509000
	m[11,2]	23.07	0.6214	9.409E-4	21.85	23.06	24.29	501001	509000
	m[11,3]	24.24	0.7286	0.001245	22.81	24.23	25.67	501001	509000
	m[11,4]	25.41	0.8587	0.001785	23.72	25.41	27.1	501001	509000
	m[12,1]	22.46	0.5472	9.879E-4	21.39	22.46	23.54	501001	509000
	m[12,2]	23.77	0.6179	9.866E-4	22.55	23.77	24.98	501001	509000
	m[12,3]	25.07	0.7247	0.001094	23.65	25.07	26.5	501001	509000
	m[12,4]	26.38	0.8541	0.001282	24.7	26.38	28.05	501001	509000
	m[13,1]	22.72	0.5487	0.001546	21.63	22.72	23.79	501001	509000
	m[13,2]	24.08	0.6201	0.001868	22.86	24.08	25.29	501001	509000
	m[13,3]	25.45	0.728	0.00229	24.02	25.45	26.87	501001	509000
	m[13,4]	26.81	0.8586	0.002767	25.13	26.81	28.5	501001	509000
	m[14,1]	22.72	0.5463	0.001078	21.65	22.72	23.8	501001	509000
	m[14,2]	24.09	0.6172	9.756E-4	22.88	24.09	25.3	501001	509000
	m[14,3]	25.45	0.7244	0.001024	24.03	25.45	26.87	501001	509000
	m[14,4]	26.82	0.8543	0.001204	25.14	26.82	28.49	501001	509000
	m[15,1]	23.5	0.5502	0.001543	22.43	23.5	24.59	501001	509000
	m[15,2]	25.05	0.6209	9.22E-4	23.83	25.05	26.27	501001	509000
	m[15,3]	26.6	0.7296	0.00129	25.17	26.6	28.03	501001	509000
	m[15,4]	28.15	0.8619	0.002196	26.46	28.15	29.84	501001	509000
	m[16,1]	21.57	0.5528	0.001733	20.48	21.57	22.65	501001	509000
	m[16,2]	22.66	0.622	0.001101	21.45	22.66	23.89	501001	509000
	m[16,3]	23.76	0.7286	0.001091	22.33	23.76	25.19	501001	509000
	m[16,4]	24.85	0.8587	0.001714	23.16	24.85	26.54	501001	509000
	m[17,1]	20.52	0.5699	0.003331	19.39	20.52	21.64	501001	509000
	m[17,2]	21.37	0.6348	0.00166	20.12	21.36	22.62	501001	509000
	m[17,3]	22.21	0.7418	0.001308	20.76	22.21	23.68	501001	509000
	m[17,4]	23.06	0.8756	0.002821	21.34	23.06	24.79	501001	509000
	m[18,1]	21.12	0.566	0.001239	20.01	21.12	22.23	501001	509000
	m[18,2]	22.43	0.6243	0.00103	21.2	22.43	23.65	501001	509000
	m[18,3]	23.73	0.726	0.00102	22.3	23.73	25.16	501001	509000
	m[18,4]	25.03	0.8556	0.001213	23.35	25.03	26.72	501001	509000
	m[19,1]	21.6	0.5692	0.001149	20.48	21.6	22.72	501001	509000
	m[19,2]	23.02	0.627	9.822E-4	21.78	23.02	24.24	501001	509000
	m[19,3]	24.43	0.728	0.001088	23.0	24.43	25.86	501001	509000
	m[19,4]	25.84	0.8569	0.001407	24.16	25.84	27.53	501001	509000
	m[20,1]	22.2	0.5758	0.002134	21.07	22.2	23.33	501001	509000
	m[20,2]	23.75	0.6315	0.001407	22.51	23.75	24.99	501001	509000
	m[20,3]	25.3	0.7313	0.001051	23.86	25.3	26.73	501001	509000
	m[20,4]	26.85	0.8599	0.001388	25.16	26.85	28.54	501001	509000
	m[21,1]	20.75	0.5661	0.001221	19.63	20.75	21.86	501001	509000
	m[21,2]	21.97	0.6246	9.223E-4	20.74	21.97	23.19	501001	509000
	m[21,3]	23.18	0.7269	0.001147	21.75	23.18	24.61	501001	509000
	m[21,4]	24.39	0.8574	0.001699	22.71	24.39	26.08	501001	509000
	m[22,1]	19.85	0.5711	0.002617	18.72	19.85	20.97	501001	509000
	m[22,2]	20.85	0.6284	0.001371	19.62	20.85	22.09	501001	509000
	m[22,3]	21.85	0.7328	0.001134	20.42	21.85	23.3	501001	509000
	m[22,4]	22.85	0.8675	0.002252	21.15	22.85	24.56	501001	509000
	m[23,1]	21.06	0.5656	0.001145	19.94	21.06	22.17	501001	509000
	m[23,2]	22.34	0.6239	9.44E-4	21.11	22.34	23.56	501001	509000
	m[23,3]	23.63	0.7256	0.001011	22.2	23.63	25.05	501001	509000
	m[23,4]	24.91	0.8553	0.001307	23.23	24.91	26.59	501001	509000
	m[24,1]	21.19	0.567	0.001128	20.08	21.19	22.3	501001	509000
	m[24,2]	22.51	0.6253	0.00108	21.28	22.51	23.73	501001	509000
	m[24,3]	23.82	0.7267	0.001259	22.39	23.82	25.25	501001	509000
	m[24,4]	25.14	0.8561	0.001589	23.45	25.14	26.82	501001	509000
	m[25,1]	19.82	0.5698	0.002544	18.7	19.83	20.94	501001	509000
	m[25,2]	20.82	0.6271	0.001284	19.59	20.82	22.06	501001	509000
	m[25,3]	21.82	0.7319	0.00119	20.38	21.81	23.26	501001	509000
	m[25,4]	22.81	0.8669	0.002404	21.11	22.81	24.52	501001	509000
	m[26,1]	18.25	0.6	0.005512	17.05	18.25	19.41	501001	509000
	m[26,2]	18.87	0.6479	0.002712	17.61	18.87	20.16	501001	509000
	m[26,3]	19.5	0.7548	0.001285	18.03	19.49	21.0	501001	509000
	m[26,4]	20.12	0.9	0.003642	18.37	20.12	21.91	501001	509000
	m[27,1]	23.18	0.5933	0.003522	22.01	23.18	24.35	501001	509000
	m[27,2]	24.96	0.6448	0.001903	23.69	24.96	26.22	501001	509000
	m[27,3]	26.74	0.7438	0.001166	25.27	26.75	28.19	501001	509000
	m[27,4]	28.53	0.8745	0.002369	26.8	28.53	30.24	501001	509000
	sigma	2.117	0.3495	8.701E-4	1.54	2.08	2.904	501001	509000
	tau	0.4849	0.07797	1.94E-4	0.3444	0.4807	0.6492	501001	509000

3b.  MVN regression with unstructured covariance.
        Response:  Distance at ages 8, 10, 12, and 14
        Covariate:  Sex, Age
        Example 10.4.1 and Table 10.5

Model 

model{ 
  for(i in 1:16){
    Distance[ i, 1:4 ] ~ dmnorm(m[ i, 1:4], precision[1:4,1:4])
    for(j in 1:4){  m[ i, j]  <- (beta[1] + beta[3]) + (beta[2] + beta[4])*Age[ j ] }   ## Male
  }
for( i in 17:27){
    Distance[ i, 1:4 ] ~ dmnorm(m[ i, 1:4], precision[1:4,1:4])
    for(j in 1:4){  m[ i, j]  <- beta[1] + beta[2]*Age[ j ]    }     ## Female
  }
covariance[1:4,1:4] <- inverse(precision[1:4,1:4])
 precision[1:4,1:4] ~ dwish(R[1:4 , 1:4 ], 4)
 for(i in 1:4){   beta[ i ] ~ dnorm(0,0.0001)  }
  for(i in 1:4){ 
    meanB[i] <- beta[1] + beta[3] + (beta[2] + beta[4])*Age[i]
    meanG[i] <- beta[1] + beta[2]*Age[i]    
    diff[i] <- meanB[i] - meanG[i]
    testdiff[i] <- step(diff[i])
  }
 rateB <- beta[2]+beta[4]
 rateG <- beta[2]
 test <- step(beta[4])
}

Inits
list(beta =c(0,0,0), precision=structure(.Data=c(.001,0,0,0,
                                        0,.001,0,0,
                                        0,0,.001,0,
                                        0,0,0,.001),.Dim=c(4,4))))

Dental Data
list(Age=c(8,10,12,14),
R=structure(.Data=c(.001,0,0,0,
                     0,.001,0,0,
                     0,0,.001,0,
                     0,0,0,.001),.Dim=c(4,4))  , 
Distance=structure(.Data=c(
26,	25,	29,	31,
21.5,	22.5,	23,	26.5,
23,	22.5,	24,	27.5,
25.5,	27.5,	26.5,	27,
20,	23.5,	22.5,	26,
24.5,	25.5,	27,	28.5,
22,	22,	24.5,	26.5,
24,	21.5,	24.5,	25.5,
23,	20.5,	31,	26,
27.5,	28,	31,	31.5,
23,	23,	23.5,	25,
21.5,	23.5,	24,	28,
17,	24.5,	26,	29.5,
22.5,	25.5,	25.5,	26,
23,	24.5,	26,	30,
22,	21.5,	23.5,	25,
21,	20,	21.5,	23,
21,	21.5,	24,	25.5,
20.5,	24,	24.5,	26,
23.5,	24.5,	25,	26.5,
21.5,	23,	22.5,	23.5,
20, 21,	21,	22.5,
21.5,	22.5,	23,	25,
23,	23,	23.5,	24,
20,	21,	22,	21.5,
16.5,	19,	19,	19.5,
24.5,	25,	28,	28),.Dim=c(27,4)) )

Results

From unstructured covariance matrix;  different lines for boys than girls
      node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	beta[1]	17.42	1.156	0.006595	15.13	17.41	19.7	50001	950000
	beta[2]	0.4771	0.09769	5.6E-4	0.2833	0.4772	0.6694	50001	950000
	beta[3]	-1.57	1.541	0.008905	-4.622	-1.565	1.464	50001	950000
	beta[4]	0.3495	0.1303	7.557E-4	0.09312	0.3491	0.6079	50001	950000
	covariance[1,1]	5.637	1.69	0.001959	3.237	5.341	9.743	50001	950000
	covariance[1,2]	2.817	1.199	0.001313	1.022	2.632	5.681	50001	950000
	covariance[1,3]	3.992	1.507	0.001625	1.77	3.75	7.603	50001	950000
	covariance[1,4]	2.823	1.283	0.00142	0.8605	2.635	5.869	50001	950000
	covariance[2,1]	2.817	1.199	0.001313	1.022	2.632	5.681	50001	950000
	covariance[2,2]	4.355	1.309	0.001486	2.496	4.126	7.54	50001	950000
	covariance[2,3]	3.093	1.289	0.001438	1.168	2.892	6.174	50001	950000
	covariance[2,4]	3.445	1.24	0.001422	1.648	3.237	6.451	50001	950000
	covariance[3,1]	3.992	1.507	0.001625	1.77	3.75	7.603	50001	950000
	covariance[3,2]	3.093	1.289	0.001438	1.168	2.892	6.174	50001	950000
	covariance[3,3]	6.509	1.911	0.002091	3.773	6.181	11.14	50001	950000
	covariance[3,4]	4.298	1.516	0.001704	2.106	4.044	7.969	50001	950000
	covariance[4,1]	2.823	1.283	0.00142	0.8605	2.635	5.869	50001	950000
	covariance[4,2]	3.445	1.24	0.001422	1.648	3.237	6.451	50001	950000
	covariance[4,3]	4.298	1.516	0.001704	2.106	4.044	7.969	50001	950000
	covariance[4,4]	5.186	1.561	0.001808	2.969	4.914	8.986	50001	950000
	diff[1]	1.226	0.803	0.00302	-0.3583	1.226	2.811	50001	950000
	diff[2]	1.925	0.743	0.001731	0.46	1.925	3.393	50001	950000
	diff[3]	2.624	0.7715	0.0012	1.102	2.623	4.147	50001	950000
	diff[4]	3.323	0.8799	0.00211	1.588	3.322	5.059	50001	950000
	meanB[1]	22.46	0.5263	0.001071	21.42	22.46	23.5	50001	950000
	meanB[2]	24.11	0.4877	6.354E-4	23.15	24.11	25.08	50001	950000
	meanB[3]	25.76	0.5068	4.618E-4	24.76	25.76	26.77	50001	950000
	meanB[4]	27.42	0.5777	7.593E-4	26.27	27.42	28.56	50001	950000
	meanG[1]	21.23	0.6022	0.002232	20.04	21.23	22.42	50001	950000
	meanG[2]	22.19	0.557	0.001274	21.08	22.19	23.28	50001	950000
	meanG[3]	23.14	0.5781	8.809E-4	22.0	23.14	24.28	50001	950000
	meanG[4]	24.09	0.6592	0.001561	22.79	24.1	25.39	50001	950000
	rateB	0.8266	0.08519	2.631E-4	0.6583	0.8266	0.9952	50001	950000
	rateG	0.4771	0.09769	5.6E-4	0.2833	0.4772	0.6694	50001	950000
	test	0.9953	0.06811	2.211E-4	1.0	1.0	1.0	50001	950000
	testdiff[1]	0.9377	0.2416	6.01E-4	0.0	1.0	1.0	50001	950000
	testdiff[2]	0.9941	0.07673	1.066E-4	1.0	1.0	1.0	50001	950000
	testdiff[3]	0.9993	0.02682	3.263E-5	1.0	1.0	1.0	50001	950000
	testdiff[4]	0.9997	0.01616	1.978E-5	1.0	1.0	1.0	50001	950000

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
				Dbar	     Dhat	      pD	     DIC	
Distance	433.498	420.188	13.310	446.808	
total		  433.498	420.188	13.310	446.808	

From unstructured covariance matrix;  different intercepts, same slopes for boys than girls
	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	beta[1]	15.38	0.9462	0.003382	13.47	15.39	17.2	50001	950000
	beta[2]	0.6744	0.07106	2.412E-4	0.5339	0.6745	0.8147	50001	950000
	beta[3]	2.045	0.823	0.001904	0.4213	2.043	3.67	50001	950000
	covariance[1,1]	5.621	1.683	0.002011	3.226	5.329	9.699	50001	950000
	covariance[1,2]	2.934	1.238	0.001477	1.105	2.737	5.897	50001	950000
	covariance[1,3]	3.994	1.521	0.001666	1.754	3.752	7.642	50001	950000
	covariance[1,4]	2.769	1.335	0.001613	0.7238	2.574	5.938	50001	950000
	covariance[2,1]	2.934	1.238	0.001477	1.105	2.737	5.897	50001	950000
	covariance[2,2]	4.418	1.347	0.001558	2.514	4.182	7.699	50001	950000
	covariance[2,3]	3.061	1.308	0.001453	1.104	2.86	6.191	50001	950000
	covariance[2,4]	3.329	1.281	0.001568	1.464	3.117	6.412	50001	950000
	covariance[3,1]	3.994	1.521	0.001666	1.754	3.752	7.642	50001	950000
	covariance[3,2]	3.061	1.308	0.001453	1.104	2.86	6.191	50001	950000
	covariance[3,3]	6.622	1.965	0.00214	3.815	6.283	11.38	50001	950000
	covariance[3,4]	4.509	1.646	0.002087	2.16	4.223	8.511	50001	950000
	covariance[4,1]	2.769	1.335	0.001613	0.7238	2.574	5.938	50001	950000
	covariance[4,2]	3.329	1.281	0.001568	1.464	3.117	6.412	50001	950000
	covariance[4,3]	4.509	1.646	0.002087	2.16	4.223	8.511	50001	950000
	covariance[4,4]	5.625	1.783	0.002489	3.131	5.301	9.989	50001	950000
	meanB[1]	22.82	0.5384	8.855E-4	21.77	22.81	23.9	50001	950000
	meanB[2]	24.16	0.5123	7.09E-4	23.16	24.16	25.18	50001	950000
	meanB[3]	25.51	0.5249	8.288E-4	24.46	25.52	26.54	50001	950000
	meanB[4]	26.86	0.5736	0.001156	25.7	26.87	27.97	50001	950000
	meanG[1]	20.77	0.6162	0.001707	19.53	20.78	21.96	50001	950000
	meanG[2]	22.12	0.5931	0.001395	20.94	22.12	23.29	50001	950000
	meanG[3]	23.47	0.6035	0.001202	22.29	23.47	24.67	50001	950000
	meanG[4]	24.82	0.6458	0.001187	23.57	24.81	26.12	50001	950000

Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
				Dbar		Dhat		pD		DIC	
Distance	439.185	426.968	12.217	451.401	
total		439.185	426.968	12.217	  451.401