BUGS CODE 4.1
# Measures of accuracy
# Binary Scores
Model;
{
# Dirichlet distribution for cell probabilities
g00~dgamma(a00,2)
g01~dgamma(a01,2)
g10~dgamma(a10,2)
g11~dgamma(a11,2)
h<-g00+g01+g10+g11
# the theta have a Dirichlet distribution
theta00<-g00/h
theta01<-g01/h
theta10<-g10/h
theta11<-g11/h
# the basic test accuracies are below
tpf<-theta11/(theta11+theta01)
se<-tpf
sp<-1-fpf
fpf<-theta10/(theta10+theta00)
tnf<-theta00/(theta00+theta10)
fnf<-theta01/(theta01+theta11)
ppv<-theta11/(theta10+theta11)
npv<-theta00/(theta00+theta01)
pdlr<-tpf/fpf
ndlr<-fnf/tnf
}
# Exercise Stress Test Pepe
# Uniform Prior (add one to each cell of the table frequencies !)
list(a00=328,a01=209,a10=116,a11=819)
# chest pain history
# Uniform prior
list(a00=198,a01=55,a10=246,a11=970)
# initial values
list(g00=1,g01=1,g10=1,g11=1)
list(g00=2,g01=2,g10=2,g11=2)
list(g00=3,g01=3,g10=3,g11=3)
BUGS CODE 4.2
# Area under the curve
# Ordinal values
# Five values
Model;
{
# generate Dirichlet distribution
g11~dgamma(a11,2)
g12~dgamma(a12,2)
g13~dgamma(a13,2)
g14~dgamma(a14,2)
g15~dgamma(a15,2)
g01~dgamma(a01,2)
g02~dgamma(a02,2)
g03~dgamma(a03,2)
g04~dgamma(a04,2)
g05~dgamma(a05,2)
g1<-g11+g12+g13+g14+g15
g0<-g01+g02+g03+g04+g05
# posterior distribution of probabilities for response of diseased patients
theta1<-g11/g1
theta2<-g12/g1
theta3<-g13/g1
theta4<-g14/g1
theta5<-g15/g1
# posterior distribution for probabilities of response of non-diseased patients
ph1<-g01/g0
ph2<-g02/g0
ph3<-g03/g0
ph4<-g04/g0
ph5<-g05/g0
# auc is area under ROC curve
#A1 is the P[Y>X]
#A2 is the P[Y=X]
auc<- A1+A2/2
A1<-theta2*ph1+theta3*(ph1+ph2)+theta4*(ph1+ph2+ph3)+
theta5*(ph1+ph2+ph3+ph4)
A2<- theta1*ph1+theta2*ph2+theta3*ph3+theta4*ph4
+theta5*ph5
}
# Mammography Example Zhou et al.
# Uniform Prior
# see Table 4.7
list(a11=2,a12=1,a13=7,a14=12,a15=13,a01=10,a02=3,a03=12,
a04=9,a05=1)
# Gallium citrate Example Zhou et al.
# Uniform Prior
list(a11=13,a12=7,a13=4,a14=2,a15=19,a01=12,a02=3,a03=4,
a04=2,a05=4)
# initial values
list(g11=1,g12=1,g13=1,g14=1,g15=1,g01=1,g02=1,g03=1,
g04=1,g05=1)
BUGS CODE 4.3a
model
{
# area under ROC curve
auc <- phi((mu2-mu1)/sqrt(1/tau1 + 1/tau2))
# This is for the non diseased patients
mu1~ dnorm(113.14, prec1)
# this is for the diseased
mu2~dnorm (122.57, prec2)
prec1<-n1*tau1
prec2<-n2*tau2
# a1=n1/2
# b1= (n1-1) (sample variance 1)/2
tau1~dgamma( a1,b1)
# a2=n2/2
# b2= (n2-1)* (sample variance 2)/2
tau2 ~dgamma(a2,b2)
# binormal parameters
a<-(mu2-mu1)*sqrt(tau2)
b<-sqrt(tau1/tau2)
# c is the proportion of auc values > .80
c<-step(auc-.80)
}
# this is the data for diabetes study
list( n1=19, n2=60, a1=9.5, b1= 745.128, a2= 30, b2= 1427.74)
# these are the initial values
list( mu1= 0, mu2= 0, tau1 =1, tau2 =1)
BUGS CODE 4.3b
model;
# Calculates posterior distribution of model parameters and the area under curve. y=test
# Based on O’Mally et al. regression method
{
# likelihood function
for(i in 1:N) {
y[i]~dnorm(mu[i],precy[d[i]+1]);
# yt[i] <- log(y[i]); # logarithmic transformation
mu[i] <- beta[1] + beta[2]*d[i] ;
}
# prior distributions - non-informative prior; similarly for informative priors
for(i in 1:P) {
beta[i] ~ dnorm(0, 0.000001);
}
for(i in 1:K) {
precy[i]~dgamma(0.001, 0.001);
vary[i] <- 1.0/precy[i];
}
# calculates area under the curve
la1 <- beta[2]/sqrt(vary[1]); # ROC curve parameters
la2 <- vary[2]/vary[1];
auc <- phi(la1/sqrt(1+la2));
}
# Diabetes data
list(K=2, P=2, N=78, y=c(123,129,115,131,119,111,129,127,118,111,
131,118,126,130,122,112,122,128,123,119,132,118,126,136,118,122,
119,117,129,120,125,115,131,123,130,113,128,138,119,118,124,127,
139,120,122,120,114,114,122,127,123,118,131,130,139,125,135,121,124,
109,106,100,88,106,108,110,111,112,94,122,110,113,106,114,101,99,128,106),
d=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,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,0,0,0,0,0,0,0,0))
# Ordinal values from mammography
list( K=2,P=2,N=70,
y=c(1,1,
2,
3,3,3,3,3,3,3,
4,4,4,4,4,4,4,4,4,4,4,4,
5,5,5,5,5,5,5,5,5,5,5,5,5,
1,1,1,1,1,1,1,1,1,1,
2,2,2,
3,3,3,3,3,3,3,3,3,3,3,3,
4,4,4,4,4,4,4,4,4,
5),
d=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,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,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
list(beta=c(0,0),precy=c(1,1))
BUGS CODE 4.4
Model;
# Clustered Data ROC Area
{
# Joint Distribution of all cell parameters
# abnormal ROI
g11~dgamma(a11,2)
g12~dgamma(a12,2)
g13~dgamma(a13,2)
g14~dgamma(a14,2)
# normal ROI
g01~dgamma(a01,2)
g02~dgamma(a02,2)
g03~dgamma(a03,2)
g04~dgamma(a04,2)
sg<-g11+g12+g13+g14+g01+g02+g03+g04
the1<-g11/sg
the2<-g12/sg
the3<-g13/sg
the4<-g14/sg
p1<- g01/sg
p2<- g02/sg
p3<- g03/sg
p4<- g04/sg
# sum of the the's
sthe<-the1+the2+the3+the4
# sum of the p's
sp<-p1+p2+p3+p4
# truncated distribution of thetas
theta1<-the1/sthe
theta2<-the2/sthe
theta3<-the3/sthe
theta4<-the4/sthe
# truncated distributions of the phi's
ph1<-p1/sp
ph2<-p2/sp
ph3<-p3/sp
ph4<-p4/sp
# area based on truncated
Area<-A1+A2/2
A1<-theta2*ph1+theta3*(ph1+ph2)+theta4*(ph1+ph2+ph3)
A2<-theta1*ph1+theta2*ph2+theta3*ph3+theta4*ph4
}
# a1j's are cell counts for abnormal ROI
# a0j's are cell counts of normal ROI
# The list below is for the clustered data of the lung cancer study.
list(a11=5,a12=27,a13=34,a14=22,a01=161,a02=29, a03=3,a04=1)
# initial values
list( g11=1,g12=1,g13=1,g14=1, g01=1,g02=1,g03=1,g04=1)
BUGS CODE 4.5
Model;
# Comparing modalities 2 by 2 table
{
# Dirichlet distribution generated
g00~dgamma(a00,2)
g01~dgamma(a01,2)
g10~dgamma(a10,2)
g11~dgamma(a11,2)
sumall<-g00+g01+g10+g11+
h00+h01+h10+h11
h00~dgamma(b00,2)
h01~dgamma(b01,2)
h10~dgamma(b10,2)
h11~dgamma(b11,2)
# cell probabilities for diseased
theta00<-g00/sumall
theta01<-g01/sumall
theta10<-g10/sumall
theta11<-g11/sumall
# cell probabilities for non diseased
ph00<-h00/sumall
ph01<-h01/sumall
ph10<-h10/sumall
ph11<-h11/sumall
# truncated distributions
stheta<-theta00+theta01+theta10+theta11
sph<-ph00+ph01+ph10+ph11
# truncated thetas
th00<-theta00/stheta
th01<-theta01/stheta
th10<-theta10/stheta
th11<-theta11/stheta
# truncated phi's
p00<-ph00/sph
p01<-ph01/sph
p10<-ph10/sph
p11<-ph11/sph
# true positives
# a designates row
# b designates column
tpfa<-th10+th11
tpfb<-th01+th11
# false positive
fpfa<-p10+p11
fpfb<-p01+p11
# ratio of tpf a to b
rtpf<-tpfa/tpfb
#ratio of fpf a to b
rfpf<-fpfa/fpfb
}
# CASS data set
# see Pepe for a description of CASS study
#Comparing est and cph
# a is est
# b is cph
list( a00=26,a01=184,a10=30,a11=787,
b00=152,b01=177,b10=47,b11=70)
# the initial values
list(g00=1,g01=1,g10=1,g11=1,
h00=1,h01=1,h10=1,h11=1)
BUGS CODE 4.6
Model;
# Comparing two modalities
# 2 ROC areas
{
# Dirichlet distribution
# modality a diseased
g11~dgamma(a11,2)
g12~dgamma(a12,2)
g13~dgamma(a13,2)
g14~dgamma(a14,2)
g15~dgamma(a15,2)
# modaltiy b diseased
g21~dgamma(a21,2)
g22~dgamma(a22,2)
g23~dgamma(a23,2)
g24~dgamma(a24,2)
g25~dgamma(a25,2)
sg<-g11+g12+g13+g14+g15+
g21+g22+g23+g24+g25
# modality a non diseased
h11~dgamma(b11,2)
h12~dgamma(b12,2)
h13~dgamma(b13,2)
h14~dgamma(b14,2)
h15~dgamma(b15,2)
# modality b non diseased
h21~dgamma(b21,2)
h22~dgamma(b22,2)
h23~dgamma(b23,2)
h24~dgamma(b24,2)
h25~dgamma(b25,2)
sh<-h11+h12+h13+h14+h15+
h21+h22+h23+h24+h25
# the thetas and phis have a Dirichlet distribution
theta11<-g11/sg
theta12<-g12/sg
theta13<-g13/sg
theta14<-g14/sg
theta15<-g15/sg
theta21<-g21/sg
theta22<-g22/sg
theta23<-g23/sg
theta24<-g24/sg
theta25<-g25/sg
ph11<-h11/sh
ph12<-h12/sh
ph13<-h13/sh
ph14<-h14/sh
ph15<-h15/sh
ph21<-h21/sh
ph22<-h22/sh
ph23<-h23/sh
ph24<-h24/sh
ph25<-h25/sh
# truncate modality a (diseased)
sumad<-theta11+theta12+theta13+theta14+theta15
th11<-theta11/sumad
th12<-theta12/sumad
th13<-theta13/sumad
th14<-theta14/sumad
th15<-theta15/sumad
# truncate modality b (diseased)
sumbd<-theta21+theta22+theta23+theta24+theta25
th21<-theta21/sumbd
th22<-theta22/sumbd
th23<-theta23/sumbd
th24<-theta24/sumbd
th25<-theta25/sumbd
# truncate modality a (non diseased)
sumand<-ph11+ph12+ph13+ph14+ph15
p11<-ph11/sumand
p12<-ph12/sumand
p13<-ph13/sumand
p14<-ph14/sumand
p15<-ph15/sumand
# truncate modality b (non diseased)
sumbnd<-ph21+ph22+ph23+ph24+ph25
p21<-ph21/sumbnd
p22<-ph22/sumbnd
p23<-ph23/sumbnd
p24<-ph24/sumbnd
p25<-ph25/sumbnd
# area for a
areaa<-areaa1+areaa2/2
areaa1<- th12*p11+th13*(p11+p12)+th14*(p11+p12+p13)+
th15*(p11+p12+p13+p14)
areaa2<-th11*p11+th12*p12+th13*p13+th14*p14+
th15*p15
# area for b
areab<-areab1+areab2/2
areab1<- th22*p21+th23*(p21+p22)+th24*(p21+p22+p23)+
th25*(p21+p22+p23+p24)
areab2<-th21*p21+th22*p22+th23*p23+th24*p24+
th25*p25
# compares areas a and b
diffarea<-areaa-areab
}
# Compares mri and ct for lung cancer lesions
# hypothetical example
# mri is modlaity a
# ct is modality b
list(a11=1,a12=4,a13=13,a14=11,a15=26,
a21=2,a22=5,a23=19,a24=12,a25=27,
b11=31,b12=19,b13=21,b14=9,b15=3,
b21=32,b22=23,b23=16,b24=6,b25=6)
BUGS CODE 4.7
model ;
# Compares two correlated ROC areas
# Compares two biomarkers
{
for(i in 1:N) {
# yt is log of y
yt[i]<-log(y[i])
yt[i]~dnorm(mu[i],precy[d[i]+1]);
# the regression of log y on logz
mu[i] <- beta[1] + beta[2]*d[i]+beta[3]*zt[i]
# zt is log of z
zt[i]<-log(z[i])
zt[i]~dnorm(vu[i],precz[d[i]+1])
# the regression of logz on logy
vu[i]<-delta[1]+delta[2]*d[i]+delta[3]*yt[i]
}
# prior distributions - non-informative prior;
for ( i in 1:3) { beta[i] ~ dnorm(0,0.0001)
delta[i]~dnorm(0,0.0001)
}
for(i in 1:K) {
precy[i]~dgamma(0.000001, 0.000001);
vary[i] <- 1.0/precy[i];
precz[i]~dgamma(0.000001,0.000001);
varz[i]<-1.0/precz[i];
}
# calculates area under the curve
la1y <- beta[2]/sqrt(vary[1]);
la2y <- vary[2]/vary[1];
aucy <- phi(la1y/sqrt(1+la2y));
la1z <- delta[2]/sqrt(varz[1]);
la2z <- varz[2]/varz[1];
aucz <- phi(la1z/sqrt(1+la2z));
diff<-aucy-aucz
}
# Wieand et al.
# two biomarkers for pancreatic cancer
# see Zhou et al.
# y is CA19-9 and z is CA125
list(K=2, N=141,
y=c(28.00,15.50,8.20,3.40,17.30,15.20,32.90,11.10,87.50,16.20,107.90,5.70,25.60,31.20,21.60,55.60,8.80,6.50,22.10,14.40,44.20,3.70,7.80,8.90,18.00,6.50,4.90,10.40,5.00,5.30,6.50,6.90,8.20,21.80,6.60,7.60,15.40,59.20,5.10,10.00,5.30,32.60,4.60,6.90,4.00,3.65,7.80,32.50,11.50,4.00,10.20,2.40,719.00,2106.67,24000.00,1715.00,3.60,521.50,1600.00,454.00,109.70,23.70,464.00,9810.00,255.00,58.70,225.00,90.10,50.00,5.60,4070.00,592.00,28.60,6160.00,1090.00,10.40,27.30,162.00,3560.00,14.70,83.30,336.00,55.70,1520.00,3.90,5.80,8.45,361.00,369.00,8230.00,39.30,43.50,361.00,12.80,18.00,9590.00,555.00,60.20,21.80,900.00,6.60,239.00,3100.00,3275.00,682.00,85.40,10290.00,770.00,247.60,12320.00,113.10,1079.00,45.60,1630.00,79.40,508.00,3190.00,542.00,1021.00,235.00,251.00,3160.00,479.00,222.00,15.70,2540.00,11630.00,1810.00,6.90,4.10,15.60,9820.00,1490.00,15.70,45.80,7.80,12.80,100.53,227.00,70.90,2500.00),
z=c(13.30,11.10,16.70,12.60,7.40,5.50,32.10,27.20,6.60,9.80,10.50,7.80,9.10,12.30,12.00,42.10,5.90,9.20,7.30,6.80,10.70,15.70,8.00,6.80,47.35,17.90,96.20,108.90,16.60,9.50,179.00,12.10,35.60,15.00,12.60,5.90,10.10,8.50,11.40,54.65,9.70,11.20,35.70,22.50,21.20,5.60,9.40,12.00,9.80,17.20,10.60,79.10,31.40,15.00,77.80,25.70,11.70,8.25,14.95,8.70,14.10,123.90,12.10,99.10,18.60,10.50,6.60,74.00,43.90,45.70,13.00,7.30,8.60,17.20,15.40,14.30,93.10,66.30,26.70,32.40,9.90,30.30,11.20,202.00,35.70,9.20,103.60,21.40,8.10,29.90,17.50,30.80,57.30,6.50,33.80,53.60,17.20,94.20,33.50,3.70,11.70,19.90,38.70,27.30,20.10,86.10,844.00,36.90,6.90,27.70,9.90,38.60,142.60,12.50,11.60,21.20,13.20,19.20,1024.00,14.10,34.80,35.30,35.00,15.50,12.10,31.60,184.80,24.80,10.40,34.50,19.40,22.20,53.90,15.40,17.30,36.80,49.80,26.57,9.70,19.20,14.20),
d= c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1))
list(beta = c( 0,0,0),delta=c(0,0,0), precy = c(1,1),precz=c(1,1))
No comments:
Post a Comment