power calculator for stepped wedge designs

steppedwedgeR.R

# Calculation of power for stepped wedge design
# Adapted from an excel-based calculator from Jim Hughes, jphughes@u.washington.edu
  # see http://faculty.washington.edu/peterg/Vaccine2006/articles/HusseyHughes.2007.pdf
  # http://faculty.washington.edu/jphughes/Power%20for%20stepped%20wedge.xls
  # http://faculty.washington.edu/jphughes/Power%20for%20stepped%20wedge%20(mean).xls
# Adapted by Eric Green, eric.green@duke.edu, @ericpgreen

# (1) enter or (2) import design matrix ===========================================================
  # skip option 1 if using option 2, and vice versa
  # power.p is the final object of interest for proportions
  # power.m is the final object of interest for means

# Option 1. enter design details ------------------------------------------------------------------
  # only use if equal allocation per round
  rounds <- 15           # number of observation rounds
  units <- 14            # number of units (clusters/individuals) to be randomized
  XatT1 <- 0             # will units receive intervention in first roundervation round (no=0; yes=1)
  # number of units receiving intervention each round (assumes equal allocation)
    units.XatTn <- ifelse(XatT1==0, 
                          units/(rounds-1),  
                          units/rounds)      

# create design matrix (do not modify) ------------------------------------------------------------
  d.mat <- data.frame(matrix(NA, nrow = units, ncol = rounds))
  d.mat[1:rounds] <- 0
  start.c <- ifelse(XatT1==0,2,1)
  start.r <- 1
  end.r <- start.r + (units.XatTn-1)
  round.loop <- ifelse(XatT1==0,rounds-1,rounds)
  l <- 1
  # begin loop
    while(l <= round.loop) {
      d.mat[start.r:end.r, (start.c):rounds] <- 1
      start.c <- start.c + 1
      start.r <- end.r + 1
      end.r <- end.r + units.XatTn
      l <- l + 1
    }

# Option 2. import specialized design matrix
  # forthcoming


# PROPORTIONS #####################################################################################

# enter more design details =======================================================================
  cluster <- 1                    # cluster-randomized=1; individual-randomized=0
  p1 <- 0.05                      # proportion (outcome) in control group
  p2 <- 0.025                     # proportion (outcome) in intervention group
  N <- units                      # do not modify: number of units
  T <- rounds                     # do not modify: number of time periods
  n <- 2                          # observations per cluster (ignore if not CRT)
  n <- ifelse(cluster==1, n, 1)   # do not modify
  k <- .15                        # coefficient of variation (usually between 0.15 and 0.40)
  k <- ifelse(cluster==1, k, 0)   # do not modify
  Za <- 1.96                      # 1.96 for alpha = 0.05

# calculations (do not modify) ====================================================================
  round.sum <- rowSums(d.mat)
  round.sum2 <- round.sum*round.sum
  unit.sum <- colSums(d.mat)
  unit.sum2 <- unit.sum*unit.sum
  U <- sum(round.sum)
  V <- sum(round.sum2)
  W <- sum(unit.sum2)
  t <- p1-p2
  d2 <- (((p1+p2)/2)*(1-(p1+p2)/2))/n
  z2 <- (p1*p1)*(k*k) 
  var.t1 <- (N*d2)*(d2+(T*z2))
  var.t2 <- ((N*U)-W)*d2
  var.t3 <- ((U*U)+(N*T*U)-(T*W)-(N*V))*z2
  var.t <- var.t1/(var.t2+var.t3)
  power1 <- sqrt((t*t)/var.t)
  power2 <- power1-Za
  power.p <- pnorm(power2, mean = 0, sd = 1)

# MEANS ###########################################################################################

# enter more design details =======================================================================
  cluster <- 1                    # cluster-randomized=1; individual-randomized=0
  m1 <- 20                        # mean (outcome) in control group
  m2 <- 10                        # mean (outcome) in intervention group
  s <- 25                         # standard deviation of outcome within site
  N <- units                      # do not modify: number of units
  T <- rounds                     # do not modify: number of time periods
  n <- 2                          # observations per cluster (ignore if not CRT)
  n <- ifelse(cluster==1, n, 1)   # do not modify
  k <- .15                        # coefficient of variation (usually between 0.15 and 0.40)
  k <- ifelse(cluster==1, k, 0)   # do not modify
  Za <- 1.96                      # 1.96 for alpha = 0.05

# calculations (do not modify) ====================================================================
  round.sum <- rowSums(d.mat)
  round.sum2 <- round.sum*round.sum
  unit.sum <- colSums(d.mat)
  unit.sum2 <- unit.sum*unit.sum
  U <- sum(round.sum)
  V <- sum(round.sum2)
  W <- sum(unit.sum2)
  t <- m1-m2
  d2 <- (s*s)/n
  z2 <- (m1*m1)*(k*k) 
  var.t1 <- (N*d2)*(d2+(T*z2))
  var.t2 <- ((N*U)-W)*d2
  var.t3 <- ((U*U)+(N*T*U)-(T*W)-(N*V))*z2
  var.t <- var.t1/(var.t2+var.t3)
  power1 <- sqrt((t*t)/var.t)
  power2 <- power1-Za
  power.m <- pnorm(power2, mean = 0, sd = 1)