[R] Bayesian estimate of prevalence with an imperfect test

[Bert Gunter]

Lian:

I doubt whether anyone on this list would be willing to wade through
your code to track this down, although some smart folks may be able to
guess what the issue is (alas for you, I am not one of them!). So now
is as good a time as any to start learning about R's debugging tools.

First, you may wish to set
options(error=recover)  ## see ?options

Then read up on:
?trace
?browser
?debug

Unless, you receive insight from a smart folk, you need to step
through your code, starting at some "suitable" entry point, and
examine the state of the variables in the error message to see
where/how the missing value occurs: Are there missings in your data?
-- Did you fail to specify a required argument in your function call?,
etc.

You may also wish to search (see ?RSiteSearch and/or the "sos"
package) on "debug" for other debugging tools (e.g. in the "debug"
package.

R generally rewards those who suffer through such efforts.

HTH

Cheers,
Bert

On Thu, Jan 5, 2012 at 6:09 AM, LianD <liandoble at hotmail.com> wrote:
> 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.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm