Introduction

This experiment recruits participants to a chatbot survey using Facebook ads. After clicking on the Facebook ad, a participant is randomized into one of 5 treatment arms (the appendix considers 7 treatment arms).

The power analysis varies along two dimensions:

  • hypothesized outcomes for each treatment arm (and hence, treatment effects)
  • number of observations (clicks)

There are three outcomes, which are correlated: Pr(unvax complete|click), Pr(unvax website click|complete), and Pr(unvax useful free text|complete). We have an estimate of 14% for Pr(unvax complete|click) from the pilots and a rough estimate of 50% for Pr(unvax useful free text|complete), but we do not have an estimate for Pr(unvax website click|complete).

In the following power calculations, the hypothesized level of Pr(unvax complete|click) ranges from 0.04 to 0.20, with treatment effects ranging from 2 to 8 percentage points when comparing adjacent treatments (i.e. treatment 1 to 2, or 2 to 3). The hypothesized level of Pr(unvax useful free text|complete) ranges from 0.35 to 0.55, with treatment effects ranging from 2 to 10 percentage points when comparing adjacent treatments. The hypothesized level of Pr(unvax website click|complete) ranges from 0.05 to 0.35, with treatment effects ranging from 3 to 10 percentage points when comparing adjacent treatments. I do not anticipate that the Pr(unvax website click|complete) outcome will necessarily achieve a high of 0.35, but the outcome we end up using in its place (based on pilot results) should.

For the number of observations, we consider 50k to 150k clicks in 25k intervals.

Hypotheses

The pre-analysis plan specifies two research questions, with five hypotheses (one in two parts) to answer them. I list the research questions and accompanying hypotheses below.

  1. Is a chatbot a more effective tool for eliciting information than a traditional survey?
    • T1 = T4
  2. What quality of a chatbot makes it a more effective tool?
    • T1 = T2
    • T2 = T3
    • T3 = T4
    • T4 = T5
    • (T2 - T1) = (T4 - T5)

Note: The current code only includes the first four hypotheses for research question 2. It will be updated soon.

Key Takeaways

  • We are well-powered to answer research question 1 with 50k clicks for all outcomes, assuming a 2-3pp difference between adjacent treatment arms (so an 8 to 11pp difference between treatments 1 and 4).
  • For Pr(unvax complete|click), we are well-powered to detect 2pp differences between adjacent treatment arms at 50k clicks.
  • For Pr(unvax useful free text|complete), we are well-powered to detect 5pp differences between adjacent treatment arms at 75k clicks.
  • For Pr(unvax website click|complete), we are well-powered to detect 5pp differences between adjacent treatment arms at 50k clicks and 3pp differences at 150k clicks.

I recommend that the main experiment aim for 75k clicks at an approximate cost of $15k at $0.146 per click (plus mobile airtime costs which would amount to no more than $1500 at the most optimistic completion rate of 20%).

Power Calculations

To assess the power for testing multiple correlated hypotheses, I simulate data, test the hypotheses, correct the p-values, and repeat for r repetitions, then calculate the proportion of repetitions in which I reject each hypothesis. The variables needed in the simulated data are the treatment assignment and the three outcomes.

The data generation process is as follows:

  • Assign each observation to a treatment group. (treatment variable)
  • Take a draw from a binomial distribution with probability equal to the hypothesized unvax completion probability for the observation’s treatment group. (completion variable)
  • For observations who completed the survey, take a draw from a binomial distribution with probability equal to the hypothesized unvax website click probability for the observation’s treatment group. (webclick variable)
    • For observations who completed the survey, take a draw from a binomial distribution with probability equal to the hypothesized free text proportion for the observation’s treatment group. (free text variable)

Power is calculated using the following parameters:

  • \(\alpha=0.05\)
  • replications \(r=1000\)
  • number of observations (clicks) \(n=\{50k, 75k,100k,125k,150k\}\) for a cost of \(\{11k,15k,22k,29k\}\) at $0.146 per click.

Results are reported for uncorrected p-values and for p-values corrected using Bonferonni, Holm, and Hochberg.

Hypothesis

Hypothesis 1. Participants in Treatment 4 provide more information, higher-quality information, and are more receptive to new information about the vaccine compared to participants in Treatment 1.

Hypothesis 2. Participants in Treatment 3 provide more information, higher-quality information, and are more receptive to new information about the vaccine compared to participants in Treatment 2.

Hypothesis 3. Participants in Treatment 4 provide more information, higher-quality information, and are more receptive to new information about the vaccine compared to participants in Treatment 3.

Hypothesis 4. Participants in Treatment 4 provide more information, higher-quality information, and are more receptive to new information about the vaccine compared to participants in Treatment 5.

Hypothesis 5. (A) Participants in Treatment 2 provide more information, higher-quality information, and are more receptive to new information about the vaccine compared to participants in Treatment 1. (B) Moreover, the difference between participants in Treatment 2 and participants in Treatment 1 is greater than the difference between participants in Treatment 4 and participants in Treatment 5. –>

set.seed(94305)

### DEFINE PARAMETERS ###

n <- seq(50000, 150000, 25000)

# outcomes
outcomes <- matrix(c(   0.08,   0.12,   0.14,   0.16,   0.14, # Pr(unvax complete|click)
                 0.14,  0.19,   0.22,   0.25,   0.22, # Pr(unvax web click|complete)
                 0.44,  0.48,   0.5,    0.52,   0.5), # Pr(useful free text| complete)
                  nrow = 3, byrow = T)


outcomes2 <-  matrix(c( 0.04,   0.12,   0.16,   0.2,    0.16, # Pr(unvax complete|click)
                 0.05,  0.15,   0.2,    0.25,   0.2, # Pr(unvax web click|complete)
                 0.35,  0.45,   0.5,    0.55,   0.5), # Pr(useful free text| complete)
                  nrow = 3, byrow = T)

# example second outcomes matrix
outcomes3 <-  matrix(c( 0.04,   0.12,   0.16,   0.2,    0.14, # Pr(unvax complete|click)
                 0.22,  0.29,   0.32,   0.35,   0.32,
                    0.3,    0.4,    0.45,   0.55,   0.5), # Pr(useful free text| complete)
                  nrow = 3, byrow = T)


outcomes_list <- list(outcomes, outcomes2, outcomes3)

# treatments
t_list <- lapply(outcomes_list, FUN=function(x){c(1:ncol(x))})

# alpha
alpha <- 0.05

# number of replications
r <- 1000

### DEFINE FUNCTION ###

trial <- function(N, outcomes, t){
  
  ### Simulate Data
  
  # create dataframe
  d <- data.frame(matrix(nrow = N, ncol = 0))
  
  # treatment assignment
  d$treat <- sample(t, N, replace=T, prob = c(3/14,3/14,2/14,4/14,2/14)) 

  
  # generate data with completion outcome
  d <- d %>% 
    mutate(complete = rbinom(n(),1,outcomes[1,treat])) 
  
  d_complete <- d %>% 
    group_by(treat) %>% 
    summarise(
      complete = sum(complete),
      clicks = n()
    )
  
  # generate data with website click outcome
  d_website <- d %>% 
    filter(complete == 1) %>% 
    mutate(webclick = rbinom(n(),1,outcomes[2,treat])) %>% 
    group_by(treat) %>% 
    summarise(
      webclick = sum(webclick),
      completes = n()
    )
  
  # generate data with useful information outcome
  d_useful <- d %>% 
    filter(complete == 1) %>% 
    mutate(useful = rbinom(n(),1,outcomes[3,treat]))  %>% 
    group_by(treat) %>% 
    summarise(
      useful = sum(useful), 
      completes = n()
    )
  
  
  ### Calculate p-values
  
  # create matrix for p-values
  p_values <- matrix(nrow = nrow(outcomes),ncol = length(t) * (length(t) - 1)/2)
  
  # calculate for completions outcome
  d_subset <- cbind(d_complete[,2], d_complete[,3])
  p_values[1,] <- na.omit(as.numeric(pairwise.prop.test(d_subset[,1], d_subset[,2], 
                                                                  p.adj = "none")$p.value))
  
  # calculate for completions outcome
  d_subset <- cbind(d_website[,2], d_website[,3])
  p_values[2,] <- na.omit(as.numeric(pairwise.prop.test(d_subset[,1], d_subset[,2], 
                                                                  p.adj = "none")$p.value))
  
   # calculate for useful outcome
  d_subset <- cbind(d_useful[,2], d_useful[,3])
  p_values[3,] <- na.omit(as.numeric(pairwise.prop.test(d_subset[,1], d_subset[,2], 
                                                                  p.adj = "none")$p.value))
  
  ### MHT corrections
  
  return(rbind(p.adjust(p_values[1,], method = "none"),
               p.adjust(p_values[1,], method = "bonferroni"),
               p.adjust(p_values[1,], method = "holm"),
               p.adjust(p_values[1,], method = "hochberg"),
               p.adjust(p_values[2,], method = "none"),
               p.adjust(p_values[2,], method = "bonferroni"),
               p.adjust(p_values[2,], method = "holm"),
               p.adjust(p_values[2,], method = "hochberg"),
               p.adjust(p_values[3,], method = "none"),
               p.adjust(p_values[3,], method = "bonferroni"),
               p.adjust(p_values[3,], method = "holm"),
               p.adjust(p_values[3,], method = "hochberg")))
}


### POWER CALCULATIONS ###

for (i in n){
  for(o in 1:length(outcomes_list)){
    treatments <- t_list[[o]]
    outcome_matrix  <- outcomes_list[[o]]
    
    # complete r replications of n trials
    results <- replicate(r,  trial(i, outcomes = outcome_matrix, t = treatments))
    
    # create matrix to store power calculations
    power_calc <- matrix(NA, 
                         nrow = nrow(outcome_matrix)*4, 
                         ncol = (length(treatments) * (length(treatments) - 1)/2)) 
    
    # for each pairwise treatment comparison
    for(j in 1:(length(treatments) * (length(treatments) - 1)/2)){
      
      # for each mht correction
      for(k in 1:4){
        
        # calculate power for completions outcome
        power_calc[k, j] <- sum(results[k, j, ] < alpha)/r
        
        # calculate power for website clicks outcome
        power_calc[4+k,j] <- sum(results[k+4,j,]<alpha)/r
        
        # calculate power for useful free text
        power_calc[8+k,j] <- sum(results[k+8,j,]<alpha)/r
        
      }
    }
    
    power_calc <- data.frame(outcome = c(rep('Unvax completes|Click',4),
                           rep('Unvax web clicks|Complete', 4),
                           rep('Unvax useful freetext|Complete', 4)), 
                         test = rep(c("None","Bonferonni", "Holm", "Hochberg"), 3),
                         power_calc)
   
    colnames(power_calc) <- c("Outcome","Correction",  'T1 v. T2', 'T1 v. T3', 'T1 v. T4','T1 v. T5','T2 v. T3','T2 v. T4','T2 v. T5','T3 v. T4','T3 v. T5','T4 v. T5')
    
    kable(power_calc, format = "html",  caption = paste("Power Calculations: ",i/1000," Thousand Clicks", "Matrix", o), digits = 3) %>%
  kable_styling()%>%
  collapse_rows(columns = 1)%>%
      scroll_box(height = "100%")%>%
      print()
  }  
}
set.seed(94305)

### DEFINE PARAMETERS ###
# parameters defined in previous code chunk

### DEFINE FUNCTION ###

trial_subset <- function(N, outcomes, t){
  
  ### Simulate Data
  
  # create dataframe
  d <- data.frame(matrix(nrow = N, ncol = 0))
  
  # treatment assignment
  d$treat <- sample(t, N, replace=T) 
  
  # generate data with completion outcome
  d <- d %>% 
    mutate(complete = rbinom(n(),1,outcomes[1,treat])) 
  
  d_complete <- d %>% 
    group_by(treat) %>% 
    summarise(
      complete = sum(complete),
      clicks = n()
    )
  
  # generate data with website click outcome
  d_website <- d %>% 
    filter(complete == 1) %>% 
    mutate(webclick = rbinom(n(),1,outcomes[2,treat])) %>% 
    group_by(treat) %>% 
    summarise(
      webclick = sum(webclick),
      completes = n()
    )

    # generate data with useful information outcome
  d_useful <- d %>% 
    filter(complete == 1) %>% 
    mutate(useful = rbinom(n(),1,outcomes[3,treat]))  %>% 
    group_by(treat) %>% 
    summarise(
      useful = sum(useful), 
      completes = n()
    )
  
  
  ### Calculate p-values
  
  # create matrix for p-values
  p_values <- matrix(nrow = nrow(outcomes),ncol = length(t) - 1)
  
  # calculate for completions outcome
  d_subset <- cbind(d_complete[,2], d_complete[,3])
  p <- pairwise.prop.test(d_subset[,1], d_subset[,2], 
                                      p.adj = "none")$p.value
  p_values[1,] <-  diag(p)
  
  # calculate for webclicks outcome
  d_subset <- cbind(d_website[,2], d_website[,3])
  p <- pairwise.prop.test(d_subset[,1], d_subset[,2], 
                                      p.adj = "none")$p.value
  p_values[2,] <- diag(p)
  
  # calculate for useful information outcome
  d_subset <- cbind(d_useful[,2], d_useful[,3])
  p <- pairwise.prop.test(d_subset[,1], d_subset[,2], 
                                      p.adj = "none")$p.value
  p_values[3,] <- diag(p)
  
  
  ### MHT corrections
  
  return(rbind(p.adjust(p_values[1,], method = "none"),
               p.adjust(p_values[1,], method = "bonferroni"),
               p.adjust(p_values[1,], method = "holm"),
               p.adjust(p_values[1,], method = "hochberg"),
               p.adjust(p_values[2,], method = "none"),
               p.adjust(p_values[2,], method = "bonferroni"),
               p.adjust(p_values[2,], method = "holm"),
               p.adjust(p_values[2,], method = "hochberg"),
               p.adjust(p_values[3,], method = "none"),
               p.adjust(p_values[3,], method = "bonferroni"),
               p.adjust(p_values[3,], method = "holm"),
               p.adjust(p_values[3,], method = "hochberg")))
}

### POWER CALCULATIONS ###

for (i in n){
  for(o in 1:length(outcomes_list)){
    treatments <- t_list[[o]]
    outcome_matrix  <- outcomes_list[[o]]
    
    # complete r replications of n trials
    results <- replicate(r,  trial_subset(i, outcomes = outcome_matrix, t = treatments))
    
    # create matrix to store power calculations
    power_calc <- matrix(NA, 
                         nrow = nrow(outcome_matrix)*4, 
                         ncol = length(treatments)-1) 
    
    # for each pairwise treatment comparison
    # for each pairwise treatment comparison
    for(j in 1:(length(treatments) - 1)){
      
      # for each mht correction
      for(k in 1:4){
        
        # calculate power for completions outcome
        power_calc[k, j] <- sum(results[k, j, ] < alpha)/r
        
        # calculate power for website clicks outcome
        power_calc[4+k, j] <- sum(results[k+4, j, ] < alpha)/r
        
        # calculate power for useful informaton  outcome
        power_calc[8+k, j] <- sum(results[k+8, j, ] < alpha)/r
        
      }
    }
    
    # Generate Graph
  
    
        power_calc <- data.frame(outcome = c(rep('Unvax completes|Click',4),
                           rep('Unvax web clicks|Complete', 4),
                           rep('Unvax useful freetext|Complete', 4)), 
                         test = rep(c("None","Bonferonni", "Holm", "Hochberg"), 3),
                         power_calc)
   
    colnames(power_calc) <- c("Outcome","Correction", 
                              'T1 v. T2', 'T2 v. T3','T3 v. T4','T4 v. T5')
 
      
    kable(power_calc,  format = "html", caption = paste("Power Calculations: ",i/1000," Thousand Clicks", "Matrix", o), digits = 3)%>%
  kable_styling()%>%
  collapse_rows(columns = 1)%>%
      scroll_box(height = "100%")%>%
      print()
  }  
}
set.seed(94305)
### DEFINE PARAMETERS ###

#n <- c(40000)
#n <- c(seq(10000, 40000, 10000), seq(50000, 150000, 25000))

n<- c(seq(50000, 150000, 25000))
# outcomes
outcomes <- matrix(c(   0.08,   0.12,   0.14,   0.16,   0.14, # Pr(unvax complete|click)
                 0.14,  0.19,   0.22,   0.25,   0.22, # Pr(unvax web click|complete)
                 0.44,  0.48,   0.5,    0.52,   0.5), # Pr(useful free text| complete)
                  nrow = 3, byrow = T)


outcomes2 <-  matrix(c( 0.04,   0.12,   0.16,   0.2,    0.16, # Pr(unvax complete|click)
                 0.05,  0.15,   0.2,    0.25,   0.2, # Pr(unvax web click|complete)
                 0.35,  0.45,   0.5,    0.55,   0.5), # Pr(useful free text| complete)
                  nrow = 3, byrow = T)

outcomes3 <-  matrix(c( 0.04,   0.12,   0.16,   0.2,    0.14, # Pr(unvax complete|click)
                 0.22,  0.29,   0.32,   0.35,   0.32,
                    0.3,    0.4,    0.45,   0.55,   0.5), # Pr(useful free text| complete)
                  nrow = 3, byrow = T)

outcomes_list <- list(outcomes, outcomes2, outcomes3)

# treatments
t_list <- lapply(outcomes_list, FUN=function(x){c(1:ncol(x))})

# alpha
alpha <- 0.05

# number of replications
r <- 1000

### DEFINE FUNCTION ###

trial_subset <- function(N, outcomes, t){
  
  ### Simulate Data
  
  # create dataframe
  d <- data.frame(matrix(nrow = N, ncol = 0))
  
  # treatment assignment
  d$treat <- sample(t, N, replace=T, prob = c(3/14, 3/14, 2/14, 4/14, 2/14)) 
  
  # generate data with completion outcome
  d <- d %>% 
    mutate(complete = rbinom(n(),1,outcomes[1,treat])) 
  
  ## Outcome 1: completion
  
  # make df with columns for treatment and completion
  # this is for diff in diff pvalues
  d_complete <- d 
  
  # aggregate completion counts and click counts by treatment
  # this is for the pairwise prop tests
  d_complete_agg <- d_complete %>% 
    group_by(treat) %>% 
    summarise(
      complete = sum(complete),
      clicks = n()
    )

  ## Outcome 2: webclicks
  
  # make df with columns for treatment and webclick
  # conditional on completes
  # this is for diff in diff pvalues
  d_website <- d %>% 
    filter(complete == 1) %>% 
    mutate(webclick = rbinom(n(),1,outcomes[2,treat]))
  
  # aggregate webclick counts and completes counts by treatment
  # this is for the pairwise prop tests  
  d_website_agg <- d_website %>% 
    group_by(treat) %>% 
    summarise(
      webclick = sum(webclick),
      completes = n()
    )

  ## Outcome 3: Useful free text
  
  # make df with columns for treatment and useful
  # conditional on completes
  # this is for diff in diff pvalues
  d_useful <- d %>% 
    filter(complete == 1) %>% 
    mutate(useful = rbinom(n(),1,outcomes[3,treat])) 
  
  # aggregate webclick counts and completes counts by treatment
  # this is for the pairwise prop tests  
  d_useful_agg <- d_useful %>% 
    group_by(treat) %>% 
    summarise(
      useful = sum(useful), 
      completes = n()
    )
  
  ### Calculate p-values
  
  # create matrix for p-values
  # number of outcomes (3) by number of hypothesis (6)
  p_values <- matrix(nrow = nrow(outcomes),ncol = 6)
  
  
  # calculate for completions outcome
  
  # first get pvalues for usual comparisons (tA vs tB)
  d_subset <- cbind(d_complete_agg[,2], d_complete_agg[,3])
  p <- pairwise.prop.test(d_subset[,1], d_subset[,2], 
                                      p.adj = "none")$p.value
  p_values[1, 1:5] <- c(p["4", "1"], 
                    p["3", "2"], 
                    p["4", "3"], 
                    p["5", "4"], 
                    p["2", "1"])
  
  # now get pvalue for diff in diff
  model = lm(complete~0+treat, d_complete %>% 
               # remove observations with treatments not relevant to diff in diff
               filter(treat != 3) %>% 
               # make numeric categories a character so lm doesn't treat them as continuous
               mutate(treat = paste("W", treat, sep = ""))) 
  
  # run hypothesis test for diff in diff
  hyp = glht(model, "(treatW2 - treatW1) - (treatW5 - treatW4) = 0")
  
  # extract pvalue
  pval <- summary(hyp)$test$pvalues
  
  p_values[1, 6] <- c(pval)
  
  # calculate for webclicks outcome
  d_subset <- cbind(d_website_agg[,2], d_website_agg[,3])
  p <- pairwise.prop.test(d_subset[,1], d_subset[,2], 
                                      p.adj = "none")$p.value
  p_values[2, 1:5] <- c(p["4", "1"], 
                    p["3", "2"], 
                    p["4", "3"], 
                    p["5", "4"], 
                    p["2", "1"])
  
  ## diff in diff
  model = lm(webclick~0+treat, d_website %>% filter(treat != 3) %>% mutate(treat = paste("W", treat, sep = "")))

  hyp = glht(model, "(treatW2 - treatW1) - (treatW5 - treatW4) = 0")
  pval <- summary(hyp)$test$pvalues
  
  p_values[2, 6] <- c(pval)
  
  
  # calculate for useful information outcome
  d_subset <- cbind(d_useful_agg[,2], d_useful_agg[,3])
  p <- pairwise.prop.test(d_subset[,1], d_subset[,2], 
                                      p.adj = "none")$p.value
  
  p_values[3, 1:5] <- c(p["4", "1"], 
                    p["3", "2"], 
                    p["4", "3"], 
                    p["5", "4"], 
                    p["2", "1"])
  
    ## diff in diff
  model = lm(useful~0+treat, d_useful %>% filter(treat != 3) %>% mutate(treat = paste("W", treat, sep = "")))

  hyp = glht(model, "(treatW2 - treatW1) - (treatW5 - treatW4) = 0")
  pval <- summary(hyp)$test$pvalues
  
  p_values[3, 6] <- c(pval)
  ### MHT corrections
  
  return(rbind(p.adjust(p_values[1,], method = "none"),
               p.adjust(p_values[1,], method = "bonferroni"),
               p.adjust(p_values[1,], method = "holm"),
               p.adjust(p_values[1,], method = "hochberg"),
               p.adjust(p_values[2,], method = "none"),
               p.adjust(p_values[2,], method = "bonferroni"),
               p.adjust(p_values[2,], method = "holm"),
               p.adjust(p_values[2,], method = "hochberg"),
               p.adjust(p_values[3,], method = "none"),
               p.adjust(p_values[3,], method = "bonferroni"),
               p.adjust(p_values[3,], method = "holm"),
               p.adjust(p_values[3,], method = "hochberg")))
}

### POWER CALCULATIONS ###

for (i in n){
  
    cat('###',i,' \n')

  
  for(o in 1:length(outcomes_list)){
    treatments <- t_list[[o]]
    outcome_matrix  <- outcomes_list[[o]]
    
    # complete r replications of n trials
    results <- replicate(r,  trial_subset(i, outcomes = outcome_matrix, t = treatments))
    
    # create matrix to store power calculations
    power_calc <- matrix(NA, 
                         nrow = nrow(outcome_matrix)*4, 
                         ncol = 6) 
    
    # for each pairwise treatment comparison
    for(j in 1:6){
      
      # for each mht correction
      for(k in 1:4){
        
        # calculate power for completions outcome
        power_calc[k, j] <- sum(results[k, j, ] < alpha)/r
        
        # calculate power for website clicks outcome
        power_calc[4+k, j] <- sum(results[k+4, j, ] < alpha)/r
        
        # calculate power for useful informaton  outcome
        power_calc[8+k, j] <- sum(results[k+8, j, ] < alpha)/r
        
      }
    }
    
    # Generate Graph
  
    
        power_calc <- data.frame(outcome = c(rep('Unvax completes|Click',4),
                           rep('Unvax web clicks|Complete', 4),
                           rep('Unvax useful freetext|Complete', 4)), 
                         test = rep(c("None","Bonferonni", "Holm", "Hochberg"), 3),
                         power_calc)
   
    colnames(power_calc) <- c("Outcome","Correction", 
                              'T1 v. T4', 'T2 v. T3','T3 v. T4','T4 v. T5', 'T1 v. T2', 'T2 - T1 v. T5 - T4')
 
      
    kable(power_calc,  format = "html", caption = paste("Power Calculations: ",i/1000," Thousand Clicks", "Matrix", o), digits = 3)%>%
  kable_styling()%>%
  collapse_rows(columns = 1)%>%
      scroll_box(height = "100%")%>%
      print()
  }  
}

50000

Power Calculations: 50 Thousand Clicks Matrix 1
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1.000 0.972 0.966 0.961 1.000 1.000
Bonferonni 1.000 0.898 0.885 0.870 1.000 1.000
Holm 1.000 0.971 0.962 0.955 1.000 1.000
Hochberg 1.000 0.971 0.964 0.957 1.000 1.000
Unvax web clicks|Complete None 1.000 0.393 0.470 0.444 0.852 0.921
Bonferonni 1.000 0.180 0.214 0.209 0.642 0.793
Holm 1.000 0.251 0.307 0.301 0.718 0.821
Hochberg 1.000 0.257 0.318 0.309 0.726 0.825
Unvax useful freetext|Complete None 0.978 0.142 0.176 0.154 0.435 0.543
Bonferonni 0.914 0.037 0.048 0.053 0.192 0.284
Holm 0.914 0.043 0.061 0.067 0.229 0.317
Hochberg 0.914 0.045 0.062 0.070 0.233 0.321
Power Calculations: 50 Thousand Clicks Matrix 2
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1.000 1.000
Bonferonni 1 1.000 1.000 1.000 1.000 1.000
Holm 1 1.000 1.000 1.000 1.000 1.000
Hochberg 1 1.000 1.000 1.000 1.000 1.000
Unvax web clicks|Complete None 1 0.907 0.927 0.921 1.000 1.000
Bonferonni 1 0.744 0.744 0.752 1.000 1.000
Holm 1 0.896 0.917 0.901 1.000 1.000
Hochberg 1 0.901 0.920 0.909 1.000 1.000
Unvax useful freetext|Complete None 1 0.675 0.801 0.808 0.943 0.997
Bonferonni 1 0.393 0.579 0.583 0.838 0.975
Holm 1 0.611 0.734 0.732 0.919 0.984
Hochberg 1 0.631 0.750 0.751 0.925 0.985
Power Calculations: 50 Thousand Clicks Matrix 3
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1.000 1.000
Bonferonni 1 1.000 1.000 1.000 1.000 1.000
Holm 1 1.000 1.000 1.000 1.000 1.000
Hochberg 1 1.000 1.000 1.000 1.000 1.000
Unvax web clicks|Complete None 1 0.351 0.452 0.398 0.815 0.912
Bonferonni 1 0.132 0.202 0.173 0.550 0.733
Holm 1 0.194 0.276 0.258 0.642 0.767
Hochberg 1 0.210 0.296 0.280 0.661 0.771
Unvax useful freetext|Complete None 1 0.677 1.000 0.772 0.952 0.998
Bonferonni 1 0.405 0.999 0.518 0.848 0.975
Holm 1 0.640 1.000 0.724 0.938 0.988
Hochberg 1 0.648 1.000 0.735 0.945 0.988

75000

Power Calculations: 75 Thousand Clicks Matrix 1
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1.000 0.998 0.995 0.997 1.000 1.000
Bonferonni 1.000 0.984 0.972 0.980 1.000 1.000
Holm 1.000 0.998 0.995 0.997 1.000 1.000
Hochberg 1.000 0.998 0.995 0.997 1.000 1.000
Unvax web clicks|Complete None 1.000 0.557 0.624 0.636 0.954 0.990
Bonferonni 1.000 0.301 0.355 0.364 0.857 0.945
Holm 1.000 0.449 0.505 0.511 0.920 0.964
Hochberg 1.000 0.482 0.543 0.538 0.929 0.966
Unvax useful freetext|Complete None 0.997 0.213 0.233 0.265 0.599 0.720
Bonferonni 0.984 0.060 0.085 0.100 0.315 0.452
Holm 0.985 0.079 0.108 0.131 0.356 0.500
Hochberg 0.985 0.086 0.112 0.137 0.359 0.504
Power Calculations: 75 Thousand Clicks Matrix 2
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1.000 1.000
Bonferonni 1 1.000 1.000 1.000 1.000 1.000
Holm 1 1.000 1.000 1.000 1.000 1.000
Hochberg 1 1.000 1.000 1.000 1.000 1.000
Unvax web clicks|Complete None 1 0.969 0.984 0.985 1.000 1.000
Bonferonni 1 0.904 0.936 0.941 1.000 1.000
Holm 1 0.969 0.984 0.985 1.000 1.000
Hochberg 1 0.969 0.984 0.985 1.000 1.000
Unvax useful freetext|Complete None 1 0.833 0.954 0.947 0.996 1.000
Bonferonni 1 0.627 0.823 0.794 0.962 0.998
Holm 1 0.824 0.945 0.932 0.993 0.999
Hochberg 1 0.826 0.948 0.937 0.993 0.999
Power Calculations: 75 Thousand Clicks Matrix 3
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1.000 1.000
Bonferonni 1 1.000 1.000 1.000 1.000 1.000
Holm 1 1.000 1.000 1.000 1.000 1.000
Hochberg 1 1.000 1.000 1.000 1.000 1.000
Unvax web clicks|Complete None 1 0.480 0.589 0.540 0.935 0.980
Bonferonni 1 0.217 0.325 0.286 0.774 0.914
Holm 1 0.343 0.451 0.408 0.854 0.941
Hochberg 1 0.366 0.478 0.430 0.867 0.942
Unvax useful freetext|Complete None 1 0.865 1.000 0.919 0.995 1.000
Bonferonni 1 0.657 1.000 0.768 0.972 0.996
Holm 1 0.860 1.000 0.907 0.994 0.999
Hochberg 1 0.863 1.000 0.910 0.994 0.999

100000

Power Calculations: 100 Thousand Clicks Matrix 1
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1.000 1.000 0.998 1.000 1.000 1.000
Bonferonni 1.000 0.997 0.997 0.997 1.000 1.000
Holm 1.000 1.000 0.998 1.000 1.000 1.000
Hochberg 1.000 1.000 0.998 1.000 1.000 1.000
Unvax web clicks|Complete None 1.000 0.669 0.755 0.744 0.990 0.997
Bonferonni 1.000 0.415 0.500 0.488 0.947 0.987
Holm 1.000 0.605 0.677 0.663 0.976 0.992
Hochberg 1.000 0.628 0.699 0.690 0.983 0.993
Unvax useful freetext|Complete None 0.999 0.265 0.305 0.299 0.732 0.828
Bonferonni 0.999 0.093 0.121 0.119 0.476 0.617
Holm 0.999 0.118 0.164 0.167 0.543 0.653
Hochberg 0.999 0.120 0.168 0.169 0.544 0.654
Power Calculations: 100 Thousand Clicks Matrix 2
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1.000 1
Bonferonni 1 1.000 1.000 1.000 1.000 1
Holm 1 1.000 1.000 1.000 1.000 1
Hochberg 1 1.000 1.000 1.000 1.000 1
Unvax web clicks|Complete None 1 0.993 0.997 0.999 1.000 1
Bonferonni 1 0.977 0.984 0.985 1.000 1
Holm 1 0.992 0.997 0.998 1.000 1
Hochberg 1 0.993 0.997 0.999 1.000 1
Unvax useful freetext|Complete None 1 0.931 0.978 0.977 0.999 1
Bonferonni 1 0.795 0.921 0.917 0.994 1
Holm 1 0.929 0.977 0.975 0.999 1
Hochberg 1 0.929 0.978 0.976 0.999 1
Power Calculations: 100 Thousand Clicks Matrix 3
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1.000 1.000
Bonferonni 1 1.000 1.000 1.000 1.000 1.000
Holm 1 1.000 1.000 1.000 1.000 1.000
Hochberg 1 1.000 1.000 1.000 1.000 1.000
Unvax web clicks|Complete None 1 0.591 0.706 0.657 0.978 0.995
Bonferonni 1 0.344 0.449 0.394 0.909 0.971
Holm 1 0.501 0.616 0.564 0.956 0.982
Hochberg 1 0.526 0.637 0.586 0.964 0.983
Unvax useful freetext|Complete None 1 0.946 1.000 0.966 1.000 1.000
Bonferonni 1 0.830 1.000 0.878 0.995 1.000
Holm 1 0.944 1.000 0.962 0.999 1.000
Hochberg 1 0.946 1.000 0.965 1.000 1.000

125000

Power Calculations: 125 Thousand Clicks Matrix 1
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1.000 1.000
Bonferonni 1 1.000 1.000 1.000 1.000 1.000
Holm 1 1.000 1.000 1.000 1.000 1.000
Hochberg 1 1.000 1.000 1.000 1.000 1.000
Unvax web clicks|Complete None 1 0.805 0.832 0.820 0.999 1.000
Bonferonni 1 0.562 0.611 0.607 0.988 0.998
Holm 1 0.766 0.799 0.776 0.998 1.000
Hochberg 1 0.775 0.812 0.791 0.998 1.000
Unvax useful freetext|Complete None 1 0.302 0.380 0.369 0.820 0.905
Bonferonni 1 0.123 0.174 0.157 0.600 0.722
Holm 1 0.154 0.238 0.232 0.670 0.772
Hochberg 1 0.159 0.248 0.241 0.677 0.773
Power Calculations: 125 Thousand Clicks Matrix 2
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1 1
Bonferonni 1 1.000 1.000 1.000 1 1
Holm 1 1.000 1.000 1.000 1 1
Hochberg 1 1.000 1.000 1.000 1 1
Unvax web clicks|Complete None 1 1.000 1.000 0.999 1 1
Bonferonni 1 0.995 0.996 0.995 1 1
Holm 1 1.000 1.000 0.999 1 1
Hochberg 1 1.000 1.000 0.999 1 1
Unvax useful freetext|Complete None 1 0.982 0.997 0.990 1 1
Bonferonni 1 0.903 0.962 0.957 1 1
Holm 1 0.982 0.997 0.990 1 1
Hochberg 1 0.982 0.997 0.990 1 1
Power Calculations: 125 Thousand Clicks Matrix 3
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1.000 1.000
Bonferonni 1 1.000 1.000 1.000 1.000 1.000
Holm 1 1.000 1.000 1.000 1.000 1.000
Hochberg 1 1.000 1.000 1.000 1.000 1.000
Unvax web clicks|Complete None 1 0.668 0.824 0.762 0.997 0.999
Bonferonni 1 0.430 0.595 0.514 0.970 0.995
Holm 1 0.612 0.762 0.693 0.995 0.998
Hochberg 1 0.628 0.776 0.712 0.995 0.998
Unvax useful freetext|Complete None 1 0.976 1.000 0.989 1.000 1.000
Bonferonni 1 0.893 1.000 0.942 1.000 1.000
Holm 1 0.975 1.000 0.988 1.000 1.000
Hochberg 1 0.976 1.000 0.989 1.000 1.000

150000

Power Calculations: 150 Thousand Clicks Matrix 1
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1.000 1.000
Bonferonni 1 1.000 1.000 1.000 1.000 1.000
Holm 1 1.000 1.000 1.000 1.000 1.000
Hochberg 1 1.000 1.000 1.000 1.000 1.000
Unvax web clicks|Complete None 1 0.847 0.892 0.896 1.000 1.000
Bonferonni 1 0.657 0.702 0.724 0.996 1.000
Holm 1 0.830 0.871 0.869 1.000 1.000
Hochberg 1 0.836 0.875 0.878 1.000 1.000
Unvax useful freetext|Complete None 1 0.346 0.449 0.438 0.879 0.950
Bonferonni 1 0.134 0.201 0.201 0.692 0.837
Holm 1 0.209 0.282 0.287 0.759 0.864
Hochberg 1 0.224 0.292 0.297 0.771 0.865
Power Calculations: 150 Thousand Clicks Matrix 2
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1 1
Bonferonni 1 1.000 1.000 1.000 1 1
Holm 1 1.000 1.000 1.000 1 1
Hochberg 1 1.000 1.000 1.000 1 1
Unvax web clicks|Complete None 1 1.000 1.000 1.000 1 1
Bonferonni 1 1.000 1.000 1.000 1 1
Holm 1 1.000 1.000 1.000 1 1
Hochberg 1 1.000 1.000 1.000 1 1
Unvax useful freetext|Complete None 1 0.991 1.000 0.997 1 1
Bonferonni 1 0.950 0.994 0.990 1 1
Holm 1 0.991 1.000 0.997 1 1
Hochberg 1 0.991 1.000 0.997 1 1
Power Calculations: 150 Thousand Clicks Matrix 3
Outcome Correction T1 v. T4 T2 v. T3 T3 v. T4 T4 v. T5 T1 v. T2 T2 - T1 v. T5 - T4
Unvax completes|Click None 1 1.000 1.000 1.000 1.000 1.000
Bonferonni 1 1.000 1.000 1.000 1.000 1.000
Holm 1 1.000 1.000 1.000 1.000 1.000
Hochberg 1 1.000 1.000 1.000 1.000 1.000
Unvax web clicks|Complete None 1 0.817 0.880 0.836 0.999 1.000
Bonferonni 1 0.562 0.681 0.603 0.989 0.998
Holm 1 0.784 0.853 0.790 0.999 1.000
Hochberg 1 0.798 0.870 0.813 0.999 1.000
Unvax useful freetext|Complete None 1 0.992 1.000 0.998 1.000 1.000
Bonferonni 1 0.951 1.000 0.987 1.000 1.000
Holm 1 0.992 1.000 0.998 1.000 1.000
Hochberg 1 0.992 1.000 0.998 1.000 1.000

Appendix

Sample size calculations – Bonferonni – pr(unvax complete|click)

The below figure uses the below equation to calculate the sample size needed for each treatment arm, given the parameters.

  • The outcome is the probability of an unvaccinated person completing the survey, conditional on clicking on the Facebook ad.
  • The sample for this outcome is the number of clicks.
  • The baseline of this outcome from the pilot is ~14%.
  • I consider treatment effects ranging from 10% to 30% of baseline (i.e. ~1.4% to 4.2%)
  • I correct for multiple hypothesis testing using Bonferroni.
  • For the Bonferroni correction, I assume the number of hypotheses = number of outcomes (3) x number of treatment comparisons (factorial(t)/2) for 5 or 7 treatment arms.

\[ n \geq \dfrac{4(1-p)}{\lambda^2p} \times (z_\alpha + \Phi^{-1}(1-\beta))^2 \]

### DEFINE PARAMETERS ###

# outcome 
p <- 0.14365 #from pilot 5 descriptive analytics analysis 595/4142

# standard deviation
sd <- sqrt(p * (1 - p))

# define the possible treatment effect sizes in *percent*
lambda <- seq(from = 0.1, to = 0.3, length.out = 100 )

# alpha
alpha <- 0.05

# beta
beta <- 0.2

# number of treatments
treat <- cbind(5,7)

### DEFINE FUNCTION ###

vax_power <- function(t, lambda) {
  4 * (qnorm(1-alpha/2/(3*factorial(t)/2))+(qnorm(1 - beta)))^2 * ((1 - p)/((lambda * lambda) * p))}

### CREATE PLOT ###

# create vector to hold number of observations for each value of lambda and number of treatments
n_obs <- matrix(NA, nrow = length(lambda), ncol = length(treat))

# calculate number of observations for each value of lambda
for (i in 1:length(lambda)) {
  for (j in 1:length(treat)){
    n_obs[i,j] <- vax_power(treat[j], lambda[i])
    }
  }

# create data frame that combines data to graph
df <- data.frame(cbind(n_obs,lambda))

# create the graph
ggplot(df, aes(x = lambda)) +
  geom_line(aes(y = df[,1], color = "5")) +
  geom_line(aes(y = df[,2], color = "7")) +
  labs(x = 'Treatment effect size (in percentage terms)',
       y = 'Number of clicks needed per treatment arm') +
  scale_color_manual(name = "Number of treatment arms", values = c("5" = "darkblue", "7" = "red"))

Then, we can calculate total costs based on $0.155 per click.

Let’s say instead that we only compare each treatment arm to one other treatment arm. This approach is what we have outlined in the experimental design document. Then, the below figure contains the same information as the above figure, but we are conducting fewer hypothesis tests.

### DEFINE PARAMETERS ###

# outcome 
p <- 0.14365 #from pilot 5 descriptive analytics analysis 595/4142

# standard deviation
sd <- sqrt(p * (1 - p))

# define the possible treatment effect sizes in *percent*
lambda <- seq(from = 0.1, to = 0.3, length.out = 100 )

# alpha
alpha <- 0.05

# beta
beta <- 0.2

# number of treatments
treat <- cbind(5,7)

### DEFINE FUNCTION ###

vax_power <- function(t, lambda) {
  4 * (qnorm(1-alpha/2/(3*(t-1)))+(qnorm(1 - beta)))^2 * ((1 - p)/((lambda * lambda) * p))}

### CREATE PLOT ###

# create vector to hold number of observations for each value of lambda and number of treatments
n_obs <- matrix(NA, nrow = length(lambda), ncol = length(treat))

# calculate number of observations for each value of lambda
for (i in 1:length(lambda)) {
  for (j in 1:length(treat)){
    n_obs[i,j] <- vax_power(treat[j], lambda[i])
    }
  }

# create data frame that combines data to graph
df <- data.frame(cbind(n_obs,lambda))

# create the graph
ggplot(df, aes(x = lambda)) +
  geom_line(aes(y = df[,1], color = "5")) +
  geom_line(aes(y = df[,2], color = "7")) +
  labs(x = 'Treatment effect size (in percentage terms)',
       y = 'Number of clicks needed per treatment arm') +
  scale_color_manual(name = "Number of treatment arms", values = c("5" = "darkblue", "7" = "red"))

Then, we can again calculate total costs based on $0.155 per click.

print(paste("For treatment effect of 15% over baseline, the treatment effect would need to be",round(0.15*p,4)," and you would need:"))
## [1] "For treatment effect of 15% over baseline, the treatment effect would need to be 0.0215  and you would need:"
print(paste(round(vax_power(treat,0.15),0),"clicks per treatment arm for", treat,"treatment arms, for a total number of ", round(treat*vax_power(treat,0.15),0),"clicks at a total cost of $",round(0.155*treat*vax_power(treat,0.15),0)))
## [1] "14563 clicks per treatment arm for 5 treatment arms, for a total number of  72813 clicks at a total cost of $ 11286" 
## [2] "15570 clicks per treatment arm for 7 treatment arms, for a total number of  108989 clicks at a total cost of $ 16893"

Sample size calculations – Bonferonni – pr(unvax web click|complete)

The below figure uses the previous equation to calculate the sample size needed for each treatment arm, given the parameters.

  • The outcome is the probability of an unvaccinated person clicking the website link, conditional on completing the survey.
  • So, the sample for this outcome is the number of completes.
  • We have not piloted this outcome, so I am guessing the baseline to be ~20% (this is probably a bit high, but it is also the only way to get costs in a reasonable range – pilot will inform if this is reasonable).
  • I consider treatment effects ranging from 10% to 30% of baseline (i.e. ~1% to 4%)
  • I correct for multiple hypothesis testing using Bonferroni.
  • For the Bonferroni correction, I assume the number of hypotheses = number of outcomes (3) x number of treatment comparisons (factorial(t)/2) for 5 or 7 treatment arms.
### DEFINE PARAMETERS ###

# outcome 
p <- 0.2

# standard deviation
sd <- sqrt(p * (1 - p))

# define the possible treatment effect sizes in *percent*
lambda <- seq(from = 0.1, to = 0.3, length.out = 100 )

# alpha
alpha <- 0.05

# beta
beta <- 0.2

# number of treatments
treat <- cbind(5,7)

### DEFINE FUNCTION ###

vax_power <- function(t, lambda) {
  4 * (qnorm(1-alpha/2/(3*factorial(t)/2))+(qnorm(1 - beta)))^2 * ((1 - p)/((lambda * lambda) * p))}

### CREATE PLOT ###

# create vector to hold number of observations for each value of lambda and number of treatments
n_obs <- matrix(NA, nrow = length(lambda), ncol = length(treat))

# calculate number of observations for each value of lambda
for (i in 1:length(lambda)) {
  for (j in 1:length(treat)){
    n_obs[i,j] <- vax_power(treat[j], lambda[i])
    }
  }

# create data frame that combines data to graph
df <- data.frame(cbind(n_obs,lambda))

# create the graph
ggplot(df, aes(x = lambda)) +
  geom_line(aes(y = df[,1], color = "5")) +
  geom_line(aes(y = df[,2], color = "7")) +
  labs(x = 'Treatment effect size (in percentage terms)',
       y = 'Number of completes needed per treatment arm') +
  scale_color_manual(name = "Number of treatment arms", values = c("5" = "darkblue", "7" = "red"))

Then, we can calculate total costs based on $1.30 per unvaccinated complete.

print(paste("For treatment effect of 25% over baseline, the treatment effect would need to be",round(0.25*p,4)," and you would need:"))
## [1] "For treatment effect of 25% over baseline, the treatment effect would need to be 0.05  and you would need:"
print(paste(round(vax_power(treat,0.25),0),"completes per treatment arm for", treat,"treatment arms, for a total number of ", round(treat*vax_power(treat,0.25),0),"clicks at a total cost of $",round(1.3*treat*vax_power(treat,0.25),0)))
## [1] "5131 completes per treatment arm for 5 treatment arms, for a total number of  25654 clicks at a total cost of $ 33350"
## [2] "7320 completes per treatment arm for 7 treatment arms, for a total number of  51241 clicks at a total cost of $ 66614"

Let’s say instead that we only compare each treatment arm to one other treatment arm. This approach is what we have outlined in the experimental design document. Then, the below figure contains the same information as the above figure, but we are conducting fewer hypothesis tests.

### DEFINE PARAMETERS ###

# outcome 
p <- 0.2 

# standard deviation
sd <- sqrt(p * (1 - p))

# define the possible treatment effect sizes in *percent*
lambda <- seq(from = 0.1, to = 0.3, length.out = 100 )

# alpha
alpha <- 0.05

# beta
beta <- 0.2

# number of treatments
treat <- cbind(5,7)

### DEFINE FUNCTION ###

vax_power <- function(t, lambda) {
  4 * (qnorm(1-alpha/2/(3*(t-1)))+(qnorm(1 - beta)))^2 * ((1 - p)/((lambda * lambda) * p))}

### CREATE PLOT ###

# create vector to hold number of observations for each value of lambda and number of treatments
n_obs <- matrix(NA, nrow = length(lambda), ncol = length(treat))

# calculate number of observations for each value of lambda
for (i in 1:length(lambda)) {
  for (j in 1:length(treat)){
    n_obs[i,j] <- vax_power(treat[j], lambda[i])
    }
  }

# create data frame that combines data to graph
df <- data.frame(cbind(n_obs,lambda))

# create the graph
ggplot(df, aes(x = lambda)) +
  geom_line(aes(y = df[,1], color = "5")) +
  geom_line(aes(y = df[,2], color = "7")) +
  labs(x = 'Treatment effect size (in percentage terms)',
       y = 'Number of completes needed per treatment arm') +
  scale_color_manual(name = "Number of treatment arms", values = c("5" = "darkblue", "7" = "red"))

Then, we can again calculate total costs based on $1.30 per unvaccinated complete.

print(paste("For treatment effect of 25% over baseline, the treatment effect would need to be",round(0.25*p,4)," and you would need:"))
## [1] "For treatment effect of 25% over baseline, the treatment effect would need to be 0.05  and you would need:"
print(paste(round(vax_power(treat,0.25),0),"completes per treatment arm for", treat,"treatment arms, for a total number of ", round(treat*vax_power(treat,0.25),0),"completes at a total cost of $",round(1.3*treat*vax_power(treat,0.25),0)))
## [1] "3518 completes per treatment arm for 5 treatment arms, for a total number of  17588 completes at a total cost of $ 22865"
## [2] "3761 completes per treatment arm for 7 treatment arms, for a total number of  26327 completes at a total cost of $ 34225"

Power calculations – Bonferonni – pr(unvax complete|click)

The below figure uses the below equation to calculate the power needed for each treatment arm, given the parameters.

  • The outcome is the probability of an unvaccinated person completing the survey, conditional on clicking on the Facebook ad.
  • The sample for this outcome is the number of clicks.
  • The baseline of this outcome from the pilot is ~14%.
  • I consider treatment effects ranging from 10% to 20% of baseline (i.e. ~1.4% to 2.8%)
  • I correct for multiple hypothesis testing using Bonferroni.
  • For the Bonferroni correction, I assume the number of hypotheses = number of outcomes (3) x number of treatment comparisons (factorial(t)/2) for 5 or 7 treatment arms.
  • I used the number of observations that seemed in a reasonable range cost-wise from the sample size analysis to set the scale for the power analysis.

\[ \boxed{1-\beta = \Phi\left( \frac{\tau}{2\sigma/\sqrt{n}} - z_{1-\alpha} \right)} \]

## DEFINE PARAMETERS ###

# define the possible total number of clicks
n <- seq(from = 50000, to = 300000, length.out = 100 )

# outcome 
p <- 0.14365 #from pilot 5 descriptive analytics analysis 595/4142

# standard deviation
sd <- sqrt(p * (1 - p))

# define the possible treatment effect sizes in *percent*
lambda <- cbind(0.1, 0.15, 0.2)

# alpha
alpha <- 0.05

# beta
beta <- 0.2

# number of treatments
treat <- cbind(5,7)

### DEFINE FUNCTION ###

vax_power2 <- function(lambda, n, t) {
  pnorm(
    ( lambda*p / (2 * sd / sqrt(n / t))) - qnorm(1-alpha/2/(3*factorial(t)/2))
  )
}

### CREATE PLOT ###

# create vector to hold number of observations for each value of lambda, number of observations, and number of treatments
power <- matrix(NA, nrow = length(n), ncol = length(treat)*length(lambda))

# calculate number of observations for each value of lambda
for (i in 1:length(lambda)) {
  for (j in 1:length(treat)){
    for (k in 1:length(n)){
      power[k,j+2*i-2] <- vax_power2(lambda[i], n[k], treat[j])
      }
    }
  }

# create data frame that combines data to graph
df <- data.frame(cbind(power,n))

# create the graph
ggplot(df, aes(x = n)) +
  geom_line(aes(y = df[,1], color = "lambda = 0.1, treat = 5")) +
  geom_line(aes(y = df[,2], color = "lambda = 0.1, treat = 7")) +
  geom_line(aes(y = df[,3], color = "lambda = 0.15, treat = 5")) +
  geom_line(aes(y = df[,4], color = "lambda = 0.15, treat = 7")) +
  geom_line(aes(y = df[,5], color = "lambda = 0.2, treat = 5")) +
  geom_line(aes(y = df[,6], color = "lambda = 0.2, treat = 7")) +
  labs(x = 'Number of Total Clicks',
       y = 'Power') +
  scale_color_manual(name = "Number of treatment arms", values = c("lambda = 0.1, treat = 5" = "darkblue", "lambda = 0.1, treat = 7" = "red", "lambda = 0.15, treat = 5" = "darkgreen", "lambda = 0.15, treat = 7" = "turquoise4", "lambda = 0.2, treat = 5" = "darkorchid4", "lambda = 0.2, treat = 7" = "springgreen2"))

Let’s say instead that we only compare each treatment arm to one other treatment arm. This approach is what we have outlined in the experimental design document. Then, the below figure contains the same information as the above figure, but we are conducting fewer hypothesis tests.

## DEFINE PARAMETERS ###

# define the possible total number of clicks
n <- seq(from = 50000, to = 300000, length.out = 100 )

# outcome 
p <- 0.14365 #from pilot 5 descriptive analytics analysis 595/4142

# standard deviation
sd <- sqrt(p * (1 - p))

# define the possible treatment effect sizes in *percent*
lambda <- cbind(0.1, 0.15, 0.2)

# alpha
alpha <- 0.05

# beta
beta <- 0.2

# number of treatments
treat <- cbind(5,7)

### DEFINE FUNCTION ###

vax_power2 <- function(lambda, n, t) {
  pnorm(
    ( lambda*p / (2 * sd / sqrt(n / t))) - qnorm(1-alpha/2/(3*(t-1)))
  )
}

### CREATE PLOT ###

# create vector to hold number of observations for each value of lambda, number of observations, and number of treatments
power <- matrix(NA, nrow = length(n), ncol = length(treat)*length(lambda))

# calculate number of observations for each value of lambda
for (i in 1:length(lambda)) {
  for (j in 1:length(treat)){
    for (k in 1:length(n)){
      power[k,j+2*i-2] <- vax_power2(lambda[i], n[k], treat[j])
      }
    }
  }

# create data frame that combines data to graph
df <- data.frame(cbind(power,n))

# create the graph
ggplot(df, aes(x = n)) +
  geom_line(aes(y = df[,1], color = "lambda = 0.1, treat = 5")) +
  geom_line(aes(y = df[,2], color = "lambda = 0.1, treat = 7")) +
  geom_line(aes(y = df[,3], color = "lambda = 0.15, treat = 5")) +
  geom_line(aes(y = df[,4], color = "lambda = 0.15, treat = 7")) +
  geom_line(aes(y = df[,5], color = "lambda = 0.2, treat = 5")) +
  geom_line(aes(y = df[,6], color = "lambda = 0.2, treat = 7")) +
  labs(x = 'Number of Total Clicks',
       y = 'Power') +
  scale_color_manual(name = "Number of treatment arms", values = c("lambda = 0.1, treat = 5" = "darkblue", "lambda = 0.1, treat = 7" = "red", "lambda = 0.15, treat = 5" = "darkgreen", "lambda = 0.15, treat = 7" = "turquoise4", "lambda = 0.2, treat = 5" = "darkorchid4", "lambda = 0.2, treat = 7" = "springgreen2"))

Power calculations – Bonferonni – pr(unvax web click|complete)

The below figure uses the previous equation to calculate the power needed for each treatment arm, given the parameters.

  • The outcome is the probability of an unvaccinated person clicking the website link, conditional on completing the survey.
  • So, the sample for this outcome is the number of completes.
  • We have not piloted this outcome, so I am guessing the baseline to be ~20% (this is probably a bit high, but it is also the only way to get costs in a reasonable range – pilot will inform if this is reasonable).
  • I correct for multiple hypothesis testing using Bonferroni.
  • For the Bonferroni correction, I assume the number of hypotheses = number of outcomes (3) x number of treatment comparisons (factorial(t)/2) for 5 or 7 treatment arms.
  • I used the number of observations that seemed in a reasonable range cost-wise from the sample size analysis to set the scale for the power analysis.

\[ \boxed{1-\beta = \Phi\left( \frac{\tau}{2\sigma/\sqrt{n}} - z_{1-\alpha} \right)} \]

## DEFINE PARAMETERS ###

# define the possible total number of completes
n <- seq(from = 10000, to = 100000, length.out = 100 )

# outcome 
p <- 0.2 #from pilot 5 descriptive analytics analysis 595/4142

# standard deviation
sd <- sqrt(p * (1 - p))

# define the possible treatment effect sizes in *percent*
lambda <- cbind(0.2, 0.25, 0.3)

# alpha
alpha <- 0.05

# beta
beta <- 0.2

# number of treatments
treat <- cbind(5,7)

### DEFINE FUNCTION ###

vax_power2 <- function(lambda, n, t) {
  pnorm(
    ( lambda*p / (2 * sd / sqrt(n / t))) - qnorm(1-alpha/2/(3*factorial(t)/2))
  )
}

### CREATE PLOT ###

# create vector to hold number of observations for each value of lambda, number of observations, and number of treatments
power <- matrix(NA, nrow = length(n), ncol = length(treat)*length(lambda))

# calculate number of observations for each value of lambda
for (i in 1:length(lambda)) {
  for (j in 1:length(treat)){
    for (k in 1:length(n)){
      power[k,j+2*i-2] <- vax_power2(lambda[i], n[k], treat[j])
      }
    }
  }

# create data frame that combines data to graph
df <- data.frame(cbind(power,n))

# create the graph
ggplot(df, aes(x = n)) +
  geom_line(aes(y = df[,1], color = "lambda = 0.2, treat = 5")) +
  geom_line(aes(y = df[,2], color = "lambda = 0.2, treat = 7")) +
  geom_line(aes(y = df[,3], color = "lambda = 0.25, treat = 5")) +
  geom_line(aes(y = df[,4], color = "lambda = 0.25, treat = 7")) +
  geom_line(aes(y = df[,5], color = "lambda = 0.3, treat = 5")) +
  geom_line(aes(y = df[,6], color = "lambda = 0.3, treat = 7")) +
  labs(x = 'Number of Total Completes',
       y = 'Power') +
  scale_color_manual(name = "Number of treatment arms", values = c("lambda = 0.2, treat = 5" = "darkblue", "lambda = 0.2, treat = 7" = "red", "lambda = 0.25, treat = 5" = "darkgreen", "lambda = 0.25, treat = 7" = "turquoise4", "lambda = 0.3, treat = 5" = "darkorchid4", "lambda = 0.3, treat = 7" = "springgreen2"))

Let’s say instead that we only compare each treatment arm to one other treatment arm. This approach is what we have outlined in the experimental design document. Then, the below figure contains the same information as the above figure, but we are conducting fewer hypothesis tests.

## DEFINE PARAMETERS ###

# define the possible total number of completes
n <- seq(from = 10000, to = 75000, length.out = 100 )

# outcome 
p <- 0.2 #from pilot 5 descriptive analytics analysis 595/4142

# standard deviation
sd <- sqrt(p * (1 - p))

# define the possible treatment effect sizes in *percent*
lambda <- cbind(0.2, 0.25, 0.3)

# alpha
alpha <- 0.05

# beta
beta <- 0.2

# number of treatments
treat <- cbind(5,7)

### DEFINE FUNCTION ###

vax_power2 <- function(lambda, n, t) {
  pnorm(
    ( lambda*p / (2 * sd / sqrt(n / t))) - qnorm(1-alpha/2/(3*(t-1)))
  )
}

### CREATE PLOT ###

# create vector to hold number of observations for each value of lambda, number of observations, and number of treatments
power <- matrix(NA, nrow = length(n), ncol = length(treat)*length(lambda))

# calculate number of observations for each value of lambda
for (i in 1:length(lambda)) {
  for (j in 1:length(treat)){
    for (k in 1:length(n)){
      power[k,j+2*i-2] <- vax_power2(lambda[i], n[k], treat[j])
      }
    }
  }

# create data frame that combines data to graph
df <- data.frame(cbind(power,n))

# create the graph
ggplot(df, aes(x = n)) +
  geom_line(aes(y = df[,1], color = "lambda = 0.2, treat = 5")) +
  geom_line(aes(y = df[,2], color = "lambda = 0.2, treat = 7")) +
  geom_line(aes(y = df[,3], color = "lambda = 0.25, treat = 5")) +
  geom_line(aes(y = df[,4], color = "lambda = 0.25, treat = 7")) +
  geom_line(aes(y = df[,5], color = "lambda = 0.3, treat = 5")) +
  geom_line(aes(y = df[,6], color = "lambda = 0.3, treat = 7")) +
  labs(x = 'Number of Total Completes',
       y = 'Power') +
  scale_color_manual(name = "Number of treatment arms", values = c("lambda = 0.2, treat = 5" = "darkblue", "lambda = 0.2, treat = 7" = "red", "lambda = 0.25, treat = 5" = "darkgreen", "lambda = 0.25, treat = 7" = "turquoise4", "lambda = 0.3, treat = 5" = "darkorchid4", "lambda = 0.3, treat = 7" = "springgreen2"))

7 Treatment arms

To assess the power for testing multiple correlated hypotheses, I simulate data, test the hypotheses, correct the p-values, and repeat for r repetitions, then calculate the proportion of repetitions in which I reject each hypothesis. The variables needed in the simulated data are the treatment assignment and the two outcomes: completion and clicking on the website link. In particular, we are interested in unvaccinated completes and unvaccinated website clicks. The sample for completion are clicks, while the sample for clicking on the website link are completes.

The data generation process is as follows:

  • Assign each observation to a treatment group. (treatment variable)
  • Take a draw from a binomial distribution with probability equal to the hypothesized unvax completion probability for the observation’s treatment group. (completion variable)
  • For observations who completed the survey, take a draw from a binomial distribution with probability equal to the hypothesized unvax website click probability for the observation’s treatment group. (webclick variable)

Power is calculated using the following parameters:

  • \(\alpha=0.05\)
  • replications \(r=1000\)
  • number of observations (clicks) \(n=\{75k,100k,150k,200k\}\) for a cost of \(\{11k,15k,22k,29k\}\) at $0.146 per click.

Results are reported for uncorrected p-values and for p-values corrected using Bonferonni, Holm, and Hochberg.

All pairwise comparison of treatments

First, we consider multiple hypothesis correction for evaluating all pairwise comparisons of treatments.

set.seed(94305)

### DEFINE PARAMETERS ###


outcomes1 <- matrix(c(0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.12, # Pr(unvax complete|click)
                 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.18, # Pr(unvax web click|complete)
                 1- 0.13, 1- 0.15, 1- 0.20, 1-0.17, 1-0.21, 1-.2, 1-.16), # Pr(useful free text| complete) 
                 nrow = 3, byrow = T)

outcomes2 <- outcomes1 + 0.05

outcomes_list <- list(outcomes1, outcomes2)
t_list <- lapply(outcomes_list, FUN=function(x){c(1:ncol(x))})

### POWER CALCULATIONS ###

for (i in n){
  for(o in 1:length(outcomes_list)){
    treatments <- t_list[[o]]
    outcome_matrix  <- outcomes_list[[o]]
    
    # complete r replications of n trials
    results <- replicate(r,  trial(i, outcomes = outcome_matrix, t = treatments))
    
    # create matrix to store power calculations
    power_calc <- matrix(NA, 
                         nrow = nrow(outcome_matrix)*4, 
                         ncol = (length(treatments) * (length(treatments) - 1)/2)) 
    
    # for each pairwise treatment comparison
    for(j in 1:(length(treatments) * (length(treatments) - 1)/2)){
      
      # for each mht correction
      for(k in 1:4){
        
        # calculate power for completions outcome
        power_calc[k, j] <- sum(results[k, j, ] < alpha)/r
        
        # calculate power for website clicks outcome
        power_calc[4+k,j] <- sum(results[k+4,j,]<alpha)/r
        
        # calculate power for useful information  outcome
        power_calc[8+k,j] <- sum(results[k+8,j,]<alpha)/r
        
      }
    }
    
    # Generate Graph
        power_calc <- data.frame(outcome = c(rep('Unvax completes|Click',4),
                           rep('Unvax web clicks|Complete', 4),
                           rep('Unvax useful freetext|Complete', 4)), 
                         test = rep(c("None","Bonferonni", "Holm", "Hochberg"), 3),
                         power_calc)
   
    colnames(power_calc) <- c("Outcome","Correction", 
                              'T1 v. T2', 'T1 v. T3', 'T1 v. T4','T1 v. T5','T1 v. T6','T1 v. T7','T2 v. T3','T2 v. T4','T2 v. T5','T2 v. T6','T2 v. T7','T3 v. T4','T3 v. T5','T3 v. T6','T3 v. T7','T4 v. T5','T4 v. T6','T4 v. T7','T5 v. T6','T5 v. T7','T6 v. T7')
 
    
    kable(power_calc, format = "html",  caption = paste("Power Calculations: ",i/1000," Thousand Clicks", "Matrix", o), digits = 3) %>%
      collapse_rows(columns = 1) %>%
      print()
  }  
}

Only comparisons of treatment \(t\) to treatment \(t-1\)

Then, I evaluate how adding structure to the hypothesis testing increases power. Specifically, each treatment \(t\) is designed to be compared to the treatment below it and the treatment above it. So, for five treatments, we make 4 comparisons: T1 v. T2, T2 v. T3, T3 v. T4, and T4 v. T5.

set.seed(94305)

### DEFINE PARAMETERS ###
# parameters defined in previous code chunk

for (i in n){
  for(o in 1:length(outcomes_list)){
    treatments <- t_list[[o]]
    outcome_matrix  <- outcomes_list[[o]]
    
    # complete r replications of n trials
    results <- replicate(r,  trial_subset(i, outcomes = outcome_matrix, t = treatments))
    
    # create matrix to store power calculations
    power_calc <- matrix(NA, 
                         nrow = nrow(outcome_matrix)*4, 
                         ncol = length(treatments)-1) 
    
    # for each pairwise treatment comparison
    for(j in 1:(length(treatments) - 1)){
      
      # for each mht correction
      for(k in 1:4){
        
        # calculate power for completions outcome
        power_calc[k, j] <- sum(results[k, j, ] < alpha)/r
        
        # calculate power for website clicks outcome
        power_calc[4+k,j] <- sum(results[k+4,j,]<alpha)/r
        
        # calculate power for useful information outcome
        power_calc[8+k,j] <- sum(results[k+8,j,]<alpha)/r
        
      }
    }
    
    # Generate Graph
        power_calc <- data.frame(outcome = c(rep('Unvax completes|Click',4),
                           rep('Unvax web clicks|Complete', 4),
                           rep('Unvax useful freetext|Complete', 4)), 
                         test = rep(c("None","Bonferonni", "Holm", "Hochberg"), 3),
                         power_calc)
   
    colnames(power_calc) <- c("Outcome","Correction", 
                              'T1 v. T2', 'T2 v. T3','T3 v. T4','T4 v. T5','T5 v. T6','T6 v. T7')
 
    kable(power_calc, format = "html", caption = paste("Power Calculations: ",i/1000," Thousand Clicks", "Matrix", o), digits = 3) %>%
      collapse_rows(columns = 1) %>%
      print()
  }  
}