source("http://myweb.uiowa.edu/pbreheny/7210/f19/notes/fun.R")

# Slide 16
CV <- qnorm(.975)
1-pnorm(CV - log(1.5)/sqrt(1/50+1/50))       # 52% power at n=50
2*(qnorm(1-.05/2) + qnorm(.8))^2/log(1.5)^2  # need n=96 for 80% power
1-pnorm(CV - log(1.5)/sqrt(1/96+1/96))       # Told you
2*(qnorm(1-.05/2) + qnorm(.8))^2/log(2/3)^2  # Doesn't matter if you flip the ratio

# Slide 17
1-pnorm(CV - log(1.5)/sqrt(1/48+1/Inf))      # Now 80% power at n=48 if n=Inf for controls

# Slide 22-24: Setup
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)

# Slide 22: 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))

# What if the accrual rate was 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))

# Slide 23: Showing variability
N <- 50
t <- seq(0, 2, len=199)
plot(0, 0, type='n', las=1, bty='n', col=pal(3)[3], xlim=range(t), ylim=c(0,30), xlab="Time (years)", ylab="Expected patient-years")
A <- rpois(N, 1.5*50)
for (i in 1:N) {
  v <- runif(A[i], 0, 1.5)
  u <- rexp(A[i], lam)
  f <- numeric(length(t))
  for (j in 1:length(t)) {
    f[j] <- sum(pmin(pmax(t[j]-v, 0), u))
  }
  lines(t, f, col=rgb(0.5,0.5,0.5,alpha=0.5), lwd=0.75)
}
lines(t, e(t, a=50, lam=lam, tstop=1.5), lwd=3, col=pal(3)[2])
abline(h=48/lam, col='gray60', lwd=2, lty=2)