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)-log10(df2$Dose);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
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.