How Can I Simulate Repeated Significance Testing in R

Question

I work in marketing, frequently running A/B tests on websites to determine which variation is best to serve to site visitors. I would like to be able to run a simulation that highlights the perils of repeated statistical significance tests (rather than a single significance test once sample size is met).

I'm new to R, as you will likely be able to tell with my code below. The issue I am having is, (p-value <= sig_level) is not returning the results I'd expect - ex, sometimes there are 0 p-values that meet this criteria. Clearly I have done something wrong with my code - Question: What have I done wrong?

Basic outline and R code are below.

Outline:

1. Use two variations, A & B, each with the same probability (in effect, an A/A test).
2. Determine sample size with minimal detectable effect of 10%, 80% power, 95% confidence
3. Calculate # of "conversions" by day (per variation), determined by a pre-set average daily visits (looped through by # of days it takes to reach sample size)
4. Run a significance test after each "day", store p-value

Code:

simulate_rep = function(n, conv_rate = 0.03, daily_visits, sig_level = 0.05) {
  control = 0
  variation1 = 0

  variation_visits = ceiling(daily_visits/2) # Split traffic in half for 2x variations, round up

  exp_length_days = ((n*2)/daily_visits) # n * 2 (# of variations) / daily visit count
  weeks_required = ceiling(exp_length_days / 7) # Rounded # of weeks required (to minimize day of week effects)
  exp_length_days_ceiling = (weeks_required * 7) # Rounded # of days required

  p_values = rep(NA, exp_length_days_ceiling)

  for (i in 1:exp_length_days_ceiling) {
    control_latest = rbinom(1, variation_visits, conv_rate)
    control = control + control_latest

    variation1_latest = rbinom(1, variation_visits, conv_rate)
    variation1 = variation1 + variation1_latest

    test = prop.test(c(control, variation1), c((i*variation_visits), (i*variation_visits)), alternative="two.sided")$p.value 

    p_values[i] = test
  }

  print(p_values)
}

MDE = 0.1 # Minimal detectable effect, 5%
conv_rate = 0.03 # Baseline conversion rate
conv_rate_diff = conv_rate  * (1 + MDE) # Base + MDE

conf_level = 0.95 # Confidence Level for sig test
sig_level = 1 - conf_level # Significance Level

# Determine sample size, experiment length
n = power.prop.test(p1 = conv_rate, p2 = conv_rate_diff, power = 0.8, alternative = "two.sided", sig.level = sig_level)$n # Determine necessary sample size

avg_daily_visits = 7945 # Sample daily traffic volume

output = rep(simulate_rep(n, conv_rate, avg_daily_visits, sig_level), 100)

summary(output)

sum(output <= 0.05)

Answer 1

It seems that the main problem from your code is a misunderstanding regarding the rep() function. All that rep does is repeat the elements of the vector the specified number of times. For example, rep( c(1,2,3), 3) will output the following results 1 2 3 1 2 3 1 2 3.

Here are some modifications to your code that should produce the desired results.

 simulate_rep = function(n, conv_rate = 0.03, daily_visits, sig_level = 0.05) {
   control = 0
   variation1 = 0

   variation_visits = ceiling(daily_visits/2) # Split traffic in half for 2x variations, round up

   exp_length_days = ((n*2)/daily_visits) # n * 2 (# of variations) / daily visit count
   weeks_required = ceiling(exp_length_days / 7) # Rounded # of weeks required (to minimize day of week effects)
   exp_length_days_ceiling = (weeks_required * 7) # Rounded # of days required

   p_values = rep(NA, exp_length_days_ceiling)

   for (i in 1:exp_length_days_ceiling) {
     control_latest = rbinom(1, variation_visits, conv_rate)
     control = control + control_latest

     variation1_latest = rbinom(1, variation_visits, conv_rate)
     variation1 = variation1 + variation1_latest

     test = prop.test(c(control, variation1), c((i*variation_visits), (i*variation_visits)), alternative="two.sided")$p.value 

     p_values[i] = test
   }


   return(p_values)
 }

 MDE = 0.1 # Minimal detectable effect, 5%
 conv_rate = 0.03 # Baseline conversion rate
 conv_rate_diff = conv_rate  * (1 + MDE) # Base + MDE

 conf_level = 0.95 # Confidence Level for sig test
 sig_level = 1 - conf_level # Significance Level

 # Determine sample size, experiment length
 n = power.prop.test(p1 = conv_rate, p2 = conv_rate_diff, power = 0.8, alternative = "two.sided", sig.level = sig_level)$n # Determine necessary sample size

 avg_daily_visits = 7945 # Sample daily traffic volume

 nsim = 100
 output = NULL
 for( i in 1:nsim){
    output = rbind( output, simulate_rep(n, conv_rate, avg_daily_visits, sig_level) )
 }

 # determine the proportion of p-values <= .05 by day, where each column of output corresponds to a day
 apply( output, 2, function(x) sum(x<=.05))