using Example R: Heart Valves Study
Goal: Show that the thrombogenicity rate (TR) is less than two times the objective performance criterion
R and WinBUGS Examples p. 1/4
Example using R: Heart Valves Study
Goal: Show that the thrombogenicity rate (TR) is less than two times the objective performance criterion Data: From both the current study and a previous study on a similar product (St. Jude mechanical valve).
R and WinBUGS Examples p. 1/4
Example using R: Heart Valves Study
Goal: Show that the thrombogenicity rate (TR) is less than two times the objective performance criterion Data: From both the current study and a previous study on a similar product (St. Jude mechanical valve). Model: Let T be the total number of patient-years of followup, and be the TR per year. We assume the number of thrombogenicity events Y P oisson(T ):
e-T (T )y f (y|) = . y!
R and WinBUGS Examples p. 1/4
Example using R: Heart Valves Study
Goal: Show that the thrombogenicity rate (TR) is less than two times the objective performance criterion Data: From both the current study and a previous study on a similar product (St. Jude mechanical valve). Model: Let T be the total number of patient-years of followup, and be the TR per year. We assume the number of thrombogenicity events Y P oisson(T ):
e-T (T )y f (y|) = . y!
Prior: Assume a Gamma(, ) prior for :
-1 e-/ p() = , > 0 . ()
R and WinBUGS Examples p. 1/4
Heart Valves Study
The gamma prior is conjugate with the likelihood, so the posterior emerges in closed form:
p(|y) y+-1 e-(T +1/)
Gamma(y + , (T + 1/)-1 ) .
The study objective is met if
P ( < 2 OP C | y) 0.95 ,
where OP C = 0 = 0.038.
R and WinBUGS Examples p. 2/4
Heart Valves Study
The gamma prior is conjugate with the likelihood, so the posterior emerges in closed form:
p(|y) y+-1 e-(T +1/)
Gamma(y + , (T + 1/)-1 ) .
The study objective is met if
P ( < 2 OP C | y) 0.95 ,
where OP C = 0 = 0.038. Prior selection: Our gamma prior has mean M = and variance V = 2 . This means that if we specify M and V , we can solve for and as
= M 2 /V and = V /M .
R and WinBUGS Examples p. 2/4
Heart Valves Study
A few possibilities for prior parameters:
R and WinBUGS Examples p. 3/4
Heart Valves Study
A few possibilities for prior parameters:
Suppose we set M = 0 = 0.038 and V = 20 (so that 0 is two standard deviations below the mean). Then = 0.25 and = 0.152, a rather vague prior.
R and WinBUGS Examples p. 3/4
Heart Valves Study
A few possibilities for prior parameters:
Suppose we set M = 0 = 0.038 and V = 20 (so that 0 is two standard deviations below the mean). Then = 0.25 and = 0.152, a rather vague prior. Suppose we set M = 98/5891 = .0166, the overall value from the St. Jude studies, and V = M (so 0 is one sd below the mean). Then = 1 and = 0.0166, a moderate (exponential) prior.
R and WinBUGS Examples p. 3/4
Heart Valves Study
A few possibilities for prior parameters:
Suppose we set M = 0 = 0.038 and V = 20 (so that 0 is two standard deviations below the mean). Then = 0.25 and = 0.152, a rather vague prior. Suppose we set M = 98/5891 = .0166, the overall value from the St. Jude studies, and V = M (so 0 is one sd below the mean). Then = 1 and = 0.0166, a moderate (exponential) prior. Suppose we set M = 98/5891 = .0166 again, but set V = M/2. This is a rather informative prior.
R and WinBUGS Examples p. 3/4
Heart Valves Study
A few possibilities for prior parameters:
Suppose we set M = 0 = 0.038 and V = 20 (so that 0 is two standard deviations below the mean). Then = 0.25 and = 0.152, a rather vague prior. Suppose we set M = 98/5891 = .0166, the overall value from the St. Jude studies, and V = M (so 0 is one sd below the mean). Then = 1 and = 0.0166, a moderate (exponential) prior. Suppose we set M = 98/5891 = .0166 again, but set V = M/2. This is a rather informative prior.
We also consider event counts that are lower (1), about the same (3), and much higher (20) than for St. Jude.
R and WinBUGS Examples p. 3/4
Heart Valves Study
A few possibilities for prior parameters:
Suppose we set M = 0 = 0.038 and V = 20 (so that 0 is two standard deviations below the mean). Then = 0.25 and = 0.152, a rather vague prior. Suppose we set M = 98/5891 = .0166, the overall value from the St. Jude studies, and V = M (so 0 is one sd below the mean). Then = 1 and = 0.0166, a moderate (exponential) prior. Suppose we set M = 98/5891 = .0166 again, but set V = M/2. This is a rather informative prior.
We also consider event counts that are lower (1), about the same (3), and much higher (20) than for St. Jude. The study objective is not met with the "bad" data unless the posterior is "rescued" by the informative prior (lower right corner, next page).
R and WinBUGS Examples p. 3/4
Heart Valves Study
Priors and posteriors, Heart Valves ADVANTAGE study, Poisson-gamma model for various prior (M, sd) and data (y) values; T = 200 , 2 OPC = 0.076
80 120
vague prior
60
20 40 60 80
P(theta < 2 OPC|y) = 1
P(theta < 2 OPC|y) = 1
60
posterior prior
posterior prior
80
posterior prior
P(theta < 2 OPC|y) = 0.153
40
20
0
0
0.0
0.05
0.10
0.15
0.20
0.0
0.05
0.10
0.15
0.20
0
0.0
20
40
0.05
0.10
0.15
0.20
M, sd = 0.038 0.076 ; Y = 1
60
M, sd = 0.038 0.076 ; Y = 3
M, sd = 0.038 0.076 ; Y = 20
moderate prior
posterior prior
posterior prior
posterior prior
80
50
60
40
40
P(theta < 2 OPC|y) = 1
P(theta < 2 OPC|y) = 1
30
40
50
P(theta < 2 OPC|y) = 0.421
30
20
20
0
0
10
0.0
0.05
0.10
0.15
0.20
0.0
0.05
0.10
0.15
0.20
0
0.0
10
20
0.05
0.10
0.15
0.20
M, sd = 0.017 0.017 ; Y = 1
M, sd = 0.017 0.017 ; Y = 3
M, sd = 0.017 0.017 ; Y = 20
informative prior
80 60 40 20 0
0.0
40
20
0
0.05
0.10
0.15
0.20
0.0
0.05
0.10
0.15
0.20
0
10
20
P(theta < 2 OPC|y) = 1
P(theta < 2 OPC|y) = 1
30
40
posterior prior
posterior prior
50
posterior prior
60
P(theta < 2 OPC|y) = 0.964
0.0
0.05
0.10
0.15
0.20
M, sd = 0.017 0.008 ; Y = 1
M, sd = 0.017 0.008 ; Y = 3
M, sd = 0.017 0.008 ; Y = 20
R and WinBUGS Examples p. 4/4
Heart Valves Study
Priors and posteriors, Heart Valves ADVANTAGE study, Poisson-gamma model for various prior (M, sd) and data (y) values; T = 200 , 2 OPC = 0.076
80 120
vague prior
60
20 40 60 80
P(theta < 2 OPC|y) = 1
P(theta < 2 OPC|y) = 1
60
posterior prior
posterior prior
80
posterior prior
P(theta < 2 OPC|y) = 0.153
40
20
0
0
0.0
0.05
0.10
0.15
0.20
0.0
0.05
0.10
0.15
0.20
0
0.0
20
40
0.05
0.10
0.15
0.20
M, sd = 0.038 0.076 ; Y = 1
60
M, sd = 0.038 0.076 ; Y = 3
M, sd = 0.038 0.076 ; Y = 20
moderate prior
posterior prior
posterior prior
posterior prior
80
50
60
40
40
P(theta < 2 OPC|y) = 1
P(theta < 2 OPC|y) = 1
30
40
50
P(theta < 2 OPC|y) = 0.421
30
20
20
0
0
10
0.0
0.05
0.10
0.15
0.20
0.0
0.05
0.10
0.15
0.20
0
0.0
10
20
0.05
0.10
0.15
0.20
M, sd = 0.017 0.017 ; Y = 1
M, sd = 0.017 0.017 ; Y = 3
M, sd = 0.017 0.017 ; Y = 20
informative prior
80 60 40 20 0
0.0
40
20
0
0.05
0.10
0.15
0.20
0.0
0.05
0.10
0.15
0.20
0
10
20
P(theta < 2 OPC|y) = 1
P(theta < 2 OPC|y) = 1
30
40
posterior prior
posterior prior
50
posterior prior
60
P(theta < 2 OPC|y) = 0.964
0.0
0.05
0.10
0.15
0.20
M, sd = 0.017 0.008 ; Y = 1
M, sd = 0.017 0.008 ; Y = 3
M, sd = 0.017 0.008 ; Y = 20
S code to create this plot is available in www.biostat.umn.edu/brad/hv.S try it yourself in S-plus or R (http://cran.r-project.org)
R and WinBUGS Examples p. 4/4
Alternate hierarchical models
One might be uncomfortable with our implicit assumption that the TR is the same in both studies. To handle this, extend to a hierarchical model:
Yi P oisson(i Ti ), i = 1, 2,
where i = 1 for St. Jude, and i = 2 for the new study.
R and WinBUGS Examples p. 5/4
Alternate hierarchical models
One might be uncomfortable with our implicit assumption that the TR is the same in both studies. To handle this, extend to a hierarchical model:
Yi P oisson(i Ti ), i = 1, 2,
where i = 1 for St. Jude, and i = 2 for the new study. Borrow strength between studies by assuming
i Gamma(, ),
iid
i.e., the two TR's are exchangeable, but not identical.
R and WinBUGS Examples p. 5/4
Alternate hierarchical models
One might be uncomfortable with our implicit assumption that the TR is the same in both studies. To handle this, extend to a hierarchical model:
Yi P oisson(i Ti ), i = 1, 2,
where i = 1 for St. Jude, and i = 2 for the new study. Borrow strength between studies by assuming
i Gamma(, ),
iid
i.e., the two TR's are exchangeable, but not identical. We now place a third stage prior on and , say
Exp(a) and IG(c, d).
R and WinBUGS Examples p. 5/4
Alternate hierarchical models
One might be uncomfortable with our implicit assumption that the TR is the same in both studies. To handle this, extend to a hierarchical model:
Yi P oisson(i Ti ), i = 1, 2,
where i = 1 for St. Jude, and i = 2 for the new study. Borrow strength between studies by assuming
i Gamma(, ),
iid
i.e., the two TR's are exchangeable, but not identical. We now place a third stage prior on and , say
Exp(a) and IG(c, d).
Fit in WinBUGS using the pump example as a guide!
R and WinBUGS Examples p. 5/4
Pump Example
Example 2.5 revisited again!
i G(, ), Exp(), IG(c, d), i = 1, . . . , k , where , c, d, and the ti are known, and Exp denotes the exponential distribution. Yi |i P oisson(i ti ),
ind ind
R and WinBUGS Examples p. 6/4
Pump Example
Example 2.5 revisited again!
i G(, ), Exp(), IG(c, d), i = 1, . . . , k , where , c, d, and the ti are known, and Exp denotes the exponential distribution. Yi |i P oisson(i ti ),
ind ind
We apply this model to a dataset giving the numbers of pump failures, Yi , observed in ti thousands of hours for k = 10 different systems of a certain nuclear power plant.
R and WinBUGS Examples p. 6/4
Pump Example
Example 2.5 revisited again!
i G(, ), Exp(), IG(c, d), i = 1, . . . , k , where , c, d, and the ti are known, and Exp denotes the exponential distribution. Yi |i P oisson(i ti ),
ind ind
We apply this model to a dataset giving the numbers of pump failures, Yi , observed in ti thousands of hours for k = 10 different systems of a certain nuclear power plant. The observations are listed in increasing order of raw failure rate ri = Yi /ti , the classical point estimate of the true failure rate i for the ith system.
R and WinBUGS Examples p. 6/4
Pump Data
i 1 2 3 4 5 6 7 8 9 10 Yi 5 1 5 14 3 19 1 1 4 22 ti 94.320 15.720 62.880 125.760 5.240 31.440 1.048 1.048 2.096 10.480 ri .053 .064 .080 .111 .573 .604 .954 .954 1.910 2.099
Hyperparameters: We choose the values = 1, c = 0.1, and d = 1.0, resulting in reasonably vague hyperpriors for and .
R and WinBUGS Examples p. 7/4
Pump Example
Recall that the full conditional distributions for the i and are available in closed form (gamma and inverse gamma, respectively), but that no conjugate prior for exists.
R and WinBUGS Examples p. 8/4
Pump Example
Recall that the full conditional distributions for the i and are available in closed form (gamma and inverse gamma, respectively), but that no conjugate prior for exists. However, the full conditional for ,
k
p(|, {i }, y)
i=1 k i=1
g(i |, ) h()
-1 i e-/ ()
can be shown to be log-concave in . Thus WinBUGS uses adaptive rejection sampling for this parameter.
R and WinBUGS Examples p. 8/4
WinBUGS code to fit this model
model { for (i in 1:k) { theta[i] ~ dgamma(alpha,beta) lambda[i] <- theta[i]*t[i] Y[i] ~ dpois(lambda[i]) } alpha ~ dexp(1.0) beta ~ dgamma(0.1, 1.0) } DATA: list(k = 10, Y = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22), t = c(94.320, 15.72, 62.88, 125.76, 5.24, 31.44, 1.048, 1.048, 2.096, 10.48)) INITS: list(theta=c(1,1,1,1,1,1,1,1,1,1), alpha=1, beta=1)
R and WinBUGS Examples p. 9/4
Pump Example Results
Results from running 1000 burn-in samples, followed by a "production" run of 10,000 samples (single chain):
node alpha beta theta[1] theta[5] theta[6] theta[10] mean 0.7001 0.929 0.0598 0.6056 0.6105 1.993 sd 0.2699 0.5325 0.02542 0.315 0.1393 0.4251 MC error 0.004706 0.00978 2.68E-4 0.003087 0.0014 0.004915 2.5% 0.2851 0.1938 0.02128 0.1529 0.3668 1.264 median 0.6634 0.8315 0.05627 0.5529 0.5996 1.958 97.5% 1.338 2.205 0.1195 1.359 0.9096 2.916
R and WinBUGS Examples p. 10/4
Pump Example Results
Results from running 1000 burn-in samples, followed by a "production" run of 10,000 samples (single chain):
node alpha beta theta[1] theta[5] theta[6] theta[10] mean 0.7001 0.929 0.0598 0.6056 0.6105 1.993 sd 0.2699 0.5325 0.02542 0.315 0.1393 0.4251 MC error 0.004706 0.00978 2.68E-4 0.003087 0.0014 0.004915 2.5% 0.2851 0.1938 0.02128 0.1529 0.3668 1.264 median 0.6634 0.8315 0.05627 0.5529 0.5996 1.958 97.5% 1.338 2.205 0.1195 1.359 0.9096 2.916
Note that while 5 and 6 have very similar posterior means, the latter posterior is much narrower (smaller sd).
R and WinBUGS Examples p. 10/4
Pump Example Results
Results from running 1000 burn-in samples, followed by a "production" run of 10,000 samples (single chain):
node alpha beta theta[1] theta[5] theta[6] theta[10] mean 0.7001 0.929 0.0598 0.6056 0.6105 1.993 sd 0.2699 0.5325 0.02542 0.315 0.1393 0.4251 MC error 0.004706 0.00978 2.68E-4 0.003087 0.0014 0.004915 2.5% 0.2851 0.1938 0.02128 0.1529 0.3668 1.264 median 0.6634 0.8315 0.05627 0.5529 0.5996 1.958 97.5% 1.338 2.205 0.1195 1.359 0.9096 2.916
Note that while 5 and 6 have very similar posterior means, the latter posterior is much narrower (smaller sd). This is because, while the crude failure rates for the two pumps are similar, the latter is based on a far greater number of hours of observation (t6 = 31.44, while t5 = 5.24). Hence we "know" more about pump 6!
R and WinBUGS Examples p. 10/4
Cross-Study (Meta-analysis) Example
^ Data: estimated log relative hazards Yij = ij obtained by fitting separate Cox proportional hazards regressions to the data from each of J = 18 clinical units participating in I = 6 different AIDS studies.
R and WinBUGS Examples p. 11/4
Cross-Study (Meta-analysis) Example
^ Data: estimated log relative hazards Yij = ij obtained by fitting separate Cox proportional hazards regressions to the data from each of J = 18 clinical units participating in I = 6 different AIDS studies.
To these data we wish to fit the cross-study model,
Yij = ai + bj + sij + ij , i = 1, . . . , I, j = 1, . . . , J,
where ai = study main effect bj = unit main effect sij = study-unit interaction term, and
2 ij N (0, ij ) iid
and the estimated standard errors from the Cox regressions are used as (known) values of the ij .
R and WinBUGS Examples p. 11/4
Cross-Study (Meta-analysis) Data
Estimated Unit-Specific Log Relative Hazards Toxo Unit A B C D E F G H I J K . . . R 0.814 -0.203 -0.133 NA -0.715 0.739 0.118 NA NA 0.271 NA . . . 1.217 NA NA NA NA -0.242 0.009 0.807 -0.511 1.939 1.079 NA . . . 0.165 ddI/ddC NuCombo ZDV+ddI -0.406 NA 0.218 NA -0.544 NA -0.047 0.233 0.218 -0.277 0.792 . . . 0.385 NuCombo ZDV+ddC 0.298 NA -2.206 NA -0.731 NA 0.913 0.131 -0.066 -0.232 1.264 . . . 0.172 0.094 NA 0.435 NA 0.600 NA -0.091 NA NA 0.752 -0.357 . . . -0.022 NA NA 0.145 NA 0.041 0.222 0.099 0.017 0.355 0.203 0.807 . . . 0.203 Fungal CMV
R and WinBUGS Examples p. 12/4
Cross-Study (Meta-analysis) Example
Note that some values are missing ("NA") since
R and WinBUGS Examples p. 13/4
Cross-Study (Meta-analysis) Example
Note that some values are missing ("NA") since not all 18 units participated in all 6 studies
R and WinBUGS Examples p. 13/4
Cross-Study (Meta-analysis) Example
Note that some values are missing ("NA") since not all 18 units participated in all 6 studies the Cox estimation procedure did not converge for some units that had few deaths
R and WinBUGS Examples p. 13/4
Cross-Study (Meta-analysis) Example
Note that some values are missing ("NA") since not all 18 units participated in all 6 studies the Cox estimation procedure did not converge for some units that had few deaths Goal: To identify which clinics are opinion leaders (strongly agree with overall result across studies) and which are dissenters (strongly disagree).
R and WinBUGS Examples p. 13/4
Cross-Study (Meta-analysis) Example
Note that some values are missing ("NA") since not all 18 units participated in all 6 studies the Cox estimation procedure did not converge for some units that had few deaths Goal: To identify which clinics are opinion leaders (strongly agree with overall result across studies) and which are dissenters (strongly disagree). Here, overall results all favor the treatment (i.e. mostly negative Y s) except in Trial 1 (Toxo). Thus we multiply all the Yij 's by 1 for i = 1, so that larger Yij correspond in all cases to stronger agreement with the overall.
R and WinBUGS Examples p. 13/4
Cross-Study (Meta-analysis) Example
Note that some values are missing ("NA") since not all 18 units participated in all 6 studies the Cox estimation procedure did not converge for some units that had few deaths Goal: To identify which clinics are opinion leaders (strongly agree with overall result across studies) and which are dissenters (strongly disagree). Here, overall results all favor the treatment (i.e. mostly negative Y s) except in Trial 1 (Toxo). Thus we multiply all the Yij 's by 1 for i = 1, so that larger Yij correspond in all cases to stronger agreement with the overall. Next slide shows a plot of the Yij values and associated approximate 95% CIs...
R and WinBUGS Examples p. 13/4
Cross-Study (Meta-analysis) Data
1: Toxo 2: ddI/ddC 3: NuCombo-ddI
Placebo better ddI better ZDV better
4
4
Estimated Log Relative Hazard
Estimated Log Relative Hazard
Estimated Log Relative Hazard
ddC better
2
2
0
0
-2
-2
-4
-4
Trt better
-4
-2
0
2
4
Combo better
A B C D E F G H
I
J K L M N O P Q R
A B C D E F G H
I
J K
L M N O P Q R
A B C D E F G H
I
J
K L M N O P Q R
Clinic
Clinic
Clinic
4: NuCombo-ddC
5: Fungal
6: CMV
ZDV better
Placebo better
Placebo better
4
4
Estimated Log Relative Hazard
Estimated Log Relative Hazard
Estimated Log Relative Hazard
Trt better
2
2
0
0
-2
-2
-4
-4
Combo better
-4
-2
0
2
4
Trt better
A B C D E F G H
I
J K L M N O P Q R
A B C D E F G H
I
J K
L M N O P Q R
A B C D E F G H
I
J
K L M N O P Q R
Clinic
Clinic
Clinic
R and WinBUGS Examples p. 14/4
Cross-Study (Meta-analysis) Example
Second stage of our model:
ai N (0, A-1 ), bj N (0, B -1 ), and sij N (0, C -1 )
iid iid iid
R and WinBUGS Examples p. 15/4
Cross-Study (Meta-analysis) Example
Second stage of our model:
ai N (0, A-1 ), bj N (0, B -1 ), and sij N (0, C -1 )
iid iid iid
Third stage of our model:
A = 0.0001 (the ai 's are fixed effects) B G(0.01, 0.01) (the bj 's are random effects)
C G(0.01, 0.01) (the sij 's are random effects)
That is, we
R and WinBUGS Examples p. 15/4
Cross-Study (Meta-analysis) Example
Second stage of our model:
ai N (0, A-1 ), bj N (0, B -1 ), and sij N (0, C -1 )
iid iid iid
Third stage of our model:
A = 0.0001 (the ai 's are fixed effects) B G(0.01, 0.01) (the bj 's are random effects)
C G(0.01, 0.01) (the sij 's are random effects)
That is, we preclude borrowing of strength across studies, but
R and WinBUGS Examples p. 15/4
Cross-Study (Meta-analysis) Example
Second stage of our model:
ai N (0, A-1 ), bj N (0, B -1 ), and sij N (0, C -1 )
iid iid iid
Third stage of our model:
A = 0.0001 (the ai 's are fixed effects) B G(0.01, 0.01) (the bj 's are random effects)
C G(0.01, 0.01) (the sij 's are random effects)
That is, we preclude borrowing of strength across studies, but encourage borrowing of strength across units
R and WinBUGS Examples p. 15/4
Cross-Study (Meta-analysis) Example
Second stage of our model:
ai N (0, A-1 ), bj N (0, B -1 ), and sij N (0, C -1 )
iid iid iid
Third stage of our model:
A = 0.0001 (the ai 's are fixed effects) B G(0.01, 0.01) (the bj 's are random effects)
C G(0.01, 0.01) (the sij 's are random effects)
That is, we preclude borrowing of strength across studies, but encourage borrowing of strength across units With I + J + IJ parameters but fewer than IJ data points, some effects must be treated as random!
R and WinBUGS Examples p. 15/4
WinBUGS code to fit this model
model { for(i in 1:I) { for(j in 1:J) { Y[i,j] ~ dnorm(theta[i,j],P[i,j]) theta[i,j] <- a[i]+b[j]+s[i,j] s[i,j] ~ dnorm(0.0,C) sigma[i,j] <- 1/sqrt(P[i,j]) } a[i] ~ dnorm(0.0,1.0E-4) } for(j in 1:J) { b[j] ~ dnorm(0.0,B) } B ~ dgamma(0.01,0.01) # unit effect precision C ~ dgamma(0.01,0.01) # interaction precision }
R and WinBUGS Examples p. 16/4
Plot of ij posterior means
Placebo better ddC better Combo better Combo better Trt better Trt better
0.4
0.6
Log Relative Hazard
P P R F IK Q A N G M BD J C H L O E I P K P Q KM R G P I PR Q A D G KMO N B F J C L H E R A DF M Q B G C J LN H E O I
G KM Q R AD J B F C LN O H E
0.2
0.0
I G DF B C N H A JL O E
I A F M BD H C J O L E
KNQ R
-0.2
Toxo
ddI/ddC
NuCombo(ddI)
NuCombo(ddC)
Fungal
CMV
1
2
3 Study
4
5
6
Unit P is an opinion leader; Unit E is a dissenter
R and WinBUGS Examples p. 17/4
Plot of ij posterior means
Placebo better ddC better Combo better Combo better Trt better Trt better
0.4
0.6
Log Relative Hazard
P P R F IK Q A N G M BD J C H L O E I P K P Q KM R G P I PR Q A D G KMO N B F J C L H E R A DF M Q B G C J LN H E O I
G KM Q R AD J B F C LN O H E
0.2
0.0
I G DF B C N H A JL O E
I A F M BD H C J O L E
KNQ R
-0.2
Toxo
ddI/ddC
NuCombo(ddI)
NuCombo(ddC)
Fungal
CMV
1
2
3 Study
4
5
6
Unit P is an opinion leader; Unit E is a dissenter Substantial shrinkage towards 0 has occurred: mostly positive values; no estimated ij greater than 0.6
R and WinBUGS Examples p. 17/4
PK Example
Wakefield et al. (1994) consider a dataset for which
Yij = plasma concentration of the drug Cadralazine xij = time elapsed since dose given
where i = 1, . . . , 10 indexes the patient, while j = 1, . . . , ni indexes the observations, 5 ni 8.
R and WinBUGS Examples p. 18/4
PK Example
Wakefield et al. (1994) consider a dataset for which
Yij = plasma concentration of the drug Cadralazine xij = time elapsed since dose given
Attempt to fit the one-compartment nonlinear pharmacokinetic (PK) model,
-1 ij (xij ) = 30i exp(-i xij /i ) .
where i = 1, . . . , 10 indexes the patient, while j = 1, . . . , ni indexes the observations, 5 ni 8.
where ij (xij ) is the mean plasma concentration at time xij .
R and WinBUGS Examples p. 18/4
PK Example
This model is best fit on the log scale, i.e.
Zij log Yij = log ij (xij ) + ij ,
where ij N i (0, ).
ind
R and WinBUGS Examples p. 19/4
PK Example
This model is best fit on the log scale, i.e.
Zij log Yij = log ij (xij ) + ij ,
where ij N (0, i ).
ind
The mean structure for the Zij 's thus emerges as
-1 log ij (xij ) = log 30i exp(-i xij /i )
= log 30 - ai - exp(bi - ai )xij ,
= log 30 - log i - i xij /i
where ai = log i and bi = log i .
R and WinBUGS Examples p. 19/4
PK Data
no. of hours following drug administration, x patient 1 2 3 4 5 6 7 8 9 10 2 1.09 2.03 1.44 1.55 1.35 1.08 1.32 1.63 1.26 1.30 4 0.75 1.28 1.30 0.96 0.78 0.59 0.74 1.01 0.73 0.70 6 0.53 1.20 0.95 0.80 0.50 0.37 0.46 0.73 0.40 0.40 8 0.34 1.02 0.68 0.62 0.33 0.23 0.28 0.55 0.30 0.25 10 0.23 0.83 0.52 0.46 0.18 0.17 0.27 0.41 0.21 0.14 24 0.02 0.28 0.06 0.08 0.02 0.03 0.01 28 0.02 0.06 32 0.02
R and WinBUGS Examples p. 20/4
PK Data, original scale
2.0
2
1.5
8 4 3 5 7 1 9 1 6
3 2 2
Y
1.0
8 4 5 1 7 9 1 6
2 3 4 8 1 5 7 9 1 6 2 3 4 8
3 4 8 7 1 9 5 6 1 2 4 3 7 5 1 8
0.5
1 5 9 7 1 6
0.0
8 7
8
5
10
15 x
20
25
30
R and WinBUGS Examples p. 21/4
PK Data, log scale
1 0
2 8 4 3 5 7 1 9 1 6 3 2 8 4 5 1 7 9 1 6 2 3 4 8 1 5 7 9 1 6 2 2 3 4 8 1 5 9 7 1 6 3 4 8 7 1 9 5 6 1 2
log Y
-2
-1
4
-3
3
8
7
-4
5 1
7
8
8
-5
5
10
15 x
20
25
30
R and WinBUGS Examples p. 22/4
PK Example
For the subject-specific random effects i (ai , bi ) , assume
i N2 (, ) , where = (a , b ) .
iid
R and WinBUGS Examples p. 23/4
PK Example
For the subject-specific random effects i (ai , bi ) , assume
i N2 (, ) , where = (a , b ) .
iid
Usual conjugate prior specification:
N2 (, C) i G(0 /2 , 0 0 /2)
iid
Wishart((R)-1 , )
R and WinBUGS Examples p. 23/4
PK Example
For the subject-specific random effects i (ai , bi ) , assume
i N2 (, ) , where = (a , b ) .
iid
Usual conjugate prior specification:
N2 (, C) i G(0 /2 , 0 0 /2)
iid
Wishart((R)-1 , )
Note that the i full conditional distributions are: not simple conjugate forms not guaranteed to be log-concave Thus, the Metropolis capability of WinBUGS is required.
R and WinBUGS Examples p. 23/4
PK Results (WinBUGS vs. Fortran)
BUGS parameter a1 a2 a7 a8 b1 b2 b7 b8 1 2 7 8 Y2,8 Y7,8 mean 2.956 2.692 2.970 2.828 1.259 0.234 1.157 0.936 362.4 84.04 18.87 2.119 0.1338 0.00891 sd 0.0479 0.0772 0.1106 0.1417 0.0335 0.0648 0.0879 0.1458 260.4 57.60 12.07 1.139 0.0339 0.00443 lag 1 acf 0.969 0.769 0.925 0.828 0.972 0.661 0.899 0.759 0.313 0.225 0.260 0.085 0.288 0.178 Sargent et al. (2000) mean 2.969 2.708 2.985 2.838 1.268 0.239 1.163 0.941 380.8 81.40 15.82 1.499 0.1347 0.00884 sd 0.0460 0.0910 0.1360 0.1863 0.0322 0.0798 0.1055 0.1838 268.8 58.41 11.12 0.931 0.0264 0.00255 lag 1 acf 0.947 0.808 0.938 0.934 0.951 0.832 0.925 0.932 0.220 0.255 0.237 0.143
R and WinBUGS Examples p. 24/4
Lip Cancer Example
Consider the spatial disease mapping model:
Yi | i Yi Ei i
ind
= = =
P o (Ei ei ) , where observed disease count, expected count (known), and x + i + i i
R and WinBUGS Examples p. 25/4
Lip Cancer Example
Consider the spatial disease mapping model:
Yi | i Yi Ei i
ind
= = =
P o (Ei ei ) , where observed disease count, expected count (known), and x + i + i i
The xi are explanatory spatial covariates; typically is assigned a flat prior.
R and WinBUGS Examples p. 25/4
Lip Cancer Example
Consider the spatial disease mapping model:
Yi | i Yi Ei i
ind
= = =
P o (Ei ei ) , where observed disease count, expected count (known), and x + i + i i
The xi are explanatory spatial covariates; typically is assigned a flat prior. Note the mean structure also contains two sets of random effects! The first, i , capture heterogeneity among the regions via
i N (0 , 1/h ) .
iid
R and WinBUGS Examples p. 25/4
Lip Cancer Example
The second set, i , capture regional clustering via a conditionally autoregressive (CAR) prior,
i | j=i N (i , 1/(c mi )) ,
where mi is the number of "neighbors" of region i, and i = m-1 ji j . i
R and WinBUGS Examples p. 26/4
Lip Cancer Example
The second set, i , capture regional clustering via a conditionally autoregressive (CAR) prior,
i | j=i N (i , 1/(c mi )) ,
where mi is the number of "neighbors" of region i, and i = m-1 ji j . i The CAR prior is translation invariant, so typically we insist I i = 0 (imposed numerically after each i=1 MCMC iteration). Still, Yi cannot inform about i or i , but only about their sum i = i + i .
R and WinBUGS Examples p. 26/4
Lip Cancer Example
The second set, i , capture regional clustering via a conditionally autoregressive (CAR) prior,
i | j=i N (i , 1/(c mi )) ,
where mi is the number of "neighbors" of region i, and i = m-1 ji j . i The CAR prior is translation invariant, so typically we insist I i = 0 (imposed numerically after each i=1 MCMC iteration). Still, Yi cannot inform about i or i , but only about their sum i = i + i . Making the reparametrization from (, ) to (, ), we have the joint posterior
p(, |y) L(; y)p()p(-).
R and WinBUGS Examples p. 26/4
Lip Cancer Example
This means that
p(i | j=i , , y) p(i ) p(i -i | {j -j }j=i ) .
R and WinBUGS Examples p. 27/4
Lip Cancer Example
This means that
p(i | j=i , , y) p(i ) p(i -i | {j -j }j=i ) .
Since this distribution is free of the data y, the i are Bayesianly unidentified (and so are the i ).
R and WinBUGS Examples p. 27/4
Lip Cancer Example
This means that
p(i | j=i , , y) p(i ) p(i -i | {j -j }j=i ) .
Since this distribution is free of the data y, the i are Bayesianly unidentified (and so are the i ). BUT this does not preclude Bayesian learning about i ; this would instead require
p(i |y) = p(i ) .
[Stronger condition: data have no impact on the marginal (not conditional) posterior.]
R and WinBUGS Examples p. 27/4
Lip Cancer Example
Dilemma: Though unidentified, the i and i are interesting in their own right, as is
sd() = , sd() + sd()
where sd() is the empirical marginal standard deviation function. Can we specify vague but proper prior values h and c that
R and WinBUGS Examples p. 28/4
Lip Cancer Example
Dilemma: Though unidentified, the i and i are interesting in their own right, as is
sd() = , sd() + sd()
where sd() is the empirical marginal standard deviation function. Can we specify vague but proper prior values h and c that lead to acceptable convergence behavior, and
R and WinBUGS Examples p. 28/4
Lip Cancer Example
Dilemma: Though unidentified, the i and i are interesting in their own right, as is
sd() = , sd() + sd()
where sd() is the empirical marginal standard deviation function. Can we specify vague but proper prior values h and c that lead to acceptable convergence behavior, and still allow Bayesian learning?
R and WinBUGS Examples p. 28/4
Lip Cancer Example
Dilemma: Though unidentified, the i and i are interesting in their own right, as is
sd() = , sd() + sd()
where sd() is the empirical marginal standard deviation function. Can we specify vague but proper prior values h and c that lead to acceptable convergence behavior, and still allow Bayesian learning? Tricky to specify a "fair" prior balance between heterogeneity and clustering (e.g., one for which 1/2) since i prior is specified marginally while the i prior is specified conditionally.
R and WinBUGS Examples p. 28/4
Dataset: Scottish lip cancer data
a) b)
> 200 150 - 200 105 - 150 95 - 105 50 - 95 < 50
> 16 12 - 16 8 - 12 7 - 8 1 - 7 0 - 1
N
60 0 60 120 Miles
a) SM Ri = 100Yi /Ei , standardized mortality ratio for lip cancer in I = 56 districts, 19751980
R and WinBUGS Examples p. 29/4
Dataset: Scottish lip cancer data
a) b)
> 200 150 - 200 105 - 150 95 - 105 50 - 95 < 50
> 16 12 - 16 8 - 12 7 - 8 1 - 7 0 - 1
N
60 0 60 120 Miles
a) SM Ri = 100Yi /Ei , standardized mortality ratio for lip cancer in I = 56 districts, 19751980 b) one covariate, xi = percentage of the population engaged in agriculture, fishing or forestry (AFF)
R and WinBUGS Examples p. 29/4
WinBUGS code to fit this model
model { for (i in 1 : regions) { O[i] ~ dpois(mu[i]) log(mu[i]) <- log(E[i]) + beta*aff[i]/10 + phi[i] + theta[i] theta[i] ~ dnorm(0.0,tau.h) eta[i] <- theta[i] + phi[i] } phi[1:regions] ~ {\red car.normal}(adj[], weights[], num[], tau.c) beta ~ dnorm(0.0, 1.0E-5) # vague prior on covariate effect # # ``fair'' prior from Best et al. (1999, Bayesian Statistics 6)
tau.h ~ dgamma(1.0E-3,1.0E-3) tau.c ~ dgamma(1.0E-1,1.0E-1)
sd.h <- sd(theta[]) # marginal SD of heterogeneity effects sd.c <- sd(phi[]) # marginal SD of clustering (spatial) effects psi <- sd.c / (sd.h + sd.c) } (See WinBUGS Map manual for DATA and INITS for this example)
R and WinBUGS Examples p. 30/4
Lip Cancer Results
posterior for priors for c , h G(1.0, 1.0), G(3.2761, 1.81) G(.1, .1), G(.32761, .181) G(.1, .1), G(.001, .001) priors for c , h G(1.0, 1.0), G(3.2761, 1.81) G(.1, .1), G(.32761, .181) G(.1, .1), G(.001, .001) mean .57 .65 .82 mean .92 .89 .90 sd .058 .073 .10 sd .40 .36 .34 l1acf .80 .89 .98 l1acf .33 .28 .31 posterior for mean .43 .41 .38 mean .96 .79 .70 sd .17 .14 .13 sd .52 .41 .35 l1acf .94 .92 .91 l1acf .12 .17 .21
posterior for 1
posterior for 56
AFF covariate is significantly = 0 under all 3 priors
R and WinBUGS Examples p. 31/4
Lip Cancer Results
posterior for priors for c , h G(1.0, 1.0), G(3.2761, 1.81) G(.1, .1), G(.32761, .181) G(.1, .1), G(.001, .001) priors for c , h G(1.0, 1.0), G(3.2761, 1.81) G(.1, .1), G(.32761, .181) G(.1, .1), G(.001, .001) mean .57 .65 .82 mean .92 .89 .90 sd .058 .073 .10 sd .40 .36 .34 l1acf .80 .89 .98 l1acf .33 .28 .31 posterior for mean .43 .41 .38 mean .96 .79 .70 sd .17 .14 .13 sd .52 .41 .35 l1acf .94 .92 .91 l1acf .12 .17 .21
posterior for 1
posterior for 56
AFF covariate is significantly = 0 under all 3 priors convergence is slow for and , but rapid for i (and i )
R and WinBUGS Examples p. 31/4
Lip Cancer Results
posterior for priors for c , h G(1.0, 1.0), G(3.2761, 1.81) G(.1, .1), G(.32761, .181) G(.1, .1), G(.001, .001) priors for c , h G(1.0, 1.0), G(3.2761, 1.81) G(.1, .1), G(.32761, .181) G(.1, .1), G(.001, .001) mean .57 .65 .82 mean .92 .89 .90 sd .058 .073 .10 sd .40 .36 .34 l1acf .80 .89 .98 l1acf .33 .28 .31 posterior for mean .43 .41 .38 mean .96 .79 .70 sd .17 .14 .13 sd .52 .41 .35 l1acf .94 .92 .91 l1acf .12 .17 .21
posterior for 1
posterior for 56
AFF covariate is significantly = 0 under all 3 priors convergence is slow for and , but rapid for i (and i ) Excess variability in the data is mostly due to clustering (E(|y) > .50), but the posterior distribution for does not seem robust to changes in the prior.
R and WinBUGS Examples p. 31/4
MAC Survival Example
Data arises from a clinical trial comparing two treatments for Mycobacterium avium complex (MAC), a disease common in late stage HIV-infected persons. Eleven clinical centers ("units") have enrolled a total of 69 patients in the trial, of which 18 have died.
R and WinBUGS Examples p. 32/4
MAC Survival Example
Data arises from a clinical trial comparing two treatments for Mycobacterium avium complex (MAC), a disease common in late stage HIV-infected persons. Eleven clinical centers ("units") have enrolled a total of 69 patients in the trial, of which 18 have died. For j = 1, . . . , ni and i = 1, . . . , k , let
tij = time to death or censoring xij =
treatment indicator for subject j in stratum i
R and WinBUGS Examples p. 32/4
MAC Survival Example
Data arises from a clinical trial comparing two treatments for Mycobacterium avium complex (MAC), a disease common in late stage HIV-infected persons. Eleven clinical centers ("units") have enrolled a total of 69 patients in the trial, of which 18 have died. For j = 1, . . . , ni and i = 1, . . . , k , let
tij = time to death or censoring xij =
treatment indicator for subject j in stratum i
Next page gives survival times (in half-days) from the MAC treatment trial, where "+" indicates a censored observation...
R and WinBUGS Examples p. 32/4
MAC Survival Example
unit A A A A B B C D D D D drug 1 2 1 2 2 1 2 2 2 2 2 time 74+ 248 272+ 344 F 4+ 156+ 100+ 20+ 64 88 148+ F F F F F F F F F 1 2 1 2 2 1 1 2 1 1 6 16+ 76 80 202 258+ 268+ 368+ 380+ 424+ I I I I I I I I I K 2 2 2 1 1 2 1 2 1 2 8 16+ 40 120+ 168+ 174+ 268+ 276 286+ 106+ unit E E E drug 1 2 2 time 214 228+ 262 unit H H H H drug 1 1 1 2 time 74+ 88+ 148+ 162
R and WinBUGS Examples p. 33/4
MAC Survival Example
With proportional hazards and a Weibull baseline hazard, stratum i's hazard is
h(tij ; xij ) = h0 (tij )i exp(0 + 1 xij ) = i ti -1 exp(0 + 1 xij + Wi ) , ij
where i > 0, = (0 , 1 ) 2 , and Wi = log i is a clinic-specific frailty term.
R and WinBUGS Examples p. 34/4
MAC Survival Example
With proportional hazards and a Weibull baseline hazard, stratum i's hazard is
h(tij ; xij ) = h0 (tij )i exp(0 + 1 xij ) = i ti -1 exp(0 + 1 xij + Wi ) , ij
where i > 0, = (0 , 1 ) 2 , and Wi = log i is a clinic-specific frailty term. The Wi capture overall differences among the clinics, while the i allow differing baseline hazards which either increase (i > 1) or decrease (i < 1) over time. We assume i.i.d. specifications for these random effects,
Wi N (0, 1/ ) and i G(, ) .
iid iid
R and WinBUGS Examples p. 34/4
MAC Survival Example
As in the mice example (WinBUGS Examples Vol 1),
ij = exp(0 + 1 xij + Wi ) ,
so that
tij W eibull(i , ij ) .
R and WinBUGS Examples p. 35/4
MAC Survival Example
As in the mice example (WinBUGS Examples Vol 1),
ij = exp(0 + 1 xij + Wi ) ,
so that
tij W eibull(i , ij ) .
We recode the drug covariate from (1,2) to (1,1) (i.e., set xij = 2drugij - 3) to ease collinearity betwee...