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.
Bert Gunter
2012-Jan-05 15:51 UTC
[R] Bayesian estimate of prevalence with an imperfect test
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
Thanks Bert I've been through my variables again and have managed to get the code working - shouldn't have tried to deal with it at the end of the day yesterday! all the best Lian -- View this message in context: http://r.789695.n4.nabble.com/Bayesian-estimate-of-prevalence-with-an-imperfect-test-tp4265595p4268526.html Sent from the R help mailing list archive at Nabble.com.