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@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.
