BUGS CODE 6.1
model;
{
# reader 2
# ROC area
# melanoma example
# code is from Congdon (2003, page 1020)
# non diseased
for(i in 1:30){for(j in 1:5){logit(ndgamma[i,j])<-ndtheta[j]-ndmu[i]}}
for(i in 1:30){ndp[i,1]<-ndgamma[i,1]}
for(i in 1:30){ndp[i,2]<-ndgamma[i,2]-
ndgamma[i,1]}
for(i in 1:30){ndp[i,3]<-ndgamma[i,3]-
ndgamma[i,2]}
for(i in 1:30){ndp[i,4]<-ndgamma[i,4]-
ndgamma[i,3]}
for(i in 1:30){ndy[i]~dcat(ndp[i,1:5])}
for(i in 1:30){ndp[i,5]<-1-ndgamma[i,4]}
# intercept depends on y
for(i in 1:30){
ndmu[i]<-ndb0[ndy[i]]}
for(i in 1:5){ndb0[i]~dnorm(0,.001)}
ndtheta[1]~dnorm(0,1)
ndtheta[2]~dnorm(0,1)
ndtheta[3]~dnorm(0,1)
ndtheta[4]~dnorm(0,1)
ndtheta[5]~dnorm(0,1)I(ndtheta[4],)
for( i in 1:5){ndq[i]<-mean(ndp[,i])}
# diseased population
for(i in 1:50){for(j in 1:5){logit(dgamma[i,j])<-dtheta[j]-dmu[i]}}
for(i in 1:50){dp[i,1]<-dgamma[i,1]}
for(i in 1:50){dp[i,2]<-dgamma[i,2]-
dgamma[i,1]}
for(i in 1:50){dp[i,3]<-dgamma[i,3]-
dgamma[i,2]}
for(i in 1:50){dp[i,4]<-dgamma[i,4]-
dgamma[i,3]}
for(i in 1:50){dy[i]~dcat(dp[i,1:5])}
for(i in 1:50){dp[i,5]<-1-dgamma[i,4]}
# intercept depends on y
for(i in 1:50){
dmu[i]<-db0[dy[i]]}
for(i in 1:5){db0[i]~dnorm(0,.001)}
dtheta[1]~dnorm(0,1)
dtheta[2]~dnorm(0,1)
dtheta[3]~dnorm(0,1)
dtheta[4]~dnorm(0,1)
dtheta[5]~dnorm(0,1)I(dtheta[4],)
for( i in 1:5){dq[i]<-mean(dp[,i])}
# roc area
area<-a1+a2/2
a1<-dq[2]*ndq[1]+dq[3]*(ndq[1]+ndq[2])+
dq[4]*(ndq[1]+ndq[2]+ndq[3])+
dq[5]*(ndq[1]+ndq[2]+ndq[3]+ndq[4])
a2<-dq[1]*ndq[1]+dq[2]*ndq[2]+dq[3]*ndq[3]+
dq[4]*ndq[4]+dq[5]*ndq[5]
}
list( ndy=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
2,2,2,2,2,2,2,
3,3,3,3,
4,4,4,
5 ),
# data for diseased
dy=c(1,1, 2,2,2,2,2,2,2,2,
3,3,3,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,5,5,5,5,5))
list(ndtheta=c(0,0,0,0,0),dtheta=c(0,0,0,0,0))
BUGS CODE 6.2
model;
# Calculates posterior distribution of model parameters and the area under curve.
# type 2 diabetes with three readers
{
# likelihood function
for(i in 1:N) {
g1[i]~dnorm(mu1[i],precy1[d[i]+1]);
mu1[i] <- beta1[1] + beta1[2]*d[i] + beta1[3]*age[i]+beta1[4]*male[i] ;
g2[i]~dnorm(mu2[i],precy2[d[i]+1]);
mu2[i] <- beta2[1] + beta2[2]*d[i] + beta2[3]*age[i]+beta2[4]*male[i] ;
g3[i]~dnorm(mu3[i],precy3[d[i]+1]);
mu3[i] <- beta3[1] + beta3[2]*d[i] + beta3[3]*age[i]+beta3[4]*male[i] ;
}
# prior distributions - non-informative prior; similarly for informative priors
for(i in 1:P) {
beta1[i] ~ dnorm(0, 0.001);
beta2[i] ~ dnorm(0, 0.001);
beta3[i] ~ dnorm(0, 0.001);
}
for(i in 1:K) {
precy1[i]~dgamma(0.001, 0.001);
vary1[i] <- 1.0/precy1[i];
precy2[i]~dgamma(0.001, 0.001);
vary2[i] <- 1.0/precy2[i];
precy3[i]~dgamma(0.001, 0.001);
vary3[i] <- 1.0/precy3[i];
}
# calculates area under the curve
# reader 1
la11 <- beta1[2]/sqrt(vary1[1]); # ROC curve parameters
la21 <- vary1[2]/vary1[1];
auc1 <- phi(la11/sqrt(1+la21));
#reader 2
la12 <- beta2[2]/sqrt(vary2[1]); # ROC curve parameters
la22 <- vary2[2]/vary2[1];
auc2 <- phi(la12/sqrt(1+la22));
# reader 3
la13 <- beta3[2]/sqrt(vary3[1]); # ROC curve parameters
la23 <- vary3[2]/vary3[1];
auc3 <- phi(la13/sqrt(1+la23));
}
list(P=4,N=307, K=2,
# g1 is blood glucose levels for reader 1
g1=c(96.63,102.98,97.72,97.82,106.94,105.32,92.61,94.99,
105.49,97.34,97.72,96.87,98.53,102.57,99.69,96.46,93.68,
97.46,104.60,98.49,107.34,96.03,105.17,96.87,98.16,
104.14,99.73,94.68,93.85,99.70,95.07,99.74,102.22,98.99,
103.72,101.55,101.55,95.54,97.47,103.37,100.31,100.55,
99.76,103.12,92.16,106.42,102.03,96.97,103.79,96.58,
113.96,100.26,95.07,104.00,101.47,105.84,103.61,98.03,
93.45,92.92,98.48,100.14,97.46,97.88,104.21,92.92,
104.49,95.51,100.49,99.46,105.03,91.78,100.75,105.68,
100.31,91.27,103.92,98.78,92.80,107.75,104.85,104.24,
93.57,100.69,97.11,101.41,84.43,101.88,94.94,94.91,
100.04,104.18,104.81,98.06,107.01,94.13,99.19,98.87,
99.01,96.42,103.26,109.30,97.20,94.74,103.36,103.82,
93.54,97.27,96.29,100.58,102.62,94.51,101.84,98.10,
102.66,99.73,96.50,104.86,100.69,97.57,101.81,98.88,
101.00,100.48,98.99,108.75,105.34,108.13,100.90,105.06,
98.10,106.16,105.64,94.18,104.07,98.64,97.82,98.49,
100.74,100.63,93.91,94.89,103.31,102.42,98.5,196.68,
109.31,95.59,99.23,102.60,104.24,103.14,109.07,103.23,
103.72,98.41,93.53,92.92,101.26,98.75,106.58,94.80,
102.49,101.80,99.97,97.73,106.66,100.91,93.13,105.04,
101.92,91.52,107.76,94.59,97.97,98.59,104.58,107.60,
98.14,101.84,101.41,92.35,99.41,99.63,96.51,100.77,
100.67,93.19,103.83,108.11,96.35,106.37,99.29,
102.72,89.20,101.92,105.87,96.66,101.85,103.92,101.38,
95.23,99.60,98.08,99.64,111.32,108.37,91.69,95.38,98.09,
92.05,106.36,93.98,102.26,103.81,98.00,99.20,106.46,
109.58,113.86,103.72,105.94,114.61,111.08,106.89,
119.51,110.30,110.00,108.31,108.68,108.98,115.01,
113.07,114.89,109.79,105.70,114.20,113.53,113.97,
110.91,110.33,115.78,111.05,108.53,111.56,110.78,
109.71,112.18,112.05,109.46,103.84,112.23,118.56,
110.60,109.54,112.31,100.78,114.07,112.14,107.85,
111.65,105.94,108.63,109.89,107.14,108.76,110.11,
104.60,107.11,112.49,113.74,103.19,105.07,109.04,
110.45,105.02,108.27,109.17,110.37,110.92,107.53,
109.22,113.01,108.74,116.74,112.10,110.88,111.08,
110.22,111.23,112.94,99.04,113.51,107.26,110.76,
108.06,97.03,109.14,105.56,111.55,108.85,98.46,110.24,
112.22,108.57,105.95,106.30),
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,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,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,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,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),
age=c(41,36,38,46,38,45,42,42,43,43,43,55,47,46,47,47,32,38,45,37,43,36,52,48,37,42,41,47,46,46,45,49,44,48,50,48,43,47,44,57,48,49,40,41,45,45,52,45,48,45,42,46,48,44,36,48,39,44,51,48,47,39,38,43,39,45,40,35,36,41,46,48,55,41,44,35,38,47,45,50,40,44,46,38,38,50,44,40,46,37,43,40,46,43,36,44,32,47,40,38,42,46,45,41,53,45,41,40,55,48,44,47,47,45,41,53,41,38,43,47,50,45,43,46,42,43,42,47,46,44,37,42,43,44,44,46,36,50,40,39,37,55,41,45,43,39,54,37,38,42,44,48,50,33,42,48,40,49,38,47,39,38,47,39,44,49,46,48,38,39,38,48,42,42,43,36,34,41,36,49,43,35,40,46,44,41,49,46,42,47,42,42,49,43,41,47,47,44,39,42,43,51,43,46,37,44,42,38,35,42,45,49,42,40,45,48,42,52,53,49,63,53,62,57,57,64,53,54,55,59,54,53,55,63,52,58,59,57,56,56,55,59,59,62,59,56,64,59,56,60,54,60,54,54,53,57,58,54,59,63,55,59,51,52,57,60,58,50,62,59,61,53,64,50,55,57,60,58,59,56,55,53,57,53,54,59,61,59,56,56,58,57,60,57,59,54,60,51,61,57,53,53,60,64,58,56,63),
male=c(1,0,1,0,0,0,1,1,1,0,0,0,0,0,0,1,0,1,1,1,0,1,0,1,1,1,1,0,1,0,0,1,1,1,0,1,0,1,1,0,1,0,1,0,0,1,0,0,1,0,1,1,0,1,0,0,1,0,0,1,1,1,0,0,1,1,1,0,0,1,1,0,0,0,1,1,0,0,0,0,0,0,0,0,1,0,1,0,1,0,0,0,0,0,0,1,1,1,1,0,1,1,1,1,0,1,1,0,1,0,1,1,0,1,1,0,0,1,1,1,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,0,0,0,0,1,0,0,1,1,1,1,1,1,0,0,0,1,0,0,1,0,0,1,0,1,1,1,0,1,0,1,0,1,1,0,0,1,0,0,1,0,0,0,0,1,0,0,1,0,0,1,1,1,0,0,0,1,1,0,1,1,0,1,0,0,0,0,1,1,0,1,1,1,1,0,1,1,0,0,1,1,1,0,0,0,0,0,0,1,0,1,1,1,0,1,1,1,0,1,1,0,1,1,1,0,0,0,1,0,0,0,0,1,1,1,0,1,0,1,1,0,1,1,1,1,1,1,1,0,0,1,1,0,0,1,1,1,1,1,0,0,0,1,0,0,1,1,1,1,0,0,1,1,0,1,1,0,0,0,0,0,1,0,1,0,1,0,0,1,1,1,0),
# glucose values of reader 2
g2=c(94.63,100.98,95.72,95.82,104.94,103.32,90.61,92.99,103.49,95.34,95.72,94.87,96.53,100.57,97.69,94.46,91.68,95.46,102.60,96.49,105.34,94.03,103.17,94.87,96.16,102.14,97.73,92.68,91.85,97.70,93.07,97.74,100.22,96.99,101.72,99.55,99.55,93.54,95.47,101.37,98.31,98.55,97.76,101.12,90.16,104.42,100.03,94.97,101.79,94.58,111.96,98.26,93.07,102.00,99.47,103.84,101.61,96.03,91.45,90.92,96.48,98.14,95.46,95.88,102.21,90.92,102.49,93.51,98.49,97.46,103.03,89.78,98.75,103.68,98.31,89.27,101.92,96.78,90.80,105.75,102.85,102.24,91.57,98.69,95.11,99.41,82.43,99.88,92.94,92.91,98.04,102.18,102.81,96.06,105.01,92.13,97.19,96.87,97.01,94.42,101.26,107.30,95.20,92.74,101.36,101.82,91.54,95.27,94.29,98.58,100.62,92.51,99.84,96.10,100.66,97.73,94.50,102.86,98.69,95.57,99.81,96.88,99.00,98.48,96.99,106.75,103.34,106.13,98.90,103.06,96.10,104.16,103.64,92.18,102.07,96.64,95.82,96.49,98.74,98.63,91.91,92.89,101.31,100.42,96.51,94.68,107.31,93.59,97.23,100.60,102.24,101.14,107.07,101.23,101.72,96.41,91.53,90.92,99.26,96.75,104.58,92.80,100.49,99.80,97.97,95.73,104.66,98.91,91.13,103.04,99.92,89.52,105.76,92.59,95.97,96.59,102.58,105.60,96.14,99.84,99.41,90.35,97.41,97.63,94.51,98.77,98.67,91.19,101.83,106.11,94.35,104.37,97.29,100.72,87.20,99.92,103.87,94.66,99.85,101.92,99.38,93.23,97.60,96.08,97.64,109.32,106.37,89.69,93.38,96.09,90.05,104.36,91.98,100.26,101.81,96.00,97.20,110.46,113.58,117.86,107.72,109.94,118.61,115.08,110.89,123.51,114.30,114.00,112.31,112.68,112.98,119.01,117.07,118.89,113.79,109.70,118.20,117.53,117.97,114.91,114.33,119.78,115.05,112.53,115.56,114.78,113.71,116.18,116.05,113.46,107.84,116.23,122.56,114.60,113.54,116.31,104.78,118.07,116.14,111.85,115.65,109.94,112.63,113.89,111.14,112.76,114.11,108.60,111.11,116.49,117.74,107.19,109.07,113.04,114.45,109.02,112.27,113.17,114.37,114.92,111.53,113.22,117.01,112.74,120.74,116.10,114.88,115.08,114.22,115.23,116.94,103.04,117.51,111.26,114.76,112.06,101.03,113.14,109.56,115.55,112.85,102.46,114.24,116.22,112.57,109.95,110.30),
# glucose levels of reader 3
g3=c(101.63,107.98,102.72,102.82,111.94,110.32,97.61,99.99,110.49,102.34,102.72,101.87,103.53,107.57,104.69,101.46,98.68,102.46,109.60,103.49,112.34,101.03,110.17,101.87,103.16,109.14,104.73,99.68,98.85,104.70,100.07,104.74,107.22,103.99,108.72,106.55,106.55,100.54,102.47,108.37,105.31,105.55,104.76,108.12,97.16,111.42,107.03,101.97,108.79,101.58,118.96,105.26,100.07,109.00,106.47,110.84,108.61,103.03,98.45,97.92,103.48,105.14,102.46,102.88,109.21,97.92,109.49,100.51,105.49,104.46,110.03,96.78,105.75,110.68,105.31,96.27,108.92,103.78,97.80,112.75,109.85,109.24,98.57,105.69,102.11,106.41,89.43,106.88,99.94,99.91,105.04,109.18,109.81,103.06,112.01,99.13,104.19,103.87,104.01,101.42,108.26,114.30,102.20,99.74,108.36,108.82,98.54,102.27,101.29,105.58,107.62,99.51,106.84,103.10,107.66,104.73,101.50,109.86,105.69,102.57,106.81,103.88,106.00,105.48,103.99,113.75,110.34,113.13,105.90,110.06,103.10,111.16,110.64,99.18,109.07,103.64,102.82,103.49,105.74,105.63,98.91,99.89,108.31,107.42,103.51,101.68,114.31,100.59,104.23,107.60,109.24,108.14,114.07,108.23,108.72,103.41,98.53,97.92,106.26,103.75,111.58,99.80,107.49,106.80,104.97,102.73,111.66,105.91,98.13,110.04,106.92,96.52,112.76,99.59,102.97,103.59,109.58,112.60,103.14,106.84,106.41,97.35,104.41,104.63,101.51,105.77,105.67,98.19,108.83,113.11,101.35,111.37,104.29,107.72,94.20,106.92,110.87,101.66,106.85,108.92,106.38,100.23,104.60,103.08,104.64,116.32,113.37,96.69,100.38,103.09,97.05,111.36,98.98,107.26,108.81,103.00,104.20,103.46,106.58,110.86,100.72,102.94,111.61,108.08,103.89,116.51,107.30,107.00,105.31,105.68,105.98,112.01,110.07,111.89,106.79,102.70,111.20,110.53,110.97,107.91,107.33,112.78,108.05,105.53,108.56,107.78,106.71,109.18,109.05,106.46,100.84,109.23,115.56,107.60,106.54,109.31,97.78,111.07,109.14,104.85,108.65,102.94,105.63,106.89,104.14,105.76,107.11,101.60,104.11,109.49,110.74,100.19,102.07,106.04,107.45,102.02,105.27,106.17,107.37,107.92,104.53,106.22,110.01,105.74,113.74,109.10,107.88,108.08,107.22,108.23,109.94,96.04,110.51,104.26,107.76,105.06,94.03,106.14,102.56,108.55,105.85,95.46,107.24,109.22,105.57,102.95,103.30))
# initial values for simulation
list(beta1=c(0,0,0,0), beta2=c(0,0,0,0),beta3=c(0,0,0,0),
precy1=c(1,1), precy2=c(1,1), precy3=c(1,1))
BUGS CODE 6.3
model
{
g[1,1]~dgamma(12,2)
g[1,2]~dgamma(3,2)
g[1,3]~dgamma(20,2)
g[2,1] ~dgamma(2,2)
g[2,2] ~dgamma(4,2)
g[2,3]~dgamma(4,2)
g[3,1]~dgamma(1,2)
g[3,2]~dgamma(9,2)
g[3,3]~dgamma(83,2)
h<- sum(g[,])
for( i in 1 : 3 ) {for( j in 1 :3 ){ theta[i,j]<-g[i,j]/h }}
theta1.<-sum(theta[1,])
theta.1<-sum(theta[,1])
theta2.<-sum(theta[2,])
theta.2<-sum(theta[,2])
theta3.<-sum(theta[3,])
theta.3<- sum(theta[,3])
kappa<- (agree-cagree)/(1-cagree)
agree<-theta[1,1]+theta[2,2]+theta[3,3]
cagree<-theta1.*theta.1+theta2.*theta.2+theta3.*theta.3
d<-agree-cagree
Kappa1<-(theta[1,1]-theta1.*theta.1)/(theta1.-theta1.*theta.1)
}
list( g = structure(.Data = c(2,2,2,2,2,2,2,2,2), .Dim = c(3,3)))
BUGS CODE 6.4
Model;
{
# generates Dirichlet Distribution
g[1,1]~dgamma(12,2); g[1,2]~dgamma(3,2); g[1,3]~dgamma(20,2);
g[2,1] ~dgamma(2,2); g[2,2] ~dgamma(4,2); g[2,3]~dgamma(4,2)
g[3,1]~dgamma(1,2); g[3,2]~dgamma(9,2); g[3,3]~dgamma(83,2);
h<- sum(g[,]) ;
for( i in 1 : 3 ) {for( j in 1 :3 ){ theta[i,j]<-g[i,j]/h }}
theta1.<-sum(theta[1,]); theta.1<-sum(theta[,1]);
theta2.<-sum(theta[2,]); theta.2<-sum(theta[,2]);
theta3.<-sum(theta[3,]); theta.3<- sum(theta[,3]);
# kappa is weighted kappa
kappa<- (agree-cagree)/(1-cagree);
agree<- w[1,1]*theta[1,1]+w[1,2]*theta[1,2]+w[1,3]*theta[1,3]+w[2,1]*theta[2,1]+w[2,2]*theta[2,2]+w[2,3]*theta[2,3] +w[3,1]*theta[3,1]+w[3,2]*theta[3,2] +w[3,3]*theta[3,3];
cagree<-w[1,1]*theta1.*theta.1+w[1,2]*theta1.*theta.2+w[1,3]*theta1.*theta.3+w[2,1]*theta2.*theta.1+w[2,2]*theta2.*theta.2+w[2,3]*theta2.*theta.3 +w[3,1]*theta3.*theta.1+w[3,2]*theta3.*theta.2 +w[3,3]*theta3.*theta.3;
d<-agree-cagree ;
}
# inital values
list(g = structure(.Data = c(2,2,2,2,2,2,2,2,2), .Dim = c(3,3)))
# data
list(w = structure(.Data = c(1,.5,0,.5,1,.5,0,.5,1), .Dim = c(3,3)))
BUGS CODE 6.5
model
{
phi~dbeta(1,1)
for( i in 1 : 5) {
for( j in 1 : 5 ) {
Y[i , j] ~ dbern(phi)}}
theta11~dbeta(20,30)
kappa<-(theta11-phi*phi)/phi *(1-phi)
}
list(Y= structure( .Data=c(1,1,1,1,1,
0,1,0,1,0,1,1,0,0,1,1,1,1,0,0,
0,0,1,1,1), .Dim=c(5,5)))
list(theta11=.5 )
BUGS CODE 6.6
Model;
{
g00~dgamma(a00,b00); g01~dgamma(a01,b01);
g10~dgamma(a10,b10); g11 ~dgamma(a11,b11);
h<- g00+g01+ g10 +g11; theta00<-g00/h;
theta01<-g01/h; theta10<-g10/h; theta11<-g11/h;
theta0.<-theta00+theta01; theta.0<-theta00+theta10;
theta1.<-theta10+theta11; theta.1<-theta01+theta11;
kappa<- (agree-cagree)/(1-cagree); g<-agree-cagree;
G<- (1/2)*(theta11/(theta11+theta01) +theta11/(theta11+theta10));
J<-theta11/(theta11+theta10+theta01); agree <-theta11+theta00;
cagree<-theta1.*theta.1+theta0.*theta.0;
}
list( a00 = 22, b00 = 2, a01 = 15, b01=2, a10= 33 , b10=2, a11=30 , b11=2)
list( g00 = 2, g01 =2, g10 =2, g11=2)
BUGS CODE 6.7
model;
{
g00~dgamma(a00,b00); g01~dgamma(a01,b01); g10~dgamma(a10,b10); g11 ~dgamma(a11,b11);
h<- g00+g01+ g10 +g11 ; theta00<-g00/h;
theta01<-g01/h; theta10<-g10/h; theta11<-g11/h;
theta0.<-theta00+theta01; theta.0<-theta00+theta10;
theta1.<-theta10+theta11; theta.1<-theta01+theta11;
p00~dgamma(c00,d00); p01~dgamma(c01,d01); p10~dgamma(c10,d10); p11 ~dgamma(c11,d11);
q<-p00+p01+p10+p11; phi00<-p00/q;
phi01<-p01/q; phi10<-p10/q; phi11<-p11/q;
phi0.<-phi00+phi01; phi.0<-phi00+phi10;
phi1.<-phi10+phi11; phi.1<-phi01+phi11 ;
kappatw<-k11/k12; k11<-2*(theta.1)*(1-theta.1)*(1-theta01/theta.1-theta10/theta.0)*(1-phi01/phi.1-phi10/phi.0);
k12<-theta1.*(1-theta1.)+ phi1.*(1-phi1.);
}
# data
list( a00 =85, b00 = 2, a01 = 25, b01=2, a10= 15, b10=2, a11=110 , b11=2,
c00 =2, d00 = 2,c01 = 85, d01=2, c10= 98, d10=2, c11=50 , d11=2)
# initial values
list( g00 = 2, g01 =2, g10 =2, g11=2, p00 = 2, p01 =2, p10 =2, p11=2)
BUGS CODE 6.8
{
# the following assumes a uniform prior
# Uses information in Table 6.30
# The 16 observations are labeled 1 to 16 from top to bottom
g[1]~dgamma(4,2); g[2]~dgamma(1,2);
g[3]~dgamma(3,2); g[4] ~dgamma(1,2);
g[5] ~dgamma(4,2); g[6]~dgamma(1,2);
g[7]~dgamma(1,2); g[8]~dgamma(1,2);
g[9]~dgamma(1,2); g[10]~dgamma(2,2);
g[11]~dgamma(2,2); g[12]~dgamma(1,2);
g[13]~dgamma(3,2); g[14]~dgamma(2,2);
g[15]~dgamma(1,2); g[16]~dgamma(8,2);
h<- sum(g[])
for( i in 1:16){ theta[i]<- g[i]/h}
# The following are the terms for chance agreement
atheta2...<-theta[5]+theta[8]+theta[10]+theta[11]+theta[12]+theta[13]+theta[14]+theta[16];
atheta.2..<-
theta[4]+theta[7]+theta[9]+theta[11]+theta[12]+theta[13]+theta[15]+theta[16];
atheta..2.<- theta[3]+theta[6]+theta[9]+theta[10]+theta[12]+theta[14]+theta[15]+theta[16];
atheta...2<-
theta[2]+theta[6]+theta[7]+theta[8]+theta[13]+theta[14]+theta[15]+theta[16];
atheta1...<-
theta[1]+theta[2]+theta[3]+theta[4]+theta[6]+theta[7]+theta[9]+theta[15];
atheta.1..<-
theta[1]+theta[2]+theta[3]+theta[5]+theta[6]+theta[8]+theta[10]+theta[14];
atheta..1.<-
theta[1]+theta[2]+theta[4]+theta[5]+theta[7]+theta[8]+theta[11]+theta[13];
atheta...1<-
theta[1]+theta[3]+theta[4]+theta[5]+theta[9]+theta[10]+theta[11]+theta[12];
kappa42<- (agree-cagree)/(1-cagree)
agree<-theta[1]+theta[16]
cagree<-atheta1...*atheta.1..*atheta..1.*atheta...1+ atheta2...*atheta.2..*atheta..2.*atheta...2
d<- agree-cagree
}
# initail values
list(g = c(2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2))
BUGS CODE 6.9
model
{
# g1 corresponds to the first case
g1~dgamma(15,2) g3~dgamma(2,2)
g5~dgamma(2,2) g6 ~dgamma(1,2)
g10~dgamma(1,2)
g14~dgamma(1,2) g16~dgamma(1,2)
g18~dgamma(1,2) g22~dgamma(1,2)
g25~dgamma(1,2) g26~dgamma(1,2)
g27~dgamma(1,2) g35~dgamma(3,2)
g36~dgamma(2,2) g37~dgamma(2,2)
g41~dgamma(2,2) g42~dgamma(1,2)
g43~dgamma(30,2)
h<- g1+ g3+ g5+ g6 + g10+ g14+ g16+ g18+ g22+ g25+g26+g27+g35+g36
+g37+ g41+g42+g43
# a1 is the probability of the outcome for case 1
a1<- g1/h; a3<- g3/h; a5<- g5/h; a6<- g6/h; a10<- g10/h;
a14<- g14/h ; a16<- g16/h; a18<- g18/h; a22<- g22/h;
a25<- g25/h ; a26<- g26/h; a27<- g27/h; a35<- g35/h;
a36<- g36/h; a37<- g37/h; a41<- g41/h; a42<- g42/h;
a43<- g43/h ;
agree<- a1+a43
GC<-agree -(1-agree)
# a1….. is the probability that pathologist A assigns a 1 to a lesion
a1.....<- a1+a3+a5+a6+a10+a14+a16+a18
a.1....<- a1+a3+a5+a6+ a10+a22+a25+a26+a27
a..1...<- a1+a3+a5+a14+a22+a25+a26+a35+a36+a37
a...1..<- a1+a3+a14+a18+a22
a....1.<- a1+a5+a10+a18+a22+a25+a35+a41
a.....1<- a1+a3+a5+a6+a10+a22+a25+a26+a35+a36+a42
a2.....<- a22+a25+a26+a27+a35+a36+a37+a41+a42+a43
a.2....<- a14+a16+a18+a35+a36+a37+a41+a42+a43
a..2...<- a6+a10+a16+a18+a27+a41+a42+a43
a...2..<- a5+a6 +a10+a16+a25+a26+a27+a35+a36+a37+a41+a42+a43
a....2.<- a3+a6+a14+a16+a26+a27+a36+a37+a42+a43
a.....2<- a14+a16+a18+a27+a37+a41+a43
# a…..2 is the probability that reader F assigns a score of 2 to a lesion
cagree<- (a1.....*a.1....*a..1...*a...1..*a....1.*a.....1)+(a2.....*a.2....*a..2...*a...2..*a....2.*a.....2)
# kappa is the index that all six agree
kappa<- (agree-cagree)/(1-cagree)
al5<- a1+a3+a5+a22+a37+a41+a42+a43
#al5 is the probability at least five agree
cal5<- a1.....*a.1....*a..1...*a...1..*a....1.*a.....1+ a1.....*a.1....*a..1...*a...1..*a....2.*a.....1 + a1.....*a.1....*a..1...*a...2..*a....1.*a.....1+a2.....*a.1....*a..1...*a...1..*a....1.*a.....1+a2.....*a.2....*a..1...*a...2..*a....2.*a.....2+a2.....*a.2....*a..2...*a...2..*a....1.*a.....2+
a2.....*a.2....*a..2...*a...2..*a....2.*a.....1+a2.....*a.2....*a..2...*a...2..*a....2.*a.....2
# cal5 is the probability at least 5 agree by chance
# kappa5 is the index that at least 5 agree
kappa5<-(al5-cal5)/(1-cal5)
GC5<- al5 -(1-al5)
}
list( g1=2, g3=2, g5 =2, g6=2, g14=2, g16=2, g18=2, g22=2, g25=2, g26=2, g27 =2, g35=2, g36 =2,g37 =2, g41=2, g42=2,g43=2).
BUGS CODE 6.10
model
{
# for site 3
# g[i] corresponds to case i of Table 6.35c.
# a uniform prior is assumed
g[1]~dgamma(101,2) g[2]~dgamma(41,2) g[3]~dgamma(23,2)
g[4]~dgamma(24,2) g[5]~dgamma(5,2) g[6]~dgamma(11,2)
g[7]~dgamma(4,2) g[8]~dgamma(4,2) g[9]~dgamma(4,2)
g[10]~dgamma(3,2) g[11]~dgamma(6,2) g[12]~dgamma(3,2)
g[13]~dgamma(2,2) g[14]~dgamma(2,2) g[15]~dgamma(.1,2)
g[16]~dgamma(3,2) g[17]~dgamma(11,2) g[18]~dgamma(3,2)
g[19]~dgamma(3,2) g[20]~dgamma(2,2) g[21]~dgamma(21,2)
g[22]~dgamma(61,2) g[23]~dgamma(12,2) g[24]~dgamma(11,2)
g[25]~dgamma(4,2) g[26]~dgamma(3,2) g[27]~dgamma(11,2)
g[28]~dgamma(2,2) g[29]~dgamma(1,2) g[30]~dgamma(11,2)
g[31]~dgamma(2,2) g[32]~dgamma(2,2) g[33]~dgamma(5,2)
g[34]~dgamma(3,2) g[35]~dgamma(2,2) g[36]~dgamma(2,2)
g[37]~dgamma(3,2) g[38]~dgamma(11,2) g[39]~dgamma(5,2)
g[40]~dgamma(5,2) g[41]~dgamma(6,2) g[42]~dgamma(6,2)
g[43]~dgamma(61,2) g[44]~dgamma(11,2) g[45]~dgamma(3,2)
g[46]~dgamma(2,2) g[47]~dgamma(2,2) g[48]~dgamma(4,2)
g[49]~dgamma(5,2) g[50]~dgamma(2,2) g[51]~dgamma(2,2)
g[52]~dgamma(3,2) g[53]~dgamma(.1,2) g[54]~dgamma(2,2)
g[55]~dgamma(11,2) g[56]~dgamma(1,2) g[57]~dgamma(1,2)
g[58]~dgamma(.01,2) g[59]~dgamma(1,2) g[60]~dgamma(.01,2)
g[61]~dgamma(3,2) g[62]~dgamma(2,2) g[63]~dgamma(2,2)
g[64]~dgamma(8,2)
h<- sum(g[])
for(i in 1:64){a[i]<- g[i]/h}
agree3<- a[1]+a[21]+a[43]+a[64]
b1..<-
a[1]+a[2]+a[3]+a[4]+a[5]+a[6]+a[7]+a[8]+a[9]+a[10]+a[11]+a[12]+a[13]+
a[14]+a[15]+a[16]
b.1.<-
a[1]+a[2]+a[3]+a[4]+a[17]+a[18]+a[19]+a[20]+a[33]+a[34]+a[35]+a[36]+a[49]+a[50]+a[51]+a[52]
b..1 <- a[1]+a[5]+a[9]+a[13]+a[17]+a[21]+a[25]+a[29]+a[33]+a[37]+a[41]+a[45]+a[49]+a[53]+a[57]+a[61]
b2..<- a[17]+a[18]+a[19]+a[20]+a[21]+a[22]+a[23]+a[24]+a[25]+a[26]+a[27]+a[28]+a[29]+a[30]+a[31]+a[32]
b.2.<-
a[5]+a[6]+a[7]+a[8]+a[21]+a[22]+a[23]+a[24]+a[37]+a[38]+a[39]+a[40]+a[53]+a[54]+a[55]+a[56]
b..2<- a[2]+a[6]+a[10]+a[14]+a[18]+a[22]+a[26]+a[30]+a[34]+a[38]+a[42]+a[46]+a[50]+a[54]+a[58]+a[62]
b3..<- a[33]+a[34]+a[35]+a[36]+a[37]+a[38]+a[39]+a[40]+a[41]+a[42]+a[43]+a[44]+a[45]+a[46]+a[47]+a[48]
b.3.<- a[9]+a[10]+a[11]+a[12]+a[25]+a[26]+a[27]+a[28]+a[41]+a[42]+a[43]+a[44]+a[57]+a[58]+a[59]+a[60]
b..3<- a[3]+a[7]+a[11]+a[15]+a[19]+a[23]+a[27]+a[31]+a[35]+a[39]+a[43]+a[47]+a[51]+a[55]+a[59]+a[63]
b4..<- a[49]+a[50]+a[51]+a[52]+a[53]+a[54]+a[55]+a[56]+a[57]+a[58]+a[59]+a[60]+a[61]+a[62]+a[63]+a[64]
b.4.<- a[13]+a[14]+a[15]+a[16]+a[29]+a[30]+a[31]+a[32]+a[45]+a[46]+a[47]+a[48]+a[61]+a[62]+a[63]+a[64]
b..4<- a[4]+a[8]+a[12]+a[16]+a[20]+a[24]+a[28]+a[32]+a[36]+a[40]+a[44]+a[48]+a[52]+a[56]+a[60]+a[64]
cagree3<- b1..*b.1.*b..1+b2..*b.2.*b..2+b3..*b.3.*b..3+b4..*b.4.*b..4
kappa3<- (agree3-cagree3)/(1-cagree3) }
list( g=c(2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2))
BUGS CODE 6.11
model
{ for (i in 1:17) {phi[i]~dbeta(1,1)}
for( i in 1 : 17) {
for( j in 1 : 10) {
Y[i , j] ~ dbern(phi[i])}}
theta11~dbeta(17,213)
kappa<-(theta11-pi*pi)/pi*(1-pi)
pi<-mean(phi[]) }
list(Y= structure( .Data=
c(1,0,0,0,0,0,0,0,0,NA,0,0,0,0,0,0,0,0,0,0,
1,0,0,0,0,0,0,NA,NA,NA,0,0,0,0,0,NA,NA,NA,NA, NA,
1,0,0,0,NA,NA,NA,NA, NA,NA,0,0,0,0,0,0,NA,NA, NA,NA,
1,0,0,NA,NA,NA,NA, NA,NA,NA,1,0,0,0,0,0,0,0,NA,NA,
1,1,0,0,0,0,0,0,NA,NA,0,0,0,0,NA,NA,NA, NA,NA,NA,
1,1,1,1,NA,NA,NA, NA,NA,NA,1,0,0,0,0,NA,NA, NA,NA,NA,
1,0,0,NA,NA,NA, NA,NA,NA,NA,1,1,1,1,0,0,0,0,NA,NA,
1,1,0,0,0,0,NA,NA,NA,NA,1,1,1,0,0,0,0,0,NA,NA,
1,0,0,0,0,0,NA,NA,NA,NA), .Dim=c(17,10)))
list(theta11=.5 )
No comments:
Post a Comment