Sampling A/B test results (revenue per visitor vectors)

Question

I have two vectors of the control version A and test version B. These vectors contain revenues by visitor. So A version has 3020 visitors who didn't purchase and B respectively 2811. Revenue data comes from different source:

A <- c(rep(0, 3020), revenue_A[, 2])
B <- c(rep(0, 2811), revenue_B[, 2])

These aren't normally distributed, but have heavy right tails. length(revenue_A[, 2]) and length(revenue_B[, 2]) are around 700 and contain values between 20 and 100.

My approach was to bootstrap these vectors 1000 times, with 10% of the values, calculate the mean revenue value and then do a t.test:

aSS <- round(0.1 * length(A))
bSS <- round(0.1 * length(B))
bootA <- c()
bootB <- c()
for (i in 1:1000) {
  tempA <- sample(A, aSS, replace = TRUE) # 10% samples of the original data
  tempB <- sample(B, bSS, replace = TRUE)
  bootA <- c(bootA, mean(tempA)) # Calculate mean of the sample
  bootB <- c(bootB, mean(tempB))
}
hist(bootA)
hist(bootB)
# --> Seem to have normal distribution, let's do t.test
t.test(bootA, bootB)

Is this the right statistical approach? I had hard time finding tutorials based on this kind of statistical calculations.

Answer 1

I believe that you are sampling the wrong distribution to calculate a p-value. Your sampling approach could be used to calculate the confidence interval (if you sample the same number of observations as the data, rather than a subset).

To calculate a p-value, you need to sample from the null distribution. In this case, the null hypothesis is no difference between the treatments, so you need to sample from that distribution. To do so, simply combine the two sets of observations, sample from that, and compare the outcome metric of interest against the observed data.

Here are some sample data (note that I modified the distribution of the non-zero data a bit):

set.seed(8675309)

A <-
  c(rep(0, 3020)
    , sample(20:100, 700
             , replace = TRUE
             , prob = 10000 - 20:100
             ))

B <-
  c(rep(0, 2811)
    , sample(20:100, 700
             , replace = TRUE
             , prob = 10000 + 20:100
             ))

We can calculate the observed difference in the means:

obs_diff <-
  mean(B) - mean(A)

For my sample, the difference is 0.95.

Then, we combine the A and B vectors to create the null distribution (that the treatment had no effect) and sample with replacement from that. Note that the sample sizes must match the data sample sizes for this to be valid.

all_values <- c(A, B)

n_rep <- 10000

null_outcomes <-
  sapply(1:n_rep, function(idx){
    mean(sample(all_values, length(B), TRUE)) -
      mean(sample(all_values, length(A), TRUE))
  })

Then, we see how often the observed data are as extreme as (or more extreme than) the null-distribution results.

mean(obs_diff >= null_outcomes)

For mine, it was as (or more) extreme than 93.97% of the null results. We can convert this to a two-tailed p-value by subtracting from 1 and doubling (skip the doubling if you are justifying a one-tailed test).

2 * (1 - mean(obs_diff >= null_outcomes) )

So, my sample data have a p-value of 0.12.

An important note: this p-value is very, very close to that calculated in a t-test of the raw data (t.test(A, B)$p.value gives 0.121). Even with skewed data, the t-test is reasonably robust when the sample size is this large. Simulation is still slightly better (it will catch the occasional weird cases), but it should give pretty similar results when the sample size is this large. More importantly, simulation allows you to use metrics other than the mean, if you need to (e.g., if the median is more important, for some reason).

Answer 2

Basically, what you are doing is a t-test on the means of resampled data.

You are estimating the bootstrapping distribution of the means of $A$ and $B$, $\bar{x}_A$ and $\bar{x}_B$, and then running a t-test between these two distributions.
This is not correct, because as you increase the number of draws, the t-test will become more and more sensitive to small differences, rejecting the null.

This is the procedure that you are using.
Assuming that the variances of the distributions of the means are different (Welch test):

$$ t=\frac{\bar{X}_A-\bar{X}_B}{\sqrt{S_A^2/N+S_B^2/N}} $$

where $S_A^2$, $S_B^2$ are the estimated variances from the distribution of the $N$ repeated resampled datasets, and $\bar{X}_A$ and $\bar{X}_B$ are their means.

The variances of the distribution of the means, $S_A^2$ and $S_B^2$, can be estimated from the variances of the original data $s_A^2$ and $s_B^2$ using the Central Limit theorem:

$$ S_A^2=s_A^2/\sqrt{n_1}\\ S_B^2=s_B^2/\sqrt{n_2} $$

where $n_1$ and $n_2$ are the original samples sizes.

This can be seen by simulating two datasets $A$ and $B$ Gamma-distributed:

set.seed(42)

n1 <- 200
n2 <- 200

# Identical distributions
x <- rgamma(n1, 0.5, 0.5)
y <- rgamma(n2, 0.5, 0.5)

t.test(x,y)  # no significant difference

    Welch Two Sample t-test

data:  x and y
t = -0.74455, df = 390.06, p-value = 0.457
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.3788515  0.1707253
sample estimates:
mean of x mean of y 
0.8889888 0.9930519 

mu_x <- mu_y <- numeric(1000)
for (i in 1:1000) {
  ind <- sample(n1, n1, replace=TRUE)
  mu_x[i] <- mean(x[ind])
  ind <- sample(n2, n2, replace=TRUE)
  mu_y[i] <- mean(y[ind])
}

# Mean of the data
(mean(x))
[1] 0.8889888
(mean(y))
[1] 0.9930519

# Mean of the means
(mean(mu_x))
[1] 0.8886424
(mean(mu_y))
[1] 0.9951539

# St. dev. of distribution of the means using central theorem formula
(sig_x <- sd(x) / sqrt(n1))
[1] 0.09151001
(sig_y <- sd(y) / sqrt(n2))
[1] 0.1056427

# St. dev of distribution of the means using the resampled means
(sig_x_1 <- sd(mu_x))
[1] 0.09012127
(sig_y_1 <- sd(mu_y))
[1] 0.1052114

# Welch test between distribution of means: now it's significant!
(res <- t.test(mu_x, mu_y))

    Welch Two Sample t-test

data:  mu_x and mu_y
t = -76.886, df = 19537, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1092270 -0.1037962
sample estimates:
mean of x mean of y 
0.8886424 0.9951539 

# Same by manual implementation

mu1 <- mean(mu_x)
mu2 <- mean(mu_y)
vx <- var(mu_x)
vy <- var(mu_y)

stderrx <- sqrt(vx/10000)
stderry <- sqrt(vy/10000)
stderr <- sqrt(stderrx^2 + stderry^2)

(tval <- (mu1-mu2) / stderr)

[1] -76.88558

# Using the central limit theorem

stderrx <- sqrt(sig_x^2/10000)
stderry <- sqrt(sig_y^2/10000)
stderr <- sqrt(stderrx^2 + stderry^2)
(tval <- (mean(x)-mean(y))/ stderr)

[1] -74.45535

which shows that you can estimate the t-statistic from resampled means without resampling!

This plot shows that as you increase the number of draws, your t-statistic becomes more and more extreme:

Bootstrapping t-test

I will basically report what you can find here on Wikipedia.
Bootstrapping can be used to perform statistical tests, but in a different way:

1. calculate the statistic: in our example the t-statistic $t$ from $x_A$ and $x_B$
2. calculate two new datasets $x'=x-\bar{x}+\bar{z}$, $y'=y-\bar{y}+\bar{z}$, where $\bar{z}$ is the global mean
3. draw a random sample $x^*$ of size $n_1$ from $x'$ and $y^*$ of size $n_2$ from $y'$, with replacement
4. calculate the t-statistic $t^*$ from $x^*$ and $y^*$
5. repeat 3 and 4 $N$ times
6. p-value: $p=\frac{\sum_{i=1}^N I[t^*\geq t]}{N}$

Some code:

mux <- mean(x)
muy <- mean(y)
muz <- mean(c(x,y))
xp <- x - mux + muz
yp <- y - muy + muz

sx <- var(x)/n1
sy <- var(y)/n2

# Observed t-statistic
t0 <- (mux - muy) / sqrt(sx+sy)

# Observed difference of means
diff_means = mux - muy

diff_mean_boot <- numeric(10000)
t_boot <- numeric(10000)
for (i in 1:10000) {
  indx <- sample(n1,n1,replace=TRUE)
  indy <- sample(n2,n2,replace=TRUE)
  diff_mean_boot[i] <- mean(xp[indx]) - mean(yp[indy])
  sxboot <- var(xp[indx])/n1
  syboot <- var(yp[indy])/n2
  t_boot[i] <- diff_mean_boot[i] / sqrt(sxboot + syboot)
}

# Estimating p-value with two-tailed test
(sum(abs(t_boot) >= abs(t0))+1)/(10000+1)
[1] 0.4645535  # similar to the original result

# Similar results for the difference of means
(sum(abs(diff_mean_boot) >= abs(diff_mean))+1)/(10000+1)
[1] 0.4566543

In case you want to avoid t-statistic, you can test for difference of means by simply calculating the statistic $\zeta=\bar{x}_A-\bar{x}_B$ (as done here).

Here you can see that the bootstrapping p-values stabilize as the number of draws increases