# Slide 12: Power
# Suppose 50 in each group
CV <- qnorm(.975)
1-pnorm(CV - log(1.5)/sqrt(1/50+1/50))
pnorm(-1.96 - log(2/3)/sqrt(1/50+1/50))  # If Delta < 1

# Slide 13: Sample size
2*(qnorm(1-.05/2) + qnorm(.8))^2/log(1.5)^2
2*(qnorm(1-.05/2) + qnorm(.8))^2/log(2/3)^2   # Same as above

# Slides 22/23

source("http://myweb.uiowa.edu/pbreheny/7210/f15/notes/fun.R")
e <- function(t, a, lam, tstop=1) {
  tt <- pmin(t, tstop)
  p <- 1 - exp(-lam*t)/(lam*tt)*(exp(lam*tt)-1)
  a*tt*p/lam
}
lam <- log(2)/(2.5*1.5/12)

# Accrual rate = 50 / year
t <- seq(0, 2, len=199)
plot(t, e(t, a=50, lam=lam, tstop=1), type='l', lwd=3, las=1, bty='n', col=pal(3)[1], ylim=c(0,30), xlab="Time (years)", ylab="Expected patient-years")
lines(t, e(t, a=50, lam=lam, tstop=1.5), lwd=3, col=pal(3)[2])
lines(t, e(t, a=50, lam=lam, tstop=2), lwd=3, col=pal(3)[3])
abline(h=48/lam, col='gray60', lwd=2, lty=2)
toplegend(legend=c("T=1", "T=1.5", "T=2"), lwd=3, col=pal(3))

# Accrual rate = 40 / year
t <- seq(0, 2, len=199)
plot(t, e(t, a=40, lam=lam, tstop=1), type='l', lwd=3, las=1, bty='n', col=pal(3)[1], ylim=c(0,30), xlab="Time (years)", ylab="Expected patient-years")
lines(t, e(t, a=40, lam=lam, tstop=1.5), lwd=3, col=pal(3)[2])
lines(t, e(t, a=40, lam=lam, tstop=2), lwd=3, col=pal(3)[3])
abline(h=48/lam, col='gray60', lwd=2, lty=2)
toplegend(legend=c("T=1", "T=1.5", "T=2"), lwd=3, col=pal(3))