#tref     time at which mortalities estimated
#n1       total sample size, stratum 1
#n2       total sample size, stratum 2
#m1c      tref-year mortality, stratum 1 control
#m2c      "          "                 2  "
#r1       % reduction in m1c by intervention, stratum 1
#r2       % reduction in m2c by intervention, stratum 2
#accrual  duration of accrual period
#tmin     minimum follow-up time
#alpha    type I error
#pr       set to T to print intermediate results

ciapower <- function(tref,   
					n1,     
					n2,     
					m1c,    
					m2c,    
					r1,     
					r2,     
					accrual,
					tmin,   
					alpha=.05,  
					pr=TRUE) { 

#Find mortality in intervention groups
if(m1c>1 | m2c>1) stop("m1c and m2c must be fractions")
m1i <- (1-r1/100)*m1c
m2i <- (1-r2/100)*m2c

if(pr) {
  cat("\nAccrual duration:",accrual,"y  Minimum follow-up:",tmin,"y\n")
  cat("\nSample size Stratum 1:",n1,"  Stratum 2:",n2,"\n")
  cat("\nAlpha=",alpha,"\n")
  d <- list(c("Stratum 1","Stratum 2"), c("Control","Intervention"))
  m <- cbind(c(m1c,m2c),c(m1i,m2i))
  dimnames(m) <- d
  cat("\n",tref,"-year Mortalities\n",sep=""); print(m)
}

#Find exponential hazards for all groups
lam1c <- -logb(1-m1c)/tref
lam2c <- -logb(1-m2c)/tref
lam1i <- -logb(1-m1i)/tref
lam2i <- -logb(1-m2i)/tref

if(pr) {
  lam <- cbind(c(lam1c,lam2c),c(lam1i,lam2i))
  dimnames(lam) <- d
  cat("\nHazard Rates\n"); print(lam)
}

#Find probability that a subject will have her event observed during
#the study, for all groups
tmax <- tmin+accrual
p1c <- 1-1/accrual/lam1c*(exp(-tmin*lam1c)-exp(-tmax*lam1c))
p2c <- 1-1/accrual/lam2c*(exp(-tmin*lam2c)-exp(-tmax*lam2c))
p1i <- 1-1/accrual/lam1i*(exp(-tmin*lam1i)-exp(-tmax*lam1i))
p2i <- 1-1/accrual/lam2i*(exp(-tmin*lam2i)-exp(-tmax*lam2i))

if(pr) {
  p <- cbind(c(p1c,p2c), c(p1i,p2i))
  dimnames(p) <- d
  cat("\nProbabilities of an Event During Study\n")
  print(p)
}

#Find expected number of events, all groups
m1c <- p1c*n1/2
m2c <- p2c*n2/2
m1i <- p1i*n1/2
m2i <- p2i*n2/2

if(pr) {
  m <- cbind(c(m1c,m2c), c(m1i,m2i))
  dimnames(m) <- d
  cat("\nExpected Number of Events\n")
  print(round(m,1))
}

#Find expected value of observed log hazard ratio
delta <- logb((lam1i/lam1c)/(lam2i/lam2c))
if(pr) cat("\nRatio of hazard ratios:",format(exp(delta)),"\n")

#Find its variance
v <- 1/m1c + 1/m2c + 1/m1i + 1/m2i
sd <- sqrt(v)
if(pr) cat("Standard deviation of log ratio of ratios:",format(sd),"\n")

z <- -qnorm(alpha/2)
#if(pr) cat("\nCritical value:",format(z),"\n")

c(Power = 1 - ( pnorm(z - abs(delta)/sd) - pnorm(-z - abs(delta)/sd) ) )
}