A/B testing ratio of sums

Context

Consider the following scenario for a company selling goods online. A user can purchase several items (i.e. basket of items), some of which are of particular importance and are tracked specifically (let's call them star items).

We wish to test a change in an algorithm (e.g. recommendation, risk assessment, add targeting, whatever...) which may influence both the number of star items sold, and the total sales.

• This a standard A/B test setup - the randomization unit is at the user level.
• The purpose of the A/B test is to compare the impact of the algorithm change: the control group has the original algorithm and the variant has the new algorithm
• One key metric of interest is defined as the ratio of star items sales over total sales. This is a sum of across all transactions of all users in the scope of each A or B group.
• This means that the analysis unit is at the transaction level, which is different than the randomization unit
• Metrics are calculated over the entire duration of the test (e.g. 2 weeks)

Detail of the metric used

Given group A and a set of users $$U_A = \{u_1,u_2,...,u_{N_A} \}$$, each user is involved in a number $$t_{u_n}$$ of transactions. The set $$T_A$$ of all transactions of all users in group A for the duration of the test is $$T_A = \{ t_{u_{11}}, t_{u_{12}}, ... t_{u_{nm} } \}$$.

The metric of interest for group A is defined over all transactions in scope of group A. The sums are at the transaction level, not the user level.

$$\text{Metric}_A = \frac{\sum_{t_{i} \in T_A} \text{sales star items }}{\sum_{t_{i} \in T_A} \text{sales } }$$

Of course, we could modify the definition to calculate the user level mean and that would simplify everything, but that is not the metric that is being used.

Questions

What statistical test could be used for such a metric ? An additional consideration is that although we can safely assume that users are i.i.d., it is most likely wrong to assume that individual purchase baskets are i.i.d. for the same user.

Here are a few ideas I came across, have any A/B testing practitioners come across similar metrics in the past ?

• z-test of proportions
• Bootstrapping and Jacknife
• Delta method
• Change the metric (last resort)

https://en.wikipedia.org/wiki/Ratio_estimator

Edit - Some clarifications

The reason behind this question is that I have often seen the z-test of proportion being used in this situation. Popular tools used for A/B testing often default on the proportion test and business users rarely check the underlying assumptions required for the test to be valid. @dnqxt's answer below is a good example: "Just use the z test of proportion !" - but I would like to see a rigorous statistical justification as to why (or why not) this test can be used in such a case.

I personally don't think that using a z-test of proportions would work here as sales from a purchase event are not Bernoulli trials. I argue that we cannot say that each dollar sold in the denominator can be viewed as a Bernoulli trial resulting in 0 or 1 star item dollar sold in the numerator. Moreover, since the randomization unit is at the user level, purchase events of the same user are not independent (but I would say that is a secondary issue). I could be wrong here, so please feel free to prove this otherwise !

We could also change the metric to make it a Bernoulli/Binomial which converges to Normal by using counts, but that would be a last resort solution

$$\frac{\# \text{sales with star items} }{\# \text{sales}}$$

• If the metric follows a normal distribution, perhaps you could use use a linear regression or UniANOVA. If the metric does not follow a normal distribution, you could try a transformation alike the one used in the logistic regression: $x = \log(metric / (1-metric))$, although it may have problems, when the metric is $0$ or $1$. Commented Mar 20, 2019 at 1:34
• Just to note something: Sales are usually very non-Normal. You would definitely need to take that into account when presenting your data. That being said, why not use a Beta regression on user level aggregated data? The analysis could be directly stratified with available demographic information and/or in the case of existing users with additional information about pre-existing user behaviour. (e.g. we know with good certainty where a user resides as we have delivery information and/or billing information.) Commented Mar 20, 2019 at 20:26
• Indeed, there is no reason to believe that the sales are normally distributed... dividing sales by number of users (i.e. mean sales per user) would give a normally distributed statistics given the number of users in scope. Perhaps take the log and then perform a test in the log transformed space ? Commented Mar 20, 2019 at 21:20
• In that case, this is an excellent question. (+1) Commented Mar 20, 2019 at 22:34
• Why do you analyse a derived (one dimensional and aggregrated) statistic, instead of trying to analyse more directly the data and raw observations (and then make conclusions about the derived statistic based on that analysis)? (This ratio of two sums has approximately an awkward Cauchy distribution: ie. the ratio of two sums of many little purchases where both sums are approximately a multivariate normal distributed variable) What do you wish to predict? What is the goal? What is the raw data that you have? Commented Mar 22, 2019 at 14:14

• z-test of proportions

This applies to a different case when you have binary outcomes. The z-test of proportions compares the proportions of those binary outcomes.

(Below some argument is made that you will be able to do a t-test, which for large numbers is approximately the same as the z-test. With proportions you can do a z-test because the binomial distribution has one parameter determining the variance and mean, unlike a normal distribution)

• Bootstrapping

This will be possible but not really neccesary because of the Delta method which provides the error of your observed statistic more straightforward.

• Delta method

You are interested in the ratio of two, possibly correlated variables, 1. the total sales and 2. the sales in star items.

These variables are likely asymptotically normal distributed since they are the sums of the sales from many individuals (the testing procedure could be considered to be a process like picking a sample of sales from individual users from a distribution of sales from individual users). Thus you can use the Delta method.

The use of the Delta method for the estimation of ratio's is described here. The result of this application of the Delta method actually coincides with an approximation of Hinkley's result, an exact expression for the ratio of two correlated normal distributed variables (Hinkley D.V., 1969, On the Ratio of Two Correlated Normal Random Variables, Biometrica vol. 56 no. 3).

(Sidenote: As noted by Xi'an in the comments an earlier description of the exact expression was given by George Marsaglia 1965 in the JASA Vol. 60, No. 309. A simple modern description is given in 2006 in Jstatsoft Volume 16 Issue 4)

For $$Z = \frac{X}{Y}$$ with $$\begin{bmatrix}X\\Y\end{bmatrix} \sim N\left(\begin{bmatrix} \mu_x \\ \mu_y \end{bmatrix} , \begin{bmatrix} \sigma_x^2 & \rho \sigma_x \sigma_y \\ \rho \sigma_x \sigma_y & \sigma_y^2 \end{bmatrix} \right)$$ The exact result is: $$f(z) = \frac{b(z)d(z)}{a(z)^3} \frac{1}{\sqrt{2\pi} \sigma_X\sigma_Y} \left[ \Phi \left( \frac{b(z)}{\sqrt{1-\rho^2}a(z)} \right) - \Phi \left( - \frac{b(z)}{\sqrt{1-\rho^2}a(z)} \right) \right] + \frac{\sqrt{1-\rho^2}}{\pi \sigma_X \sigma_Y a(z)^2} \exp \left( -\frac{c}{2(1-\rho^2)}\right)$$ with $$\begin{array}{} a(z) &=& \left( \frac{z^2}{\sigma_X^2} - \frac{2 \rho z}{\sigma_X \sigma_Y} + \frac{1}{\sigma_Y^2} \right) ^{\frac{1}{2}} \\ b(z) &=& \frac{\mu_X z}{ \sigma_X^2} - \frac{\rho (\mu_X+ \mu_Y z)}{ \sigma_X \sigma_Y} + \frac{\mu_Y}{\sigma_Y^2} \\ c &=& \frac{\mu_X^2}{\sigma_Y^2} - \frac{2 \rho \mu_X \mu_Y + }{\sigma_X \sigma_Y} + \frac{\mu_Y^2}{\sigma_Y^2}\\ d(z) &=& \text{exp} \left( \frac {b(z)^2 - c a(z) ^2}{2(1-\rho^2)a(z)^2}\right) \end{array}$$ And an approximation based on an assymptotic behaviour is: (for $$\mu_Y/\sigma_Y \to \infty$$): $$F(z) \to \Phi\left( \frac{z - \mu_X/\mu_Y}{\sigma_X \sigma_Y a(z)/\mu_Y} \right)$$ You end up with the Delta method result when you insert the approximation $$a(z) = a(\mu_X/\mu_Y)$$ $$a(z) \sigma_X \sigma_Y /\mu_Y \approx a(\mu_X/\mu_Y) \sigma_X \sigma_Y /\mu_Y = \left( \frac{\mu_X^2\sigma_Y^2}{\mu_Y^4} - \frac{2 \mu_X \rho \sigma_X \sigma_Y}{\mu_Y^3} + \frac{\sigma_X^2}{\mu_Y^2} \right) ^{\frac{1}{2}}$$

The values for $$\mu_X,\mu_Y,\sigma_X,\sigma_Y,\rho$$ can be estimated from your observations which allow you to estimate the variance and mean of the distribution for single users and related to this the variance and mean for the sample distribution of the sum of several users.

• Change the metric

I believe that it is interresting to do at least an intitial plot of the distribution of the sales (not the ratios) from the single users. Eventually you might end up with a situation that there is a difference between the users in group A and B, but it just happens to be not significant when you regard the single variable of the ratio (this is a bit similar to MANOVA being more powerfull than single ANOVA tests).

While the knowledge of a difference between groups, without a significant difference in the metric that you are interrested in, may not help you much in making decisions, it does help you in understanding the underlying theory and possibly design better changes/experiments next time.

Illustration

Below is a simple illustration:

Let the hypothetical distribution of sales from users be distributed as fractions $$a,b,c,d$$ which indicate how many user are of a particular case (in reality this distribution will be more complex):

                           star item sales
0$$40$$

other item sales  0$$a b 10$$      c               d


Then the sample distribution for totals from a groups with 10000 users, with for one algorithm $$a=0.190,b=0.001,c=0.800,d=0.009$$ and the other algorithm $$a=0.170,b=0.001,c=0.820,d=0.009$$ will look like:

Which shows 10000 runs drawing new users and computing the sales and ratios. The histogram is for the distribution of the ratios. The lines are computations using the function from Hinkley.

• You can see that the distribution of the two total sales numbers is approximately a multivariate normal. The isolines for the ratio show that you can estimate the ratio very well as a linear sum (as in the previous mentioned/linked linearized Delta method) and that an approximation by a Gaussian distribution should work well (and then you can use a t-test which for large numbers will be just like a z-test).
• You can also see that a scatterplot like this might provide you with more information and insight in comparison to using only the histogram.

R-Code for computing the graph:

set.seed(1)
#
#
# function to sampling hypothetic n users
# which will buy star items and/or regular items
#
#                                star item sales
#                             0$$40$$
#
#  regular item sales  0$$a b # 10$$      c               d
#
#
sample_users <- function(n,a,b,c,d) {
# sampling
q <- sample(1:4, n, replace=TRUE, prob=c(a,b,c,d))
# total dolar value of items
dri = (sum(q==3)+sum(q==4))*10
dsi = (sum(q==2)+sum(q==4))*40
# output
list(dri=dri,dsi=dsi,dti=dri+dsi, q=q)
}

#
# function for drawing those blocks for the tilted histogram
#
block <- function(phi=0.045+0.001/2, r=100, col=1) {
if (col == 1) {
bgs <- rgb(0,0,1,1/4)
cols <- rgb(0,0,1,1/4)
} else {
bgs <- rgb(1,0,0,1/4)
cols <- rgb(1,0,0,1/4)
}
angle <- c(atan(phi+0.001/2),atan(phi+0.001/2),atan(phi-0.001/2),atan(phi-0.001/2))
rr <- c(90000,90000+r,90000+r,90000)
x <- cos(angle)*rr
y <- sin(angle)*rr
polygon(x,y,col=cols,bg=bgs)
}
block <- Vectorize(block)

#
# function to compute Hinkley's density formula
#
fw <- function(w,mu1,mu2,sig1,sig2,rho) {
#several parameters
aw <- sqrt(w^2/sig1^2 - 2*rho*w/(sig1*sig2) + 1/sig2^2)
bw <- w*mu1/sig1^2 - rho*(mu1+mu2*w)/(sig1*sig2)+ mu2/sig2^2
c <- mu1^2/sig1^2 - 2 * rho * mu1 * mu2 / (sig1*sig2) + mu2^2/sig2^2
dw <- exp((bw^2 - c*aw^2)/(2*(1-rho^2)*aw^2))

# output from Hinkley's density formula
out <- (bw*dw / ( sqrt(2*pi) * sig1 * sig2 * aw^3)) * (pnorm(bw/aw/sqrt(1-rho^2),0,1) - pnorm(-bw/aw/sqrt(1-rho^2),0,1)) +
sqrt(1-rho^2)/(pi*sig1*sig2*aw^2) * exp(-c/(2*(1-rho^2)))

out
}
fw <- Vectorize(fw)

#
# function to compute
# theoretic distribution for sample with parameters (a,b,c,d)
# lazy way to compute the mean and variance of the theoretic distribution
fwusers <- function(na,nb,nc,nd,n=10000) {
users <- c(rep(1,na),rep(2,nb),rep(3,nc),rep(4,nd))
dsi <- c(0,40,0,40)[users]
dri <- c(0,0,10,10)[users]
dti <- dsi+dri

sig1 <- sqrt(var(dsi))*sqrt(n)
sig2 <- sqrt(var(dti))*sqrt(n)
cor <- cor(dti,dsi)
mu1 <- mean(dsi)*n
mu2 <- mean(dti)*n

w <- seq(0,1,0.001)
f <- fw(w,mu1,mu2,sig1,sig2,cor)
list(w=w,f=f,sig1 = sig1, sig2=sig2, cor = cor, mu1= mu1, mu2 = mu2)
}

# sample many ntr time to display sample distribution of experiment outcome
ntr <- 10^4

# sample A
dsi1 <- rep(0,ntr)
dti1 <- rep(0,ntr)
for (i in 1:ntr) {
users <- sample_users(10000,0.19,0.001,0.8,0.009)
dsi1[i] <- users$$dsi dti1[i] <- users$$dti
}

# sample B
dsi2 <- rep(0,ntr)
dti2 <- rep(0,ntr)
for (i in 1:ntr) {
users <- sample_users(10000,0.19-0.02,0.001,0.8+0.02,0.009)
dsi2[i] <- users$$dsi dti2[i] <- users$$dti
}

# hiostograms for ratio
ratio1 <- dsi1/dti1
ratio2 <- dsi2/dti2
h1<-hist(ratio1, breaks = seq(0, round(max(ratio2+0.04),2), 0.001))
h2<-hist(ratio2, breaks = seq(0, round(max(ratio2+0.04),2), 0.001))

# plotting

plot(0, 0,
xlab = "sum of total sales", ylab = "sum of star item sales",
xlim = c(82000,92000),
ylim = c(2500,6000),
pch=21, col = rgb(0,0,1,1/10), bg = rgb(0,0,1,1/10))
title("sample distribution for sum of 10 000 users")

# isolines
brks <- seq(0, round(max(ratio2+0.02),2), 0.001)
for (ls in 1:length(brks)) {
col=rgb(0,0,0,0.25+0.25*(ls%%5==1))
lines(c(0,10000000),c(0,10000000)*brks[ls],lty=2,col=col)
}

# scatter points
points(dti1, dsi1,
pch=21, col = rgb(0,0,1,1/10), bg = rgb(0,0,1,1/10))
points(dti2, dsi2,
pch=21, col = rgb(1,0,0,1/10), bg = rgb(1,0,0,1/10))

# diagonal axis
phi <- atan(h1breaks) r <- 90000 lines(cos(phi)*r,sin(phi)*r,col=1) # histograms phi <- h1$$mids r <- h1$$density*10 block(phi,r,col=1) phi <- h2$$mids r <- h2$$density*10 block(phi,r,col=2) # labels for histogram axis phi <- atan(h1$$breaks)[1+10*c(1:7)] r <- 90000 text(cos(phi)*r-130,sin(phi)*r,h1$$breaks[1+10*c(1:7)],srt=-87.5,cex=0.9) text(cos(atan(0.045))*r-400,sin(atan(0.045))*r,"ratio of sum of star items and sum of total items", srt=-87.5,cex=0.9) # plotting functions for Hinkley densities using variance and means estimated from theoretic samples distribution wf1 <- fwusers(190,1,800,9,10000) wf2 <- fwusers(170,1,820,9,10000) rf1 <- 90000+10*wf1$$f phi1 <- atan(wf1$$w) lines(cos(phi1)*rf1,sin(phi1)*rf1,col=4) rf2 <- 90000+10*wf2$$f phi2 <- atan(wf2$$w) lines(cos(phi2)*rf2,sin(phi2)*rf2,col=2)  • Fantastic answer Martijn, as always. So to recap, once I use the Delta method approximation, or the more accurate Hinkley result which includes correlation, I obtain a new r.v. which is asymptotically normal. I can then perform a standard z-test of difference in means on this new r.v. for control and variant ? i.e. H0 be both control and variant have the same mean, or mean = 0 ? I have actually never seen an example of using the delta method for a test of difference of means... do you have an example ? Commented Mar 28, 2019 at 19:39 • The Delta method also includes correlation. Hinkley is only more accurate but the difference is not so important when the variance of your variables is small. You could say that$$x/y-\mu_x/\mu_y\approx(x-\mu_x)\frac{1}{\mu_y}+(y-\mu_y)\frac{\mu_x}{\mu_y^2}$$ie you approximate the ratio error by a linear sum for which the variance can be easily computed (also with correlation). When the distribution spreads out over a larger area then those iso-lines are not any more approximately equidistant and parallel, then the difference between the Delta method and Hinkley's formula becomes larger. Commented Mar 28, 2019 at 20:23 • Lets continue this discussion in chat ? chat.stackexchange.com/rooms/91683/ab-test-ratio-of-sums thanks ! Commented Mar 28, 2019 at 20:30 • @usεr11852 The distribution of the sales from single users does not need to be normal. Like in the example where it is some sort of multivariate Bernoulli distribution. What you do need is that the sample distribution of the sum of sales of users needs to be (approximately) normal. This will be the case when the distribution for single users has finite moments (then they will vanish for the sum). They will indeed vanish faster when you do not have very high skew or fat tails, and you will be needing a sufficient sample size (but the same is true for bootstrapping). Commented Mar 29, 2019 at 9:34 • As a minor remark (on a great answer), George Marsaglia had an earlier publication(1969) on the ratio density, as “Ratios of Normal Variables and Ratios of Sums of Uniform Variables.” Journal of the American Statistical Association, 60, 193–204. Hinkley (1969) wrongly pointed out a mistake in that paper, corrected in Hinley (1970 Biometrika). Commented Mar 25, 2021 at 8:14 What is described is a classical case of an A/B experiment where we have dependence between users and items (purchases here); we need to account for this because otherwise we will have a biased estimate of the associated variance. To counter that, we either bootstrap by taking account user/item clusters or we go fully parametric using a mixed model. The paper by Bakshy & Eckles (2013) Uncertainty in Online Experiments with Dependent Data: An Evaluation of Bootstrap Methods is an invaluable reference on that matter that focuses on online A/B tests. Looking into the things in more detail, in some ways the B&E paper is a extended use-case of the Owen's (2007) paper The pigeonhole bootstrap. Similarly the mixed effect approach is based on the uber-classic paper by Bayeen et al. (2009) Mixed-effects modeling with crossed random effects for subjects and items. To briefly comment in the methods you mention: as you correctly identify the $$z$$-test on proportion is oversimplifying; it will assume IID and as the B&E paper presents, that assumption can be seriously misleading. The same limitation extends to bootstrap, if we ignore the structure of the data. Regarding the Delta method: deviations from normality and/or noisy data usually make the Delta method suboptimal in comparison to bootstrap approaches (e.g. Hole 2007 A comparison of approaches to estimating confidence intervals for willingness to pay measures) but I have seen some recent papers (e.g. Deng et al. 2017 Trustworthy analysis of online A/B tests: Pitfalls, challenges and solutions and Deng et al. 2018 Applying the Delta Method in Metric Analytics: A Practical Guide with Novel Ideas) that seem promising; notice that there is the implicit assumption that the average treatment effect will be normal. Finally, the change of metric is also an excellent idea when reasonable. We should not be afraid to push for changes that are more mathematically coherent just because a metric is already in place. In conclusion: If there is an over-arching scheme on all the papers I cited above is that we need to ensure that our unit of analysis and our unit of randomisation are aligned with our research questions. • +1 fantastic answer thanks a lot - would you argue that if the IID assumptions were held, then the z-proportion test would work in this case ? How would you demonstrate that this is a sum of Bernoulli R.V's ? Commented Mar 21, 2019 at 0:53 • I am glad I could help. I would not usez\$-proportions test in an A/B test setting because I do think that IID assumptions are realistically impossible. Research is booming (as you can see) so they are multiple better alternatives! Commented Mar 21, 2019 at 0:56
• I am following up on your answer - indeed the i.i.d assumptions for purchase events of the same user are strong. But more fundamentally - how could we justify using the proportion test on such a metric - wouldn't this require each event to be a 0-1 Bernoulli trial ? That would require considering each dollar of sales as an independent event ? Commented Mar 26, 2019 at 18:37
• Yes, it would but that clearly is a very strong assumption/awkward simplification. That being said, A/B tests often examine binary outcomes (eg. churned/not churned) so the use of proportions tests might be appropriate. As a matter of convenience, sometimes one might use a test for continious outcomes on discrete data (or vice versa) and the results may not outright misleading. As you correctly note, that does not make them sound. Commented Mar 28, 2019 at 0:54

A very simple approach would be to use a permutation test. This is is distribution-free test so you don't have to worry about distributions.

The idea is simple. You randomly shuffle the labels and count how many times the measured difference of your metric of interest is larger than the difference you got from your real data. The ratio you get is your p-value.

Why does it work ? Well, if the null hypothesis is True, then randomly shuffling the labels (i.e being in group A or B) would very often yield better values then the one you measured. However, if your recommendation algorithm works, then randomly shuffling would rarely yield better result than the one you got.

You may also use bootstrapping to get confidence intervals on your metric of both your A and B group. This, two, is allowed without assumptions on your distribution. That's not equivalent of a statistical test (even if the CI intervals do not cross), but the visual aspect of "bars + error" bars can be interesting for your team.

I have answered a very similar question (where you indeed found me). How to test the difference in two proportions when the outcomes aren't binary? . I think the code I suggested there apply here too.

p1 <- sum(sales_data[target_control==1,"final_value"])/sum(sales_data[target_control==1,"initial_value"])
p2 <- sum(sales_data[target_control==0,"final_value"])/sum(sales_data[target_control==0,"initial_value"])
yourGap<-abs(p1-p2)
L<-sales_data["target_control"]==1

count=0
for ( i in 1:10000) {
Lperm=sample(L)
if (abs(p1_perm-p2_perm)>=yourGap) {
count=count+1
}
}
pvalue=count/10000

• I think that using a straightforward permutation test approach as the one you describe will require IID measurements. From my understanding of the problem, this does not hold and therefore a (vanilla) permutation test will lead to anti-conservative tests. As there is user-related (and probabably item) clustering the IID assumption can lead to variance under-estimation as the individual treatment effect variation is large. Please see some of the references (B&E 2013 or D et al. 2017) in my answer for more details. (They deal with bootstrap by the general ideas transfer to permutation tests.) Commented Mar 22, 2019 at 22:10
• Thanks for this comment. I am not totally sure about this i.i.d problem. Unless I missunderstood, with the metric he suggested, the problem can be boiled down to a single value by user, which dissolves the i.i.d problem. Commented Mar 23, 2019 at 6:04
• Unfortunately that is not the case. The OP mentions it at least twice: once "... we could modify the definition to calculate the user level mean..." (i.e. this is not a user level) and later on when saying: "purchase events of the same user are not independent". Permuting the purchase event labels will thus be incoherent as they are not IID. (Which is the point of the papers I cite. :) ) Commented Mar 24, 2019 at 0:55
• I see. This makes this answer quite irrelevant and unless I find some useful spin I think of deleting it. I distanced myself from statistics during the two last years, but these papers you references sound interesting. Definitely on my reading list in case I find interest again in the field. Commented Mar 24, 2019 at 18:37
• Thank you for your answers - what about bootstrapping or permuting at the user level - i.e. ensuring that for each permutation , all or none of the sales for a given user are included. This could help reduce the issues of independence below the user level. Would that work better ? Commented Mar 26, 2019 at 18:39

Bootstrap confidence intervals would be my choice of technique for this scenario. I would like to outline an approach with some example choice numbers that you can use and the reasoning behind the approach:

1. You have two pots/bags and each bag contains people from the control group and variant group: $$U_{ctr}$$ and $$U_{var}$$ with respective sizes $$N_{ctr}$$ and $$N_{var}$$. I changed your notation slightly, I hope that's OK.
2. You select a random sample of $$k$$ people from both groups with replacement. If both of your populations are "sufficiently large" (say at least 2000 users for example) you can choose $$k\leq N_{ctr}$$ and $$k \leq N_{var}$$. Rule of thumb: I usually select $$k = \frac{min(N_{ctr},N_{var})}{5}$$ for more flexible results but in most bagging (bootstrap aggregating) algorithms, the default option is to sample from the whole population. If your populations are smaller then you can still do that but make sure you select a "sufficiently large" $$k$$ (say at least 400 users), again by sampling with replacement. Let's say note them $$SU_{ctr}$$ and $$SU_{var}$$ both of size $$k$$
3. You calculate your metric by getting all transactions that each person in each group made by looking at the original $$T$$ transactions dataset for each user in $$SU_{ctr}$$ and $$SU_{var}$$. You will then end up with $$Metric_{ctr1}$$ and $$Metric_{ctr2}$$. Store these values. Important notice: you should calculate these values by summing the total star items sales and dividing by the total sales. Don't take the average of each person's individual sales basket. This is very important as this is the metric you are looking at.
4. Go back to point 2 and iterate. The optimal number of bootstrap samples, $$B$$, that you can choose depends on many factors but again, a good rule of thumb would be around 1000 times.
5. You now have a $$B$$ amount of $$Metric_{ctr}$$ and same for $$Metric_{var}$$. You can now choose to compare their means using many of the usual techniques. I would personally choose to construct confidence intervals and see if they overlap or an independent sample t-test AND complete the analysis with some histograms/density plots & boxplots.

Off-topic personal opinion: Always choose to viz things such as distributions whenever possible, we have the power to do that nowadays. The tests above are totally fine but there are cases where they might go wrong. For instance if you choose $$B$$ to be extremely high, say 1000000, then even the smallest difference between the means is more likely to be flagged as significant.

The above is robust because no matter what underlying distribution is, the central limit theorem ensures that if $$B$$ is sufficiently large, both means of $$Metric_{var}$$ and $$Metric_{ctr}$$ across the samples will be normally distributed and the tests will be valid. You will witness that from the visuals as well. Any concerns about underlying distributions of different user spends etc will be dealt by the CLT.

There are plenty of references and good reads by the users before me. Furthermore, there is a lot of research conducted on the optimal example numbers I mentioned above you can look into it to. I just wanted to give you a more empirical and easy to understand outline of approach that is robust. You can start with that and see whether things change by changing the above example numbers.

• Excellent ! This is exactly what I was looking for thanks a lot ! Using a different approach, any thoughts on how one could also apply the Delta method to this problem ? stats.stackexchange.com/questions/291594/… Commented Mar 27, 2019 at 1:24
• Also - would you have some resources on the appropriate choices for B and SU ? In this (fictional) example we could imagine that there are 1M users in control and say 100k users ij variant Commented Mar 27, 2019 at 1:46
• General references: Bible of Bootstrap is amazon.co.uk/… but there is an excellent summary here: galton.uchicago.edu/~eichler/stat24600/Handouts/bootstrap.pdf Commented Mar 27, 2019 at 10:14
• Regarding the SU: In bagging (bootstrap aggregating) the state of the art (and default option for many bagging algorithms) is to sample from both populations keeping the original sizes untouched. There is a theorem that states that we expect 63.2% unique and the rest duplicated users so feel free to select N_ctr and N_var sizes just sample from there. I select 1/5 as I usually want more flexible estimates than stable ones. It's a bias-variance trade-off question once again! Commented Mar 27, 2019 at 10:14
• Regarding the B: Here things are even more unclear. In wikipedia's page, they go for an example of B=1000, in the elem. of stat. learn web.stanford.edu/~hastie/Papers/ESLII.pdf they try 100, 200 in 2000 observations. In this ensemble modelling case onlinelibrary.wiley.com/doi/full/10.1111/coin.12198 they try 500 and 1000. and here below 100 statistics.berkeley.edu/sites/default/files/tech-reports/… I'd say that this number depends on your machine capabilities. In practice, numbers above 1000 don't make huge difference but any choice above 500 I'd say is fair. Commented Mar 27, 2019 at 10:17

regression with cluster robust standard errors will be a solution.

assume raw data df like this:

randomization_id,  analysis_id, treat, y
1,            1,     1, 0
1,            2,     1, 1
1,            3,     1, 0
2,            1,     0, 1
2,            2,     0, 0

library(sandwich)
library(lmtest)
mod <- lm(y ~ treat, df)
coeftest(mod, vcov = vcovCL, cluster = ~randomization_id)


this will give exactly the same value as Delta method.

The advantage of this approach is that it can easily handle more complex situations, for example, control for the imbalance of covariates to improve statistical power.

reference: Deng, A., Lu, J., & Qin, W. (2021). The equivalence of the Delta method and the cluster-robust variance estimator for the analysis of clustered randomized experiments. arXiv preprint arXiv:2105.14705.