# bangladesh cholera vaccine trial  jh 2015.11.09

# attrition rate

# exp(-730 * rate ) = 9/16

a.rate = -log(9/16)/730

N=50000 # trial and error

n.c = N * exp( - ??? *(1:730) )
plot(1:730,n.c, ylim=c(0,N))
polygon(c(1:730,730,1),c(n.c,0,0),col="yellow" )

P1000Y=365*1000
100000*1.6/P1000Y

observed.attack.rate=106/(????????)

PD = sum(n.c); PD    # person days

assumed.attack.rate=1.6/(1000*365)

c(assumed.attack.rate, observed.attack.rate)

E.attacks = assumed.attack.rate * ?????????
E.attacks

##  assuming for now, Y_0 and Y_1 are Poisson
##  i.e. no extra-Poisson (between clsuter) variation

d.crit = -1.96 * sqrt(?????? + ?????? )

#  alt 

Efficacy = 0.65; Coverage=0.65

alt.Rate.Ratio = 1-Efficacy*Coverage;
alt.Rate.Ratio

alt.E.attacks = ???????? * E.attacks

c.crit.z.alt=(d.crit-
       (alt.E.attacks-E.attacks) )/
         sqrt(alt.E.attacks+E.attacks)
c.crit.z.alt
pnorm(c.crit.z.alt)

## in a diagram

E.diff.n.cases.null = E.attacks - E.attacks
V.diff.n.cases.null = E.attacks + E.attacks

E.diff.n.cases.alt  = alt.E.attacks - E.attacks
V.diff.n.cases.alt  = alt.E.attacks + E.attacks

possible.diffs = floor( E.diff.n.cases.alt - 
         3*sqrt(V.diff.n.cases.alt)
                      ):
                 ceiling( E.diff.n.cases.null + 
         3*sqrt(V.diff.n.cases.null)
                      )

prob.mass.fn.null = dnorm(possible.diffs,
                          E.diff.n.cases.null,
                      sqrt(V.diff.n.cases.null)
                         )

prob.mass.fn.alt = dnorm(possible.diffs,
                          E.diff.n.cases.alt,
                      sqrt(V.diff.n.cases.alt)
                         )

plot(possible.diffs,prob.mass.fn.null,
  type="h",col="green",
  ylim=c(-max(prob.mass.fn.alt),
          max(prob.mass.fn.null)),
   xlab="N.cases.vaccine.arm - N.cases.ctl.arm",
   ylab="Prob. of observing.This difference")
abline(v=0-1.96*sqrt(V.diff.n.cases.null),
col="green",lwd=3)
text(1*sqrt(V.diff.n.cases.null),
     2/3*max(prob.mass.fn.null),
     "NULL",col="green",adj=c(-0.5,0))
SD=sqrt(V.diff.n.cases.null)
h = exp(-1/2)*max(prob.mass.fn.null)
segments( 0,h,floor(SD),h,
         col="green",lwd=2)
text(SD/2,h,toString(round(SD,1)),
   col="green",adj=c(0.5,1.25) )
text(0,0,paste("0 = ",
               toString(round(E.attacks,1)),
               " - ",
               toString(round(E.attacks,1)),
               sep=""),
col="green",
     adj=c(0,1.25) )
points(possible.diffs+0.5,-prob.mass.fn.alt,
  type="h", col="red")
text(E.diff.n.cases.alt-1*sqrt(V.diff.n.cases.alt),
     -2/3*max(prob.mass.fn.null),
     "ALT",col="red",adj=c(1.5,0))

text(E.diff.n.cases.alt,0,
     paste(
     toString(alt.Rate.Ratio),
     " 
     x ",
     toString(round(E.attacks,1)),
     " 
     -     
      ",
     toString(round(E.attacks,1)),
     " 
     = ",
     toString(round(E.diff.n.cases.alt,1)),
     sep=""),
     col="red",
     adj=c(1,-0.1) )
text(E.diff.n.cases.alt,0,"|",col="red")     
     






c=c(65,106); PD=c(41809947,39327744)
Rate=c/PD
round(Rate[1]/Rate[2],2)
efficacy=round(100*(1-Rate[1]/Rate[2]),0)
log.ratio=log(Rate[1]/Rate[2]) 

SE.log = sqrt(sum(1/c))

round(100*(1-exp(log.ratio+c(-1.96,0,1.96)*SE.log)),0)

round(100*(1-exp(log.ratio+c(-1.645,0,1.645)*SE.log)),0)




N1=94675
a.rate = -log(38866/N1)/730
n.c = N1 * exp( - a.rate*(1:730) )
prob = n.c/sum(n.c)
#plot(1:730,prob, ylim=c(0,max(prob)))
t.events1 = sort( sample(1:730,65, prob=prob) )
n.riskset1 = n.c[t.events1]
Sum = sum(1/(n.riskset1^2))
Var.logR1 = (1/0.0011^2) * (0.9989^2) * Sum

N0=80056 
a.rate = -log(37303/N0)/730
n.c = N0 * exp( - a.rate*(1:730) )
prob = n.c/sum(n.c)
t.events0 = sort( sample(1:730,106, prob=prob) )
n.riskset0 = n.c[t.events0]
Sum = sum(1/(n.riskset0^2))
Var.logR0 = (1/0.0011^2) * (0.9989^2) * Sum


Var.logRR = Var.logR1 + Var.logR0

CI.RR = round( 100*exp( log(0.0011/0.0019)+
            c(-1.96,0,1.96)*sqrt(Var.logRR) 
          ), 0)
100-CI.RR