# Bayesian confidence interval of ED50

I want to calculate Bayesian confidence interval of PD 50 (median protective dose) values as shown in this paper DOI: 10.1016/j.vaccine.2006.12.049. My data is

library(dplyr)
Dose=c(1,1/4,1/16)
Protected=c(5,3,3)
Total=c(rep(5,3))
df=as.data.frame(cbind(Dose,Protected,Total));df


Then, I calculated Log10 (PD50) using the Karber formula

The r code used for calculating the PD50 is follows

df2=df %>% mutate (PropProt=Protected/Total);df2
d=log10(df2$Dose[1])-log10(df2$Dose[2]);d
X0=log10(min(df2$Dose[df2$PropProt==1.0]));X0
sum=sum(df2$PropProt);sum logPD=X0-d/2+d*sum;logPD PD=10^logPD;PD  For my data the Log10 PD50 is calculated as 1.02 and PD 50 as 10.56. Then, i did random sampling from posterior distribution with uniform prior distribution for each proportion library(LearnBayes) set.seed(27021979) a=1; b=1 # parameters of beta prior n=5; y=5 # sample size and number of yes’s in sample a1=a+y; b1=b+n-y # parameters of beta posterior a1/(a1+b1)# a point estimate is given by the posterior mean p1=rbeta(5000,a1,b1) a=1; b=1 # parameters of beta prior n=5; y=3 # sample size and number of yes’s in sample a2=a+y; b2=b+n-y # parameters of beta posterior a2/(a2+b2)# a point estimate is given by the posterior mean p2=rbeta(5000,a2,b2) a=1; b=1 # parameters of beta prior n=5; y=3 # sample size and number of yes’s in sample a3=a+y; b3=b+n-y # parameters of beta posterior a3/(a3+b3)# a point estimate is given by the posterior mean p3=rbeta(5000,a3,b3)  then for each iteration i calculated the PD50 value using the formula shown above. P=data.frame(p1,p2,p3);head(P) P$Prop=with(P, p1+p2+p3);head(P)
P$X0=X0;head(P) P$d=d;head(P)
P$c=d/2;head(P) P$e=with(P,d*Prop);head(P)
P$logPD=with(P,X0-c+e);head(P) P$PD=with(P,10**logPD);head(P)


from these PD50 values, when i calculated the 95% confidence interval

t.test(P$PD, conf.level = 0.95)$conf.int


or

library(TeachingDemos)
z.test(P$PD,mu = mean(P$PD),stdev=sd(P\$PD))


95 percent confidence interval: 8.513967 - 8.692223, the interval is very narrow compared to wider interval ( 4.27- 21.82) as mentioned in the publication.

The method followed in the publication is "estimation of uncertainty around [ri/ni] by combining uninformative beta (1,1) prior with the data, (likelihood) posterior distributions were calculated. The uncertainty distribution of the PD50 value was obtained by repeatedly drawing randomly from these posterior distributions followed by, at each iteration (n=5000), the calculation of PD50 according to the formula described above. the histograms of all these PD values, calculated at each iteration, represent the uncertainty of the estimates and summary statistics like mean and percentiles were given.

I am not able to decode what went wrong in this analysis.

I highly appreciate the help in this method.

## migrated from stackoverflow.comSep 14 '17 at 12:31

This question came from our site for professional and enthusiast programmers.

• Not a direct answer to your question, but if you're going the bayesian route anyway, wouldn't credible intervals be a better fit than confidence intervals? – Thies Heidecke Sep 14 '17 at 13:13
• the median values are calculated using a karber formula, which is having proportion as one of parameter (ri/ni), these proportions are derived by sampling the posterior. but i think credible intervals are still more narrower than confidence interval, though i have not attempted to calculate(perhaps i dont know how to do that in r). i just want to know the methodology how it could be reproduced as seen in the publication. – R.P. Tamil Selvan Sep 15 '17 at 9:59