Estimation of prevalence distribution: 2 tests, multiple populations
WinBUGS 1.4 code to accompany the paper entitled "Branscum AJ, Gardner IA, Johnson WO. Bayesian modeling of animal- and herd-level prevalences. Prev Vet Med. 2004 Dec 15;66(1-4):101-112".
Code prepared by Adam Branscum, November 18, 2003
branscum@ms.uky.edu
Departments of Biostatistics and Statistics
University of Kentucky
Example from section 3.1.1.
Estimate the prevalence distribution for Brucellosis in cattle in Mexico.
Data source: Hanson TE, Johnson WO, Gardner IA. Hierarchical models for estimating herd prevalence and test accuracy in the absence of a gold standard. J Agric Biol Environ Statist. 2003 Jun 1; 8(2), 223–239.
An example of a joint sequential test.
Test 1 is BAPA so Sebapa and Spbapa are the sensitivity and specificity for BAPA and test 2 is Rivanol so Seriv and Spriv are the sensitivity and specificity of Rivanol.
model
{
for(i in 1:K)
{
x[i, 1:3] ~ dmulti(p[i, 1:3], n[i])
p[i,1] <- pi[i]*(Sebapa*Seriv+covDp) + (1-pi[i])*((1-Spbapa)*(1-Spriv)+covDn)
p[i,2] <- pi[i]*(Sebapa*(1-Seriv)-covDp) + (1-pi[i])*((1-Spbapa)*Spriv-covDn)
p[i,3] <- 1-p[i,1]-p[i,2]
pi[i] ~ dbeta(alpha, beta)
}
Sebapa ~ dbeta(88.3,1.9) ## Mode=0.99, 95% sure > 0.95
Spbapa ~ dbeta(13.3,6.3) ## Mode=0.70, 95% sure > 0.50
Seriv ~ dbeta(130.7,15.4) ## Mode=0.90, 95% sure > 0.85
Spriv ~ dbeta(99.7,6.2) ## Mode=0.95, 95% sure > 0.90
alpha <- mu*psi
beta <- psi*(1-mu)
mu ~ dbeta(16.9,48.6) ## Mode=0.25; 95% sure < 0.35
psi ~ dgamma(7.23, 1.28)
ls <- (Sebapa-1)*(1-Seriv)
lc <- (Spbapa-1)*(1-Spriv)
us <- min(Sebapa,Seriv)-Sebapa*Seriv
uc <- min(Spbapa, Spriv) - Spbapa*Spriv
covDn ~ dunif(lc, uc)
covDp ~ dunif(ls, us)
rhoDp <- covDp / sqrt(Sebapa*(1-Sebapa)*Seriv*(1-Seriv))
rhoDn <- covDn / sqrt(Spbapa*(1-Spbapa)*Spriv*(1-Spriv))
pi21 ~ dbeta(alpha,beta)
a[1] <- 1-step(pi21-0.50)
a[2] <- 1-step(pi21-0.10)
a[3] <- step(pi21-0.90)
a[4] <- step(pi21-0.75)
}
list( Sebapa=0.99, Spbapa=0.70, Seriv=0.90, Spriv=0.95, mu=0.25, psi=7)
list(K=20, n=c(67,12,22,71,18,80,147,45,37,118,12,62,11,20,56,60,26,54,17,17),
x=structure(.Data=c(4,2,61,
1,1,10,
2,2,18,
2,0,69,
2,0,16,
2,0,78,
80,5,62,
4,5,36,
10,0,27,
10,12,96,
1,0,11,
1,1,60,
3,0,8,
3,7,10,
38,0,18,
1,0,59,
3,0,23,
13,3,38,
7,0,10,
1,1,15),.Dim=c(20,3)))
node |
mean |
sd |
MC error |
2.50% |
median |
97.50% |
start |
sample |
Sebapa |
0.972 |
0.01888 |
2.77E-04 |
0.9255 |
0.9761 |
0.9966 |
5001 |
50000 |
Seriv |
0.9176 |
0.01896 |
1.19E-04 |
0.8772 |
0.9188 |
0.9513 |
5001 |
50000 |
Spbapa |
0.9296 |
0.01477 |
1.92E-04 |
0.8996 |
0.93 |
0.9574 |
5001 |
50000 |
Spriv |
0.9426 |
0.02129 |
1.43E-04 |
0.8926 |
0.946 |
0.9754 |
5001 |
50000 |
a[1] |
0.9148 |
0.2792 |
0.001232 |
0 |
1 |
1 |
5001 |
50000 |
a[2] |
0.4314 |
0.4953 |
0.002647 |
0 |
0 |
1 |
5001 |
50000 |
a[3] |
0.00208 |
0.04556 |
2.10E-04 |
0 |
0 |
0 |
5001 |
50000 |
a[4] |
0.01478 |
0.1207 |
5.01E-04 |
0 |
0 |
0 |
5001 |
50000 |
alpha |
0.6935 |
0.2334 |
0.003044 |
0.3263 |
0.6641 |
1.23 |
5001 |
50000 |
beta |
2.961 |
0.9713 |
0.01053 |
1.393 |
2.849 |
5.177 |
5001 |
50000 |
mu |
0.1918 |
0.03431 |
2.81E-04 |
0.1294 |
0.1902 |
0.2636 |
5001 |
50000 |
psi |
3.654 |
1.161 |
0.01339 |
1.77 |
3.526 |
6.284 |
5001 |
50000 |
rhoDn |
0.3981 |
0.168 |
0.00198 |
0.05055 |
0.4086 |
0.6859 |
5001 |
50000 |
rhoDp |
0.3119 |
0.2146 |
0.002489 |
-0.01402 |
0.2886 |
0.774 |
5001 |
50000 |