[R] Bayesian estimate of prevalence with an imperfect test

[LianD]

Hi all!

I'm new to this forum so please excuse me if I don't conform perfectly to
the protocols on this board!

I'm trying to get an estimate of true prevalence based upon results from an
imperfect test. I have various estimates of se/sp which could inform my
priors (at least upper and lower limits even if with a uniform distribution)
and found the following code on this website..

http://www.lancs.ac.uk/staff/diggle/prevalence_estimation.R/

(the folllowing code has been cut and pasted directly from the web resource
- the only change I have made is to fill in my values for T, n, low/high
se/sp and the alpha beta for the distributions)

prevalence.bayes<-function(theta,T,n,lowse=0.5,highse=1.0,
   sea=1,seb=1,lowsp=0.5,highsp=1.0,spa=1,spb=1,ngrid=20,coverage=0.95) {


ibeta<-function(x,a,b) {
      pbeta(x,a,b)*beta(a,b)
      }
   ntheta<-length(theta)
   bin.width<-(theta[ntheta]-theta[1])/(ntheta-1)
   theta<-theta[1]+bin.width*(0:(ntheta-1))
   integrand<-array(0,c(ntheta,ngrid,ngrid))
   h1<-(highse-lowse)/ngrid
   h2<-(highsp-lowsp)/ngrid
   for (i in 1:ngrid) {
      se<-lowse+h1*(i-0.5)
      pse<-(1/(highse-lowse))*dbeta((se-lowse)/(highse-lowse),sea,seb)
      for (j in 1:ngrid) {
         sp<-lowsp+h2*(j-0.5)
         psp<-(1/(highsp-lowsp))*dbeta((sp-lowsp)/(highsp-lowsp),spa,spb)
         c1<-1-sp
         c2<-se+sp-1
         f<-(1/c2)*choose(n,T)*(ibeta(c1+c2,T+1,n-T+1)-ibeta(c1,T+1,n-T+1))
         p<-c1+c2*theta
         density<-rep(0,ntheta)
         for (k in 1:ntheta) {
            density[k]<-dbinom(T,n,p[k])/f
            }
         integrand[,i,j]<-density*pse*psp
         }
      }
   post<-rep(0,ntheta)
   for (i in 1:ntheta) {
      post[i]<-h1*h2*sum(integrand[i,,])
      }
   ord<-order(post,decreasing=T)
   mode<-theta[ord[1]]   
   take<-NULL
   prob<-0
   i<-0
   while ((prob<coverage/bin.width)&(i<ntheta)) {
      i<-i+1
      take<-c(take,ord[i])
      prob<-prob+post[ord[i]]
      }
   if (i==ntheta) {
      print("WARNING: range of values of theta too narrow")
      }
   interval<-theta[range(take)] 
  
list(theta=theta,post=post,mode=mode,interval=interval,coverage=prob*bin.width)
   }


n<-383
T<-97
ngrid<-25
lowse<-0.527
highse<-0.847
lowsp<-0.446
highsp<-0.851
sea<-4.38
seb<-2.93
spa<-3.18
spb<-3.54
theta<-0.001*(1:400)
coverage<-0.95
result<-prevalence.bayes(theta,T,n,lowse,highse,
sea,seb,lowsp,highsp,spa,spb,ngrid,coverage)

result$mode # 0.115
result$interval # 0.011 0.226
plot(result$theta,result$post,type="l",xlab="theta",ylab="p(theta)")

however, when I try to run the code I get the following error "Error in
while ((prob < coverage/bin.width) & (i < ntheta)) { : 
  missing value where TRUE/FALSE needed"

I have to admit I don't really understand what this error is telling me -
has anyone else ever used this code and would you mind letting me know what
I need to do to run it?

thanks everyone so much in advance for your help!

all the best and a happy new year!

Lian








--
View this message in context: http://r.789695.n4.nabble.com/Bayesian-estimate-of-prevalence-with-an-imperfect-test-tp4265595p4265595.html
Sent from the R help mailing list archive at Nabble.com.