1
$\begingroup$

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

X0=Log10 of the reciprocal of lowest vaccine dose group *i* at which all animals are protected; *d* = difference between the Log10 dilution intervals;ni is the number of animals used at each vaccine dose group *i* and *ri* is the number of protected animals out of ni

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.

$\endgroup$

migrated from stackoverflow.com Sep 14 '17 at 12:31

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

  • $\begingroup$ 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? $\endgroup$ – Thies Heidecke Sep 14 '17 at 13:13
  • $\begingroup$ 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. $\endgroup$ – R.P. Tamil Selvan Sep 15 '17 at 9:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.