WinBUGS examples for Bayesian prevalence estimation for imperfect tests with multiple pool sizes

Example 1 (Truckee River data)
---------------------------------------------------------------------------

A first approach directly parametrizes the prevalence ...

model
   {
      for (i in 1:N) 
        {
         y[i] ~ dbern(q[i])
         q[i] <- Se*(1 -pow(1-Pi,n[i])) + (1-Sp)*(pow(1-Pi,n[i]))
        }
      Se ~ dbeta(48, 12) 
      Sp ~ dbeta(54, 6)
      Pi ~ dbeta(1, 1)
   }


# Truckee data
list(
y = c(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),
n = c(3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 2, 4, 4, 5),
N = 27)


# inits
list(Se=.5,Sp=.5,Pi=.5)


--------------------------------

This second way uses a logit transformation that allows more general models.  

model
   {
      for (i in 1:N) 
        {
         y[i] ~ dbern(q[i])
         q[i] <- Se*(1 -pow(1-Pi[i],n[i])) + (1-Sp)*(pow(1-Pi[i],n[i]))
         logit(Pi[i]) <- alpha0 
        }
      Se ~ dbeta(48, 12) 
      Sp ~ dbeta(54, 6)
      alpha0 ~ dlogis(0, 1.0)
      Pi0 <- exp(alpha0)/(1 +exp(alpha0))
   }


# Truckee data
list(
y = c(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),
n = c(3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 2, 4, 4, 5),
N = 27)


# inits
list(Se=.5,Sp=.5,alpha0=0.)

---------------------------------------------------------------------------

Example 2 (Henry's Lake data)

Here four sets of initialization parameters are used so that multiple MCMC chains can be run.

model
   {
      for (i in 1:N) 
        {
         y[i] ~ dbern(q[i])
         q[i] <- Se*(1 -pow(1-Pi[i],n[i])) + (1-Sp)*(pow(1-Pi[i],n[i]))
         logit(Pi[i]) <- alpha0 
        }
      Se ~ dbeta(450, 50) 
      Sp ~ dbeta(96, 24)
      alpha0 ~ dlogis(0, 1.0)
      Pi0 <- exp(alpha0)/(1 +exp(alpha0))
   }


# Henry's Lake 1996 
list(
y = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1),
n = c(5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5),
N = 12)

# inits
list(Se=.2,Sp=.2,alpha0=-1.)
list(Se=.2,Sp=.8,alpha0=1.)
list(Se=.8,Sp=.2,alpha0=1.)
list(Se=.8,Sp=.8,alpha0=-1.)

---------------------------------------------------------------------------

Example 3 (Multiple year hypothetical data where prevalence varies between years)

model
   {
      for (i in 1:N) 
        {
         y[i] ~ dbern(q[i])
         q[i] <- Se*(1 -pow(1-Pi[i],n[i])) + (1-Sp)*(pow(1-Pi[i],n[i]))
         logit(Pi[i]) <- alpha[1]*x1[i] + alpha[2]*x2[i] + alpha[3]*x3[i]
        }
      Se ~ dbeta(90, 10) 
      Sp ~ dbeta(90, 10)
      alpha[1] ~ dlogis(0, 1.0)
      alpha[2] ~ dlogis(0, 1.0)
      alpha[3] ~ dlogis(0, 1.0)
      Prob[1] <- exp(alpha[1])/(1 +exp(alpha[1]))
      Prob[2] <- exp(alpha[2])/(1 +exp(alpha[2]))
      Prob[3] <- exp(alpha[3])/(1 +exp(alpha[3]))

      Pmean   <- (Prob[1]+Prob[2]+Prob[3])/3

      diff05 <- step(abs(Prob[1]-Pmean)-.05) +step(abs(Prob[2]-Pmean)-.05) +step(abs(Prob[3]-Pmean)-.05)
      diff10 <- step(abs(Prob[1]-Pmean)-.10) +step(abs(Prob[2]-Pmean)-.10) +step(abs(Prob[3]-Pmean)-.10)
      diff15 <- step(abs(Prob[1]-Pmean)-.15) +step(abs(Prob[2]-Pmean)-.15) +step(abs(Prob[3]-Pmean)-.15)
      d05gt0 <- step(diff05-.01)
      d10gt0 <- step(diff10-.01)
      d15gt0 <- step(diff15-.01)

   }


# Hypothetical data, all years
list(
y = c(1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1,
      0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0,
      0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1),
n = c(5, 5, 5, 4, 5, 4, 1, 2, 5, 5, 6, 5,
      5, 2, 5, 5, 5, 5, 5, 5, 4, 5, 5, 6,
      8, 5, 5, 3, 5, 4, 1, 3, 5, 4, 6),
x1 = c(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),
x2 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
x3 = 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,
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
N = 35)


# inits
list(Se=.8,Sp=.8,alpha=c(-1.,1.,1.))


---------------------------------------------------------------------------

Example 4 (Multiple year hypothetical data where sensitivity, specificity, and prevalence vary between years)

model
   {
      for (i in 1:N) 
        {
         y[i] ~ dbern(q[i])
         Se[i] <- Sei[1]*x1[i] + Sei[2]*x2[i] + Sei[3]*x3[i]
         Sp[i] <- Spi[1]*x1[i] + Spi[2]*x2[i] + Spi[3]*x3[i]
         q[i] <- Se[i]*(1 -pow(1-Pi[i],n[i])) + (1-Sp[i])*(pow(1-Pi[i],n[i]))
         logit(Pi[i]) <- alpha[1]*x1[i] + alpha[2]*x2[i] + alpha[3]*x3[i]
        }

      Sei[1] ~ dbeta(70, 30) 
      Sei[2] ~ dbeta(80, 20) 
      Sei[3] ~ dbeta(90, 10) 
      Spi[1] ~ dbeta(70, 30)
      Spi[2] ~ dbeta(80, 20)
      Spi[3] ~ dbeta(90, 10)
      alpha[1] ~ dlogis(0, 1.0)
      alpha[2] ~ dlogis(0, 1.0)
      alpha[3] ~ dlogis(0, 1.0)
      Prob[1] <- exp(alpha[1])/(1 +exp(alpha[1]))
      Prob[2] <- exp(alpha[2])/(1 +exp(alpha[2]))
      Prob[3] <- exp(alpha[3])/(1 +exp(alpha[3]))

      Pmean   <- (Prob[1]+Prob[2]+Prob[3])/3

      diff05 <- step(abs(Prob[1]-Pmean)-.05) +step(abs(Prob[2]-Pmean)-.05) +step(abs(Prob[3]-Pmean)-.05)
      diff10 <- step(abs(Prob[1]-Pmean)-.10) +step(abs(Prob[2]-Pmean)-.10) +step(abs(Prob[3]-Pmean)-.10)
      diff15 <- step(abs(Prob[1]-Pmean)-.15) +step(abs(Prob[2]-Pmean)-.15) +step(abs(Prob[3]-Pmean)-.15)
      d05gt0 <- step(diff05-.01)
      d10gt0 <- step(diff10-.01)
      d15gt0 <- step(diff15-.01)

   }


# Hypothetical data, all years
list(
y = c(1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1,
      0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0,
      0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1),
n = c(5, 5, 5, 4, 5, 4, 1, 2, 5, 5, 6, 5,
      5, 2, 5, 5, 5, 5, 5, 5, 4, 5, 5, 6,
      8, 5, 5, 3, 5, 4, 1, 3, 5, 4, 6),
x1 = c(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),
x2 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
x3 = 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,
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
N = 35)


# inits
list(Sei=c(.7,.8,.9),Spi=c(.9,.8,.7),alpha=c(-1.,1.,1.))