BUGS CODE 7.1
model;
{
# conditional distribution of p
p~dbeta(alp,bep)
alp<- y11+y10+y01+y00+ap
bep<- sn -ap+bp
sn<- n11+n10+n01+n00
# conditional distribution of s1
s1~dbeta(als1,bes1)
als1<-y11+y01+as1
bes1<-y10+y00+bs1
#conditional distribution of s2
s2~dbeta(als2,bes2)
als2<-y11+y10+as2
bes2<-y01+y00+bs2
# conditional distribution of c1
c1~dbeta(alc1,bec1)
alc1<-n10+n00-y10-y00+ac1
bec1<-n11+n01-y11-y01+bc1
# conditional distribution of c2
c2~dbeta(alc2,bec2)
alc2<- n01+n00-y01-y00+ac2
bec2<- n11+n10-y11-y10+bc2
# conditional distribution of y11
y11~dbin( py11,n11)
py11<-p*s1*s2/(p*s1*s2+(1-c1)*(1-c2)*(1-p))
# conditional distribution of y10
y10~dbin(py10,n10)
py10<-p*(1-s1)*s2/(p*(1-s1)*s2+c1*(1-c2)*(1-p))
# conditional distribution of y01
y01~dbin(py01,n01)
py01<-p*s1*(1-s2)/(p*s1*(1-s2)+c2*(1-c1)*(1-p))
# conditional distribution of y00
y00~dbin(py00,n00)
py00<-p*(1-s1)*(1-s2)/(p*(1-s1)*(1-s2)+c1*c2*(1-p))
}
# hyperparameters for uniform prior Joseph
list(ap=1,bp=1,as1=1,bs1=1,as2=1,bs2=1,
ac1=1,bc1=1,ac2=1,bc2=1, n11=38,n10=87,n01=2,n00=35)
# hyperparameters for informative prior joseph
list(ap=80,bp=80,as1=4.44,bs1=13.31,as2=21.96,
bs2=5.49,ac1=71.25,bc1=3.75,ac2=4.1,bc2=1.76, n11=38,n10=87,n01=2,n00=35)
# initial values
list(p=.5,c1=.95,c2=.9,s1=.9,s2=.9,y11=1,y10=1,
y01=1,y00=00)
BUGS CODE 7.2
# accuracy parameters for two binary tests R and T
# the sensitivity and specificity of R and T
# T is the new and R is the reference
model;
{
p~dbeta(y..,m)
# p and d are the disease rate.
d<- p
m<-n-y..+1
y..<-y11+y10+y01+y00+ 1
y11~dbin(a11,b11)
a11<- theta11*d/(theta11*d+ph11*(1-d))
y10~dbin(a10,b10)
a10<-theta10*d/(theta10*d+ph10*(1-d))
y01~dbin(a01,b01)
a01<- theta01*d/(theta01*d+ph01*(1-d))
y00~dbin(a00,b00)
a00<-theta00*d/(theta00*d+ph00*(1-d))
# Dirichlet distribution for the thetas
g11~dgamma(r11,2)
g10~dgamma(r10,2)
g01~dgamma(r01,2)
g00~dgamma(r00,2)
r11<-y11+1
r10<-y10+1
r01<-y01+1
r00<-y00+1
sg<-g11+g10+g01+g00
theta11<-g11/sg
theta10<-g10/sg
theta01<-g01/sg
theta00<-g00/sg
# Dirichlet distribution for the phs
h11~dgamma(s11,2)
h10~dgamma(s10,2)
h01~dgamma(s01,2)
h00~dgamma(s00,2)
s11<-z11+1
s10<-z10+1
s01<-z01+1
s00<-z00+1
sh<-h11+h10+h01+h00
ph11<-h11/sh
ph10<-h10/sh
ph01<-h01/sh
ph00<-h00/sh
z11<-b11-y11
z10<-b10-y10
z01<-b01 -y01
z00<-b00-y00
# tests for conditional independence of the thetas
theta11c<- (theta11+theta10)*(theta11+theta01)-theta11
theta10c<- (theta11+theta10)*(theta10+theta00)-theta10
theta01c<- (theta11+theta01)*(theta01+theta00)-theta01
theta00c<- (theta10+theta00)*(theta01+theta00)-theta00
# tests for conditional independence for phs
ph11c<- (ph11+ph10)*(ph11+ph01)-ph11
ph10c<- (ph11+ph10)*(ph10+ph00)-ph10
ph01c<- (ph11+ph01)*(ph01+ph00)-ph01
ph00c<- (ph10+ph00)*(ph01+ph00)-ph00
#sensitivity and specificity
# sensitivity and specificity of the stool exam
s1<- (theta11+ theta01)
c1<- ph10+ph00
# sensitivity and specificity of the serology exam
s2<-theta11+theta10
c2<- ph01+ph00
}
# Strongyloides example taken from Joseph
list(n=163, b11=38,b10=87,b01=2,b00=35)
# starting values
list( g11=1,g10=1,g01=1, g00=1, p=.5, h11=1,h10=1,h01=1,h00=1, y11=1,y10=1,y01=1,
y00=1))
BUGS CODE 7.3
model;
{
# conditional distribution of p
#p0 is an informative prior for p
# note that the distribution for p0 has be deactivated by the comment sign #!
# p0~dbeta(ap0,bp0)
# the following are hyperparamerers for the disease rate of the second site
ap0<- 990
bp0<-10
# p1 depends on the data only
p1~dbeta(alp,bep)
# p can be expressed as a mixture of p0 and p1
p<- p1
alp<- y11+y10+y01+y00+ap
bep<- sn -alp+bp
sn<- n11+n10+n01+n00+1
# conditional distribution of s1
s1~dbeta(als1,bes1)
als1<-y11+y01+as1
bes1<-y10+y00+bs1
#conditional distribution of s2
s2~dbeta(als2,bes2)
als2<-y11+y10+as2
bes2<-y01+y00+bs2
# conditional distribution of c1
c1~dbeta(alc1,bec1)
alc1<-n10+n00-y10-y00+ac1
bec1<-n11+n01-y11-y01+bc1
# conditional distribution of c2
c2~dbeta(alc2,bec2)
alc2<- n01+n00-y01-y00+ac2
bec2<- n11+n10-y11-y10+bc2
# conditional distribution of y11
y11~dbin( py11,n11)
py11<-p*s1*s2/(p*s1*s2+(1-c1)*(1-c2)*(1-p))
# conditional distribution of y10
y10~dbin(py10,n10)
py10<-p*(1-s1)*s2/(p*(1-s1)*s2+c1*(1-c2)*(1-p))
# conditional distribution of y01
y01~dbin(py01,n01)
py01<-p*s1*(1-s2)/(p*s1*(1-s2)+c2*(1-c1)*(1-p))
# conditional distribution of y00
y00~dbin(py00,n00)
py00<-p*(1-s1)*(1-s2)/(p*(1-s1)*(1-s2)+c1*c2*(1-p))
}
# hyperparameters for Tuberculosis example
# Page 371 of Zhou site 1
# Test 1 is Tine Test
# Test 2 is Mantour test
# uniform prior for accuracy parameters
list(ap= 1, bp=1, as1=1,bs1=1,as2=1,bs2=1,
ac1=1,bc1=1,ac2=1,bc2=1, n11=14,n10=4,n01=9,n00=528)
# hyper parameters for Tuberculosis example
# Page 371 of Zhou site 2
# Test 1 is Tine Test
# Test 2 is Mantour test
# uniform prior for accuracy parameters
list(ap=1,bp=1,as1=1,bs1=1,as2=1,bs2=1,
ac1=1,bc1=1,ac2=1,bc2=1, n11=887,n10=31,n01=37,n00=367)
# initial values
list(c1=.5,c2=.5,s1=.5,s2=.5,y11=1,y10=1,
y01=1,y00=1)
BUGS CODE 7.4
model;
{
p~dbeta(y..,m)
# note the mixture is not used for d
p0~dbeta(ap0,bp0)
# d is a mixture of the uniform prior and an
# informative prior p0
# d is the posterior distribution of disease rate
d<- .10*p+.90*p0
m<-n-y..+1
y..<-y11+y10+y01+y00+ 1
y11~dbin(a11,b11)
a11<- theta11*d/(theta11*d+ph11*(1-d))
y10~dbin(a10,b10)
a10<-theta10*d/(theta10*d+ph10*(1-d))
y01~dbin(a01,b01)
a01<- theta01*d/(theta01*d+ph01*(1-d))
y00~dbin(a00,b00)
a00<-theta00*d/(theta00*d+ph00*(1-d))
g11~dgamma(r11,2)
g10~dgamma(r10,2)
g01~dgamma(r01,2)
g00~dgamma(r00,2)
r11<-y11+1
r10<-y10+1
r01<-y01+1
r00<-y00+1
sg<-g11+g10+g01+g00
theta11<-g11/sg
theta10<-g10/sg
theta01<-g01/sg
theta00<-g00/sg
h11~dgamma(s11,2)
h10~dgamma(s10,2)
h01~dgamma(s01,2)
h00~dgamma(s00,2)
s11<-z11+1
s10<-z10+1
s01<-z01+1
s00<-z00+1
sh<-h11+h10+h01+h00
ph11<-h11/sh
ph10<-h10/sh
ph01<-h01/sh
ph00<-h00/sh
z11<-b11-y11
z10<-b10-y10
z01<-b01 -y01
z00<-b00-y00
# test for conditional independence
theta11c<- (theta11+theta10)*(theta11+theta01)-theta11
theta10c<- (theta11+theta10)*(theta10+theta00)-theta10
theta01c<- (theta11+theta01)*(theta01+theta00)-theta01
theta00c<- (theta10+theta00)*(theta01+theta00)-theta00
# test for conditional independence for phs
ph11c<- (ph11+ph10)*(ph11+ph01)-ph11
ph10c<- (ph11+ph10)*(ph10+ph00)-ph10
ph01c<- (ph11+ph01)*(ph01+ph00)-ph01
ph00c<- (ph10+ph00)*(ph01+ph00)-ph00
#sensitivity and specificity
s1<- (theta11+ theta01)
c1<- ph10+ph00
s2<-theta11+theta10
c2<- ph01+ph00
}
# Strongyloides example taken from Joseph
list(n=163, b11=38,b10=87,b01=2,b00=35)
# Tine test & Mantour test site 1
list(n=556, ap0=30, bp0=970,b11=14,b10=4,b01=9,b00=528)
# Tine test & Mantour test site 2
list( n=1323, ap0=980, bp0=20, b11=887,b10=31,b01=37,b00=367)
# starting values
list( g11=1,g10=1,g01=1, g00=1, p=.5, h11=1,h10=1,h01=1,h00=1, y11=1,y10=1,y01=1,
y00=1))
BUGS CODE 7.5
model;
{
p1~dbeta(alphap1,betap1)
p0~dbeta(100,900)
# note p is not expressed as a mixture for this example
p<-p1
alphap1<- sy+1
betap1<- sn -alphap1
sn<-n111+n101+n110+n000+n011+n001+n010+n000
sy<-y111+y101+y110+y000+y011+y001+y010+y000
s1~dbeta(alphas1,betas1)
alphas1<-y111+y101+y110+y100+as1
betas1<- y011+y001+y010+y000+ bs1
s2~dbeta(alphas2,betas2)
alphas2<-y111+y110+y011+y010+as2
betas2<- y101+y100+y001+y000+bs2
s3~dbeta(alphas3,betas3)
alphas3<- y111+y101+y011+y001+as3
betas3<- y110+y100+y010+y000+bs3
c1~dbeta(alphac1,betac1)
alphac1<-n011+n001+n010+n000-y011-y001
-y010-y000+ac1
betac1<- n111+n101+n110+n100-y111-y101
-y110-y100+bc1
c2~dbeta(alphac2,betac2)
alphac2<-n101+n100+n001+n000-y101-y100-
y001-y000 +ac2
betac2<- n111+n110+n011+n010-y111-y110
-y011-y010+bc2
c3~dbeta(alphac2,betac3)
alphac3<-n110+n100+n010+n000-y110+y100+y010+y000+ac3
betac3<- n111+n101+n011+n001-y111+y101+y011+y001+bc3
y111~dbin( py111,n111)
py111<-p*s1*s2*s3/(p*s1*s2*s3+(1-c1)*(1-c2)*(1-c3)*(1-p))
y101~dbin(py101,n101)
py101<-p*s1*(1-s2)*s3/(p*s1*(1-s2)*s3+(1-c1)*c2*(1-c3*(1-p)))
y110~dbin(py110,n110)
py110<-p*s1*s2*(1-s3)/(p*s1*s2*(1-s3)+c3*(1-c2)*(1-c1)*(1-p))
y100~dbin(py100,n100)
py100<-p*s1*(1-s2)*(1-s3)/(p*s1*(1-s2)*(1-s3)+(1-c1)*c2*c3*(1-p))
y011~dbin(py011,n011)
py011<- p*(1-s1)*s2*s3/(p*(1-s1)*s2*s3+(1-p)*c1*(1-c2)*(1-c3))
y001~dbin(py001,n001)
py001<- p*(1-s1)*(1-s2)*s3/(p*(1-s1)*(1-s2)*s3+(1-p)*c1*c2*(1-c3))
y010~dbin(py010,n010)
py010<- p*(1-s1)*s2*(1-s3)/(p*(1-s1)*s2*(1-s3)+(1-p)*c1*(1-c2)*c3)
y000~dbin(py000,n000)
py000<- p*(1-s1)*(1-s2)*(1-s3)/(p*(1-s1)*(1-s2)*(1-s3)+(1-p)*c1*c2*c3)
}
# hyperparameters for pepe example
# Test1 is culture, 2 is Elisa, 3 is PCR
list(n111=20,n101=4,n110=2,n100=2,
n011=4,n001=3,n010=2,n000=292,
as1=1,bs1=1,as2=1,bs2=1,
ac1=1,bc1=1,ac2=1,bc2=1,
as3=1,bs3=1,ac3=1,bc3=1)
# hyperparameters for prostae enlargement
list(n111=46,n101=27,n110=20,n100=31,n011=31,
n001=41,n010=79,n000=1533,
as1=1,bs1=1,as2=1,bs2=1,
ac1=1,bc1=1,ac2=1,bc2=1,
as3=1,bs3=1,ac3=1,bc3=1)
# hyperparameters Pepe table 7.14
# Page 204
list(n111=70,n101=25,n110=5,n100=10,
n011=110,n001=150,n010=100,n000=530,
as1=1,bs1=1,as2=1,bs2=1,
ac1=1,bc1=1,ac2=1,bc2=1,
as3=1,bs3=1,ac3=1,bc3=1)
# initial values
list(c1=.5,c2=.5,s1=.5,s2=.5,s3=.5,c3=.5,
y111=1,y101=1,y110=1,y100=1,y011=1,y001=1,
y010=1,y000=1)
BUGS CODE 7.6
model;
# three tests wo cia
{
p~dbeta(y...,m)
# p0~dbeta(ap0,bp0)
# d is a mixture of the uniform prior and an
# informative prior p0
# if d is used as a mixture ap0 and bp0 must be defined in the list statement for the data
# d is the posterior distribution of disease rate
d<- p
m<-n...-y...+1
y...<- y111+y101+y110+y100+y011+y001+y010+y000+1
n...<- n111+n101+n110+n100+n011+n001+n010+n000
y111~dbin(a111,n111)
a111<- theta111*d/(theta111*d+ph111*(1-d))
y110~dbin(a110,n110)
a110<-theta110*d/(theta110*d+ph110*(1-d))
y101~dbin(a101,n101)
a101<- theta101*d/(theta101*d+ph101*(1-d))
y100~dbin(a100,n100)
a100<-theta100*d/(theta100*d+ph100*(1-d))
y011~dbin(a011,n011)
a011<-theta011*d/(theta011*d+ph011*(1-d))
y001~dbin(a001,n001)
a001<-theta001*d/(theta001*d+ph001*(1-d))
y010~dbin(a010,n010)
a010<-theta010*d/(theta010*d+ph010*(1-d))
y000~dbin(a000,n000)
a000<-theta000*d/(theta000*d+ph000*(1-d))
g111~dgamma(r111,2)
g110~dgamma(r110,2)
g101~dgamma(r101,2)
g100~dgamma(r100,2)
g011~dgamma(r011,2)
g010~dgamma(r010,2)
g001~dgamma(r001,2)
g000~dgamma(r000,2)
r111<-y111+1
r110<-y110+1
r101<-y101+1
r100<-y100+1
r011<-y011+1
r010<-y010+1
r001<-y001+1
r000<-y000+1
sg<-g111+g110+g101+g100+g011+g010+g001+g000
theta111<-g111/sg
theta110<-g110/sg
theta101<-g101/sg
theta100<-g100/sg
theta011<-g011/sg
theta010<-g010/sg
theta001<-g001/sg
theta000<-g000/sg
h111~dgamma(s111,2)
h110~dgamma(s110,2)
h101~dgamma(s101,2)
h100~dgamma(s100,2)
h011~dgamma(s011,2)
h010~dgamma(s010,2)
h001~dgamma(s001,2)
h000~dgamma(s000,2)
s111<-z111+1
s110<-z110+1
s101<-z101+1
s100<-z100+1
s011<-z011+1
s010<-z010+1
s001<-z001+1
s000<-z000+1
sh<-h111+h110+h101+h100 + h011+h010+h001+h000
ph111<-h111/sh
ph110<-h110/sh
ph101<-h101/sh
ph100<-h100/sh
ph011<-h011/sh
ph010<-h010/sh
ph001<-h001/sh
ph000<-h000/sh
z111<-n111-y111
z110<-n110-y110
z101<-n101 -y101
z100<-n100-y100
z011<-n011-y011
z010<-n010-y010
z001<-n001 -y001
z000<-n000-y000
# test for conditional independence
theta111c<- (theta111+theta110)*(theta111+theta101)-
theta111
theta110c<- (theta111+theta110)*(theta110+theta100)-
theta110
theta101c<- (theta111+theta101)*(theta101+theta100)-
theta101
theta100c<- (theta110+theta100)*(theta101+theta100)-
theta100
theta011c<- (theta011+theta010)*(theta011+theta001)-
theta011
theta010c<- (theta011+theta010)*(theta010+theta000)-
theta010
theta001c<- (theta011+theta001)*(theta001+theta000)-
theta001
theta000c<- (theta010+theta000)*(theta001+theta100)-
theta000
# test for conditional independence for phs
ph111c<- (ph111+ph110)*(ph111+ph101)-ph111
ph110c<- (ph111+ph110)*(ph110+ph100)-ph110
ph101c<- (ph111+ph101)*(ph101+ph100)-ph101
ph100c<- (ph110+ph100)*(ph101+ph100)-ph100
ph011c<- (ph011+ph010)*(ph011+ph001)-ph011
ph010c<- (ph011+ph010)*(ph010+ph000)-ph010
ph001c<- (ph011+ph001)*(ph001+ph000)-ph001
ph000c<- (ph010+ph000)*(ph001+ph000)-ph000
#sensitivity and specificity
s1<- theta111+ theta101+theta110+theta100
c1<- ph011+ph001+ph010+ph000
s2<-theta111+theta110+theta011+theta010
c2<- ph101+ph100+ph001+ph000
s3<-theta111+theta101+theta011+theta001
c3<-ph110+ph100+ph010+ph000
}
# prostate enlaragement with 3 readers
list( n111=46,n110=20,n101=27,n100=31,
n011=31,n010=79,n001=41,n000=1533)
# hyperparameters for pepe example
# Test1 is culture, 2 is Elisa, 3 is PCR
list(ap0=50,bp0=1000,n111=20,n101=4,n110=2,
n100=2,n011=4,n001=3,n010=2,n000=292)
# starting values
list(p=.5, h111=1,h110=1,h101=1,h100=1,
h011=1,h010=1,h001=1,h000=1,
y111=1,y110=1,y101=1,y100=1,
y011=1,y010=1,y001=1,y000=1)
BUGS CODE 7.7
model;
{
# Two tests T1 and T2 with ordinal values
# Dirichlet for theta
for(i in 1:3){ for (j in 1:3){h[i,j]~dgamma( zh[i,j],2)}}
sh<-sum(h[,])
for(i in 1:3){ for (j in 1:3){theta[i,j]<-h[i,j]/sh}}
for(i in 1:3){ for (j in 1:3){zh[i,j]<-y[i,j]+1}}
# Dirichlet for ph
for(i in 1:3){ for (j in 1:3){g[i,j]~dgamma( zg[i,j],2)}}
sg<-sum(g[,])
for(i in 1:3){ for (j in 1:3){ph[i,j]<-g[i,j]/sg}}
for(i in 1:3){ for (j in 1:3){zg[i,j]<- n[i,j] - y[i,j]+1}}
# dist of augmented data y[i,j]
for(i in 1:3){ for (j in 1:3){y[i,j]~dbin(w[i,j],n[i,j])}}
for(i in 1:3){ for (j in 1:3){w[i,j]<-d*theta[i,j]/(d*theta[i,j]+(1-d)*ph[i,j])}}
# posterior distribution of p
# the disease rate d is a mixure of p and p0
# p0 is the prior distribution of the disease rate
# p is the disease rate with a uniform prior
d<-.1*p+.9*p0
p0~dbeta(ap,bp)
p~dbeta(alphap,betap)
sy<-sum(y[,])
sn<-sum(n[,])
alphap<- sy+1
betap<- sn-sy+1
# ROC area T1
# Prob T1 = 1,2,3 given D=1
theta.1<-
theta[1,1]+theta[2,1]+theta[3,1]
theta.2<-
theta[1,2]+theta[2,2]+theta[3,2]
theta.3<-
theta[1,3]+theta[2,3]+theta[3,3]
# Prob T1=1,2,3 given D=0
ph.2<- ph[1,2]+ph[2,2]+ph[3,2]
ph.3<- ph[1,3]+ph[2,3]+ph[3,3]
ph.1<- ph[1,1]+ph[2,1]+ph[3,1]
T1<-T11+T12/2
T11<- theta.2*ph.1+theta.3*(ph.1+ph.2)
# T12 is the prob of a tie
T12<- theta.1*ph.1+theta.2*ph.2+theta.3*ph.3
# ROC area T2
# Prob T2=1,2,3 given D=1
theta1.<- theta[1,1]+theta[1,2]+theta[1,3]
theta2.<-
theta[2,1]+theta[2,2]+theta[2,3]
theta3.<-
theta[3,1]+theta[3,2]+theta[3,3]
# Prob T2=1,2,3 given D=0
ph2.<-
ph[2,1]+ph[2,2]+ph[2,3]
ph3.<-
ph[3,1]+ph[3,2]+ph[3,3]
ph1.<-
ph[1,1]+ph[1,2]+ph[1,3]
T2<- T21+T22/2
T21<-theta2.*ph1.+theta3.*(ph1.+ph2.)
# T22 is the Prob of a tie
T22<- theta1.*ph1.+theta2.*ph2.+theta3.*ph3.
}
# staging for melanoma
list(ap=20,bp=80, n=structure(.Data=c(45,26,28,25,20,20,46,30,46),
.Dim=c(3,3)))
# hypothetical
list(ap=64,bp=36, n=structure(.Data=c(30,40,50,1,30,60,1,1,70),
.Dim=c(3,3)))
# initial values
list(p=.5,
g=structure(.Data=c(1,1,1,1,1,1,1,1,1),.Dim=c(3,3)),
h=structure(.Data=c(1,1,1,1,1,1,1,1,1),.Dim=c(3,3)))
