### The appendix contains R functions which may be used to generate data, 
### apply a treatment effect, and calculate conditional power.  Two
### examples are provided.

#######################################################################
### EXAMPLE 1 -- Average Conditional Power of 2-Sided Tests --      ###
###              m = 2, n = 3, normal null data, constant treatment ###
###              effect = 1.5, size of the test = 2/10, 25 reps     ###
#######################################################################
### Generate null data
null.1 <- nulldata(m = 2, n = 3, reps = 25, dist = 'normal')

### Add a constant treatment effect of 1.5 to the treatment group
trtd.1 <- constte(nulldata = null.1, m = 2, delta = 1.5)

### Calculate the exact conditional power for each replication for 
### the permutation t test, Student's t test, and the WMW test
library(coin)
power.perm.1 <- cp.perm(trtd.1, c = 2, alpha = 2/10)
power.stud.1 <- cp.stud(trtd.1, c = 2, alpha = 2/10)
power.wmw.1 <- cp.wmw(trtd.1, c = 2, alpha = 2/10)

### Average conditional power across the replications
c(mean(power.perm.1), mean(power.stud.1), mean(power.wmw.1))

#######################################################################
### EXAMPLE 2 -- Average Conditional Power of 2-Sided Tests --      ###
###              m = 5, n = 5, skewed-bimodal null data, exp trtmt  ###
###              effect = 2.0, size of the test = 14/252, 5 reps    ###
#######################################################################
### Generate null data
null.2 <- nulldata(m = 5, n = 5, reps = 5, dist = 'mixture')

### Add an exponentially dist'd (mean = 2.0) trtmt effect to trtmt gp
trtd.2 <- expte(nulldata = null.2, m = 5, delta = 2.0)

### Calculate the exact conditional power for each replication for 
### the permutation t test, Student's t test, and the WMW test
library(coin)
power.perm.2 <- cp.perm(trtd.2, c = 5, alpha = 14/252)
power.stud.2 <- cp.stud(trtd.2, c = 5, alpha = 14/252)
power.wmw.2 <- cp.wmw(trtd.2, c = 5, alpha = 14/252)

### Average conditional power across the replications
c(mean(power.perm.2), mean(power.stud.2), mean(power.wmw.2))

############################## FUNCTIONS ##############################
nulldata <- function(m, n, reps, dist)
{
  # generate the null data for use in a conditional power study
  # additional replications are appended column-wise
  # m    ... size of the control group
  # n    ... size of the treatment group
  # reps ... number of replications
  # dist ... distribution of null data (all have mean = 0, var = 1)
  
  # define function for the null distribution
  distrib <- switch(dist,
                    exponential = function(x) rexp(x,1) - 1,
                    mixture = function(x) {			    
                      ind <- runif(x)				    
                      a <- ifelse(ind <= 0.8,1,0)		    
                      b <- rep(1,x)-a				    
                      v <- a*rnorm(x,0,1) + b*rnorm(x,5,2.25) 
                      v <- (v-1)*sqrt(1/5.81)			    
                      return(v)					    
                    },						    
                    normal = function(x) rnorm(x),		    
                    uniform = function(x) runif(x, min=-sqrt(12)/2, max=sqrt(12)/2))
  
  N <- m + n						
  cnm <- choose(N,m)					
  nullorder <- rbind(combn(N,m),combn(N,n)[,cnm:1])	
  nulldata <- matrix(distrib(N*reps),N,reps)		
  mc <- matrix(numeric(cnm*reps*N),N,cnm*reps)		
  for (i in 1:ncol(nulldata)) {			
    mc[,(cnm*(i-1)+1):(cnm*i)] <- 	
      matrix(nulldata[,i][nullorder], N, cnm)	
  }                                                    
  return(mc)						
}

constte <- function(nulldata, m, delta)
{
  # apply a constant additive shift treatment effect to the null data
  # for the treatment group
  # nulldata ... a matrix generated by function 'nulldata'
  # m        ... size of the control group
  # delta    ... constant shift
  
  n <- dim(nulldata)[1] - m               
  N <- n + m                              
  reps <- dim(nulldata[2])/choose(N,m)    
  nulldata[(m+1):N,] <- nulldata[(m+1):N,] + delta
  return(nulldata)                                 
}


expte <- function(nulldata, m, delta)
{
  # apply a non-constant exponentially distributed treatment effect to 
  # the null data for the treatment group
  # nulldata ... a matrix generated by function 'nulldata'
  # m        ... size of the control group
  # delta    ... 
  
  n <- dim(nulldata)[1] - m                              # trtmt group size
  N <- n + m                                             # total sample size
  cnm <- choose(N,m)                                     # size of perm dist
  reps <- dim(nulldata)[2]/cnm                           # number of reps
  trtef <- matrix(rexp(N*reps,rate=(1/delta)), N, reps)  # generate trtmt eff 
  nullorder <- rbind(combn(N,m),combn(N,n)[,cnm:1])      # perm distribution
  trtef.rep <- matrix(numeric(cnm*reps*N),N,cnm*reps)    # matrix of 0s
  for (i in 1:ncol(trtef)) {                             # fill trtef.rep 
    trtef.rep[,(cnm*(i-1)+1):(cnm*i)] <-                 # with trtmt effect 
      matrix(trtef[,i][nullorder], N, cnm)	
  }      
  nulldata[(m+1):N,] <- nulldata[(m+1):N,] + trtef.rep[(m+1):N,]  # add trtmt
  return(nulldata)
}

fun_p <- function(z) 
{
  # calculate the exact p-value from the permutation test
  # this function requires package coin
  # z ...two-column dataframe (the columns are (1) the data
  #      and (2) a factor indicating group membership)
  
  pvalue(independence_test(x~f, data=z, distribution = "exact"))
}

fun_df <- function(matrix, m, n)
{
  # given a matrix fun_df returns a list of length = ncols
  # such that each element is a data frame with m "a"s and n "b"s
  # matrix ...data matrix
  # m      ...size of control group
  # n      ...size of treatment group
  
  apply(matrix, 2, function(x)
  {
    data.frame(x, f=factor(c(rep("a",m),rep("b",n))))
  })
}

fun_st <- function(x, m, n)
{
  # run Student's t-test on a vector where the first m elements were
  # assigned to the control group
  # x ...vector of data
  # m ...size of control group
  # n ...size of treatment group
  
  t.test(x[1:m], x[m+1:(m+n)], var.equal=TRUE)$p.val
}

cp.perm <- function(data, c, alpha)
{
  # run and calculate conditional power for the permutation t test
  # data  ...null data with treatment effect applied via constte or expte
  # c     ...size of control group
  # alpha ...nominal size of the test
  
  t <- dim(data)[1] - c   
  N <- t + c              
  cnm <- choose(N,c)        
  reps <- dim(data)[2]/cnm    
  
  dat.df <- fun_df(matrix = data, m = c, n = t)
  perm.t <- lapply(dat.df, fun_p)
  permt <- unlist(perm.t)
  permt.mat <- matrix(permt, cnm, reps)
  power <- function(pvals)
  {
    length(which(pvals <= alpha + 1e-8))/cnm
  }    # 1e-8 is added to prevent rounding error
  power.perm <- apply(permt.mat, 2, power)
  return(power.perm)
}

cp.wmw <- function(data, c, alpha)
{
  # run and calculate conditional power for the WMW test
  # data  ...null data with treatment effect applied via constte or expte
  # c     ...size of control group
  # alpha ...nominal size of the test
  
  t <- dim(data)[1] - c       
  N <- t + c                  
  cnm <- choose(N,c)          
  reps <- dim(data)[2]/cnm    
  
  rank.data <- apply(data, 2, rank)
  rank.data.df <- fun_df(rank.data, m = c, n = t)
  wmw.t <- lapply(rank.data.df, fun_p)
  wmw <- unlist(wmw.t)
  wmw.mat <- matrix(wmw, cnm, reps)
  power <- function(pvals)
  {
    length(which(pvals <= alpha + 1e-8))/cnm
  }    # 1e-8 is added to prevent rounding error
  power.wmw <- apply(wmw.mat, 2, power)
  return(power.wmw)
}

cp.stud <- function(data, c, alpha)
{
  # run and calculate conditional power for Student's t test
  # data  ...null data with treatment effect applied via constte or expte
  # c     ...size of control group
  # alpha ...nominal size of the test
  
  t <- dim(data)[1] - c            
  N <- t + c                       
  cnm <- choose(N,c)               
  reps <- dim(data)[2]/cnm         
  
  studt <- apply(data, 2, fun_st, m = c, n = t)
  studt.mat <- matrix(studt, cnm, reps)
  power <- function(pvals)
  {
    length(which(pvals <= alpha + 1e-8))/cnm
  }    # 1e-8 is added to prevent rounding error
  power.studt <- apply(studt.mat, 2, power)
  return(power.studt)
}