BUGS CODE 9.1
model;
# Binormal model for area under ROC
# Calculates posterior distribution of model parameters and the area under curve.
# calculates the coordinates of the optimal psa threshold
# yt= log transformed outcome, d=disease status, # var y[N], yt[N], d[N], mu[N], beta[P], vary[K], precy[K], auc, la1, la2;
{
# likelihood function
for(i in 1:N) {
y[i]~dnorm(mu[i],precy[d[i]+1]);
#yt[i] <- log(y[i]);
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
a <- beta[2]/sqrt(vary[2]);
# ROC curve parameters
la2 <- vary[1]/vary[2];
# ROC area
auc <- phi(a/sqrt(1+la2));
b<- sqrt(la2) ;
ka<-R*(1-p)/p
nfpf<- a*b-sqrt(a*a+2*(1-b*b)*log(ka/b));
dfpf<-1-b*b;
fpf <- phi(nfpf/dfpf) ;
ntpf<- a - b*sqrt(a*a+(1-b*b)*log(ka/b));
tpf<- phi(ntpf/dfpf);
}
# psa data
list(K=2, P=2, N=683, p=.5, R=1,
y=c(.03,
.09,.23,.27,.27,.29,.29,.29,.30,.31,.33,.35,.37,.37,.42,.43,.44,.45,
.45,.46,.46,.47,.47,.48,.49,.49,.50,.50,.50,.51,.51,.55,.55,.56,.57,
.57,.58,.58,.58,.58,.59,.59,.59,.61,.61,.62,.62,.63,.63,.64,.64,.64,
.64,.65,.65,.65,.66,.66,.66,.66,.66,.66,.67,.67,.67,.67,.67,.68,.68,
.69,.69,.69,.69,.69,.70,.71,.72,.72,.73,.74,.74,.75,.75,.75,.75,.75,
.76,.76,.77,.77,.77,.77,.77,.77,.78,.78,.78,.78,.78,.78,.79,.79,.79,
.79,.80,.80,.80,.81,.81,.81,.81,.82,.83,.83,.84,.85,.86,.87,.87,.87,
.87,.87,.88,.89,.89,.89,.89,.89,.92,.92,.92,.93,.93,.93,.93,.93,.93,
.94,.94,.95,.95,.95,.95,.96,.96,.97,.97,.98,.98,.98,.98,.98,.99,1.00,1.00,
1.00,1.01,1.01,1.02,1.03,1.03,1.03,1.03,1.03,1.03,1.04,1.04,1.04,
1.04,1.04,1.05,1.05,1.05,1.05,1.06,1.06,1.06,1.06,1.07,1.07,1.07,
1.08,1.08,1.08,1.11,1.11,1.12,1.12,1.13,1.13,1.13,1.14,1.15,1.15,
1.15,1.15,1.15,1.15,1.15,1.15,1.16,1.16,1.16,1.17,1.17,1.17,1.17,
1.18,1.18,1.18,1.18,1.18,1.19,1.19,1.19,1.20,1.20,1.21,1.22,1.22,
1.22,1.23,1.23,1.24,1.24,1.24,1.25,1.25,1.25,1.25,1.25,1.25,1.26,
1.26,1.26,1.27,1.27,1.27,1.27,1.27,1.27,1.28,1.28,1.29,1.30,1.30,
1.31,1.31,1.32,1.32,1.33,1.34,1.35,1.35,1.35,1.35,1.35,1.35,1.35,
1.36,1.37,1.37,1.37,1.38,1.39,1.39,1.40,1.40,1.40,1.40,1.41,1.41,
1.41,1.41,1.41,1.41,1.43,1.43,1.43,1.43,1.44,1.44,1.45,1.46,1.46,
1.47,1.47,1.47,1.48,1.48,1.49,1.49,1.50,1.50,1.50,1.50,1.51,
1.51,1.51,1.51,1.53,1.54,1.54,1.55,1.55,1.56,1.57,1.57,1.58,1.58,
1.58,1.61,1.62,1.62,1.62,1.62,1.64,1.67,1.67,1.67,1.67,1.67,1.68,
1.69,1.69,1.70,1.70,1.70,1.71,1.71,1.71,1.71,1.71,1.71,1.71,1.71,
1.73,1.73,1.73,1.74,1.79,1.80,1.80,1.83,1.85,1.85,1.88,1.88,1.88,
1.89,1.89,1.89,1.91,1.91,1.91,1.92,1.93,1.93,1.94,1.95,1.96,2.01,
2.01,2.03,2.03,2.03,2.04,2.04,2.05,2.05,2.06,2.07,2.08,2.08,2.10,
2.11,2.13,2.13,2.14,2.16,2.17,2.19,2.19,2.19,2.22,2.22,2.23,2.24,
2.27,2.27,2.27,2.28,2.28,2.29,2.29,2.30,2.30,2.33,2.34,2.34,2.35,
2.36,2.36,2.37,2.40,2.41,2.42,2.43,2.43,2.43,2.43,2.46,2.50,2.50,
2.51,2.51,2.52,2.53,2.55,2.55,2.56,2.56,2.57,2.58,2.61,2.62,2.62,
2.63,2.63,2.63,2.66,2.69,2.70,2.71,2.73,2.77,2.79,2.82,2.82,2.82,
2.83,2.84,2.84,2.85,2.86,2.86,2.87,2.88,2.88,2.90,2.92,2.92,2.93,
2.95,2.96,2.96,2.96,2.97,2.98,3.03,3.03,3.04,3.05,3.05,3.08,3.10,
3.11,3.13,3.17,3.17,3.18,3.20,3.21,3.24,3.25,3.25,3.29,3.30,3.30,
3.32,3.32,3.33,3.34,3.35,3.38,3.41,3.42,3.43,3.45,3.51,3.55,3.57,
3.57,3.58,3.58,3.61,3.65,3.65,3.66,3.68,3.69,3.70,3.73,3.77,3.78,
3.78,3.78,3.80,3.84,3.88,3.89,3.95,3.97,3.97,4.00,4.03,4.03,4.04,
4.05,4.08,4.12,4.15,4.19,4.20,4.20,4.20,4.30,4.33,4.34,4.38,4.39,
4.40,4.41,4.44,4.47,4.47,4.48,4.52,4.54,4.60,4.62,4.64,4.70,4.75,
4.75,4.76,4.78,4.90,4.90,4.93,4.94,4.98,5.02,5.09,5.10,5.11,5.12,
5.13,5.13,5.25,5.28,5.37,5.39,5.44,5.44,5.53,5.54,5.64,5.65,5.67,
5.73,5.75,5.81,5.85,6.07,6.07,6.16,6.18,6.27,6.29,6.31,6.41,6.48,6.48,
6.50,6.52,6.52,6.54,6.54,6.56,6.56,6.77,6.92,6.93,7.09,7.19,7.21,
7.23,7.24,7.28,7.29,7.42,7.43,7.53,7.59,7.64,7.78,7.90,8.04,8.15,
8.31,8.37,8.57,8.62,8.69,9.07,9.11,9.15,9.15,9.17,9.24,9.30,9.33,
9.76,9.94,9.96,9.97,10.11,10.60,10.71,10.92,11.33,11.40,11.54,
11.62,11.65,12.69,12.69,13.61,13.94,14.82,15.41,15.84,15.84,15.89,
16.18,16.48,16.70,16.81,17.10,17.17,17.57,19.35,20.10,20.24,20.47,
20.53,21.48,22.50,23.81,24.63,25.06,26.67,27.68,29.31,31.46,33.02,
35.93,37.63,37.66,38.39,43.30,48.80,49.16,51.72,61.16,72.07,79.21,
90.66,99.97,99.98,99.98,99.98),
d=c(.00,1.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,
1.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,
.00,.00,.00,.00,.00,.00,1.00,.00,.00,1.00,.00,.00,.00,.00,.00,.00,.00,
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,1.00,
.00,1.00,.00,.00,.00,.00,.00,.00,1.00,.00,.00,.00,.00,.00,.00,.00,
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,
.00,.00,1.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,
.00,.00,.00,.00,.00,.00,.00,.00,1.00,1.00,.00,.00,1.00,.00,.00,.00,
.00,1.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,
.00,.00,1.00,.00,.00,.00,.00,.00,.00,.00,1.00,1.00,.00,.00,1.00,.00,
1.00,1.00,.00,.00,.00,.00,.00,.00,.00,1.00,.00,.00,.00,.00,.00,.00,
.00,.00,1.00,.00,.00,.00,.00,.00,1.00,.00,.00,.00,.00,.00,.00,.00,
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,1.00,.00,.00,.00,.00,.00,.00,
.00,.00,.00,.00,1.00,.00,.00,.00,.00,1.00,1.00,.00,.00,1.00,1.00,.00,
.00,.00,.00,.00,.00,.00,.00,.00,1.00,.00,.00,.00,.00,1.00,1.00,.00,.00,
.00,.00,.00,.00,.00,.00,.00,.00,.00,1.00,.00,.00,1.00,.00,.00,.00,.00,
.00,.00,.00,.00,.00,.00,1.00,.00,.00,.00,.00,.00,1.00,.00,.00,.00,.00,
.00,.00,1.00,1.00,.00,.00,1.00,.00,.00,.00,.00,1.00,.00,1.00,.00,.00,
.00,.00,1.00,1.00,.00,.00,.00,.00,.00,1.00,.00,1.00,.00,.00,.00,.00,.00,
.00,.00,.00,.00,1.00,1.00,1.00,.00,1.00,.00,.00,.00,1.00,1.00,.00,.00,.00,
.00,1.00,.00,.00,.00,.00,.00,1.00,1.00,.00,1.00,.00,1.00,.00,1.00,.00,
.00,.00,1.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,1.00,1.00,.00,1.00,.00,
1.00,.00,1.00,.00,1.00,.00,.00,1.00,.00,1.00,.00,.00,.00,1.00,1.00,.00,.00,
.00,.00,.00,1.00,1.00,.00,.00,.00,1.00,.00,1.00,1.00,.00,1.00,1.00,.00,
1.00,.00,.00,.00,1.00,.00,.00,.00,.00,1.00,1.00,.00,1.00,.00,.00,.00,
1.00,.00,1.00,.00,.00,1.00,.00,1.00,.00,.00,.00,1.00,.00,1.00,1.00,.00,
1.00,.00,.00,1.00,1.00,.00,.00,1.00,.00,1.00,.00,.00,.00,.00,.00,.00,
.00,1.00,1.00,.00,1.00,1.00,.00,1.00,.00,1.00,1.00,.00,.00,.00,1.00,1.00,
.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,.00,1.00,1.00,1.00,
1.00,.00,1.00,1.00,1.00,1.00,.00,.00,.00,.00,1.00,1.00,1.00,1.00,.00,
1.00,1.00,.00,1.00,.00,1.00,1.00,1.00,.00,1.00,.00,.00,.00,1.00,1.00,
.00,1.00,1.00,1.00,1.00,.00,1.00,.00,1.00,1.00,1.00,.00,1.00,1.00,.00,
.00,.00,1.00,.00,1.00,1.00,1.00,.00,1.00,.00,1.00,.00,1.00,1.00,.00,
1.00,1.00,1.00,.00,1.00,.00,1.00,.00,1.00,.00,.00,1.00,1.00,1.00,.00,
1.00,1.00,1.00,.00,1.00,.00,1.00,1.00,1.00,1.00,.00,1.00,1.00,.00,1.00,
1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,
1.00,1.00,1.00,1.00,.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,
1.00,1.00,1.00,1.00,1.00,.00,1.00,1.00,.00,1.00,.00,1.00,1.00,1.00,
1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,.00,1.00,1.00,
1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00))
# diabetes data
list( N= 78,K=2,P=2, p=.5,R=1,y=c(123,129,115,131,119,111,129,127,118,111,131,118,126,
130,122,112,122,128,123,119,119,132,118,126,136,118,122,
119,117,129,120,125,115,131,123,130,113,128,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))
list(beta=c(0,0),precy=c(1,1))
BUGS CODE 9.2.
model;
# Pepe Caret PSA Etzioni et al.
# Uniform prior
{
# This is the calculation for the critical benefit
# method of Vickers
# uniform prior
# pt is threshold probability
# the paramaeters for the gamma variables are determined seperately by logistic regression
# pt=.01
g01[1,1]~dgamma(1,2)
g01[1,2]~dgamma(1,2)
g01[2,1]~dgamma(455,2)
g01[2,2]~dgamma(230,2)
h01<- sum(g01[,])
for( i in 1 :2 ) {for( j in 1 :2 ){ theta01[i,j]<-g01[i,j]/h01}}
cb01<-theta01[2,2]-theta01[2,1]*(.01/.99)
#pt=.05
g05[1,1]~dgamma(1,2)
g05[1,2]~dgamma(1,2)
g05[2,1]~dgamma(455,2)
g05[2,2]~dgamma(230,2)
h05<- sum(g05[,])
for( i in 1 :2 ) {for( j in 1 :2 ){ theta05[i,j]<-g05[i,j]/h05}}
cb05<-theta05[2,2]-theta05[2,1]*(.05/.95)
#pt=.1
g1[1,1]~dgamma(1,2)
g1[1,2]~dgamma(1,2)
g1[2,1]~dgamma(455,2)
g1[2,2]~dgamma(230,2)
h1<- sum(g1[,])
for( i in 1 :2 ) {for( j in 1 :2 ){ theta1[i,j]<-g1[i,j]/h1}}
cb1<-theta1[2,2]-theta1[2,1]*(1/9)
#pt=.15
g15[1,1]~dgamma(16,2)
g15[1,2]~dgamma(2,2)
g15[2,1]~dgamma(440,2)
g15[2,2]~dgamma(229,2)
h15<- sum(g15[,])
for( i in 1 :2 ) {for( j in 1 :2 ){ theta15[i,j]<-g15[i,j]/h15}}
cb15<-theta15[2,2]-theta15[2,1]*(15/85)
# pt=.2
g2[1,1]~dgamma(252,2)
g2[1,2]~dgamma(26,2)
g2[2,1]~dgamma(204,2)
g2[2,2]~dgamma(205,2)
h2<- sum(g2[,])
for( i in 1 :2 ) {for( j in 1 :2 ){ theta2[i,j]<-g2[i,j]/h2}}
cb2<-theta2[2,2]-theta2[2,1]/4
# pt=.25
g25[1,1]~dgamma(338,2)
g25[1,2]~dgamma(51,2)
g25[2,1]~dgamma(118,2)
g25[2,2]~dgamma(180,2)
h25<- sum(g25[,])
for( i in 1 :2 ) {for( j in 1 :2 ){ theta25[i,j]<-g25[i,j]/h25 }}
cb25<-theta25[2,2]-theta25[2,1]/3
# pt=.3
g3[1,1]~dgamma(379,2)
g3[1,2]~dgamma(76,2)
g3[2,1]~dgamma(77,2)
g3[2,2]~dgamma(155,2)
h3<- sum(g3[,])
for( i in 1 :2 ) {for( j in 1 :2 ){ theta3[i,j]<-g3[i,j]/h3 }}
cb3<-theta3[2,2]-theta3[2,1]*3/7
# pt=.35
g35[1,1]~dgamma(406,2)
g35[1,2]~dgamma(91,2)
g35[2,1]~dgamma(50,2)
g35[2,2]~dgamma(140,2)
h35<- sum(g35[,])
for( i in 1 :2 ) {for( j in 1 :2 ){ theta35[i,j]<-g35[i,j]/h35 }}
cb35<-theta35[2,2]-theta35[2,1]*35/65
# pt=.4
g4[1,1]~dgamma(416,2)
g4[1,2]~dgamma(109,2)
g4[2,1]~dgamma(40,2)
g4[2,2]~dgamma(122,2)
h4<- sum(g4[,])
for( i in 1 :2 ) {for( j in 1 :2 ){ theta4[i,j]<-g4[i,j]/h4 }}
cb4<-theta4[2,2]-theta4[2,1]*2/3
# differences in clinical benefit for a minus that for scenario b
d01<- cb01-clinb01
d05<-cb05-clinb05
d1 <- cb1-clinb1
d15<-cb15-clinb15
d2 <-cb2-clinb2
d3 <-cb3-clinb3
d25<- cb25-clinb25
d35<-cb35-clinb35
d4 <- cb4-clinb4
# probability that clinical benefit of a is greater than that for scenario b
p01<-step(d01)
p05<-step(d05)
p1<- step(d1)
p15<-step(d15)
p2<- step(d2)
p25<-step(d25)
p3<- step(d3)
p35<-step(d35)
p4<-step(d4)
p<- p01*p05*p1*p15*p2*p25*p3*p35*p4
# The following is the clinical benefit where all patients are biopsed and based on the
# Pepe total PSA data and the method of Vickers
# scenario b
# pt=.01
clinb01<- (theta01[1,2]+theta01[2,2]) -(1-(theta01[1,2]+theta01[2,2]))*(.01/.99)
# pt=.-5
clinb05<- (theta05[1,2]+theta05[2,2]) -(1-(theta05[1,2]+theta05[2,2]))*(.05/.95)
# pt=.1
clinb1<- (theta1[1,2]+theta1[2,2]) -(1-(theta1[1,2]+theta1[2,2]))*(1/9)
# pt=.15
clinb15<- (theta15[1,2]+theta15[2,2]) -(1-(theta15[1,2]+theta15[2,2]))*(15/85)
# pt = .2
clinb2<- (theta2[1,2]+theta2[2,2]) -(1-(theta2[1,2]+theta2[2,2]))*(1/4)
# pt = .25
clinb25<- (theta25[1,2]+theta25[2,2]) -(1-(theta25[1,2]+theta25[2,2]))*(1/3)
# pt=.3
clinb3<- (theta3[1,2]+theta3[2,2]) -(1-(theta3[1,2]+theta3[2,2]))*(3/7)
# pt=.35
clinb35<- (theta35[1,2]+theta35[2,2]) -(1-(theta35[1,2]+theta35[2,2]))*(35/65)
#pt=.4
clinb4<- (theta4[1,2]+theta4[2,2]) -(1-(theta4[1,2]+theta4[2,2]))*(2/3)
}
# data appears in above statements
# initial values are initiated with specification tool
BUGS CODE 9.3
model {
for( i in 1 : N ) {
y[i] ~ dbern(p[i])
logit( p[i]) <- beta[1] + beta[2]*n[i]+beta[3]*r[i]
}
phat <- mean(p[])
for (i in 1:3 ){
beta[i] ~ dnorm(0.0,0.0001)}
}
No comments:
Post a Comment