# AB testing with underlying segments in population

I run ab tests with the variants applied to multiple segments, say, countries. In each country, my baseline target is very different from each other.

I wonder if I should ideally model these countries in my analysis somehow, or I can just consider the ab test one global thing and proceed with a classical z/t-test.

Intuitively, if the number of countries is close to the number of samples, the variance between countries will dominate the target, so it seems clear to me that I should consider that. From the other side, if N_countries << N_samples it seems to me that we random effect of the different countries will cancel out.

This rational was very handwaving and I'm not sure where to go. I've been looking into mixed linear models but I'm not sure if that is the way to go.

• (+1) What do you mean by "...it seems to me that we random effect of the different countries will cancel out." ? I don't know what you man by "cancel out". If the response from subjects differs across countries, then there is non-independence of measurements and you need to account for that - random intercepts in a mixed model is often a good way to do that. How many countries and how many samples per country do you have ? Oct 15, 2020 at 17:45
• Your argument makes sense, but then intuitively I keep thinking all data comes from some underlying groups with some hidden non-independence, right? For example we could be talking about gender, the device or many other things. My thinking is that if we have enough data and the assignment of the variants between these groups is random, the non-independence doesn't matter, but I can't really justify that in a non-handwaiving way
– jcp
Oct 15, 2020 at 17:54
• I could have something like 10k-100k samples and 10-100 countries
– jcp
Oct 15, 2020 at 17:55
• The non independence still matters if you are interested in inference, which I assume you are, and not prediction. If you are talking about random assignment to the A/B groups that doesn't matter because subjects in one country might respond differently to subjects in another country and that non-independence still needs to be accounted for. Oct 15, 2020 at 18:33
• I'm interested in inference, but not in that depth I believe. My goal is to in the end is to find out if variant B performs better than A globally. Maybe one performs better in some countries and the other in some other countries, but one of them will be better when you average everything out. In that case should I still be interested in the non-independence? In that case, why would I not be interested in even more granularity? Like in a city level instead of country for example. Do you see where I am going?
– jcp
Oct 15, 2020 at 21:59

In this situation we have repeated measures within countries.

The research question concerns whether

variant B performs better than A globally

-- from comments in the OP

This is a question of inference.

When the response differs accross countries this gives rise to non-independence. This must be accounted for otherwise inferences will be wrong.

Maybe one performs better in some countries and the other in some other countries, but one of them will be better when you average everything out.

-- from comments in the OP

...this is incorrect. The point estimate may be unbiased but the standard errors will be wrong. This means, if there is truly no difference between A and B, then you will not know the correct probability of having observed these data, or data more extreme, and this is precisely the probability that you will be interested in.

We can see this from a simple simulation. We simulate 10,000 observations in total, from 25 countries (C), with random responses to the treatment factor (trt), where we use an effect size of 0.05 for trt

library(tidyverse)
library(lmerTest)

set.seed(15)
N <- 10000

dt <- data.frame(C = rbinom(N, 50, 0.4), trt = rbinom(N, 1, 0.5))
dt$Y <- 1 X <- model.matrix(~ trt, dt) myFormula <- "Y ~ trt + (1|C)" foo <- lFormula(eval(myFormula), dt) Z <- t(as.matrix(foo$$reTrms$$Zt)) betas <- c(0, 0.05) b <- rnorm(length(table(dt$C)), 0, 3)

dt$Y <- X %*% betas + Z %*% b + rnorm(nrow(dt)) lmerTest::lmer(myFormula, dt) %>% summary() ## Fixed effects: ## Estimate Std. Error df t value Pr(>|t|) ## (Intercept) -0.04930 0.62200 26.04703 -0.079 0.937437 ## trt 0.07032 0.01980 9972.83472 3.552 0.000384 ***  So we find that if there is actually no difference between A and B, the probability of finding the effect size of 0.07 or larger is 0.00038. However, if we ignore the non-independence, we obtain: lm(Y ~ trt, dt) %>% summary() ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.16261 0.04705 3.456 0.00055 *** ## trt 0.01087 0.06640 0.164 0.86997  which is a vastly different, and incorrect, conclusion ! In order to satisfy ourselves that this is not just a result of random variation, we can run the above code 200 times, with a different seed each time, extract the point estimates and p-values: # Monte Carlo simulation of the above n.sim <- 200 # vectors to store the point estimates and p values pe.lmm <- as.numeric(n.sim) pv.lmm <- as.numeric(n.sim) pe.lm <- as.numeric(n.sim) pv.lm <- as.numeric(n.sim) for(i in 1:n.sim) { set.seed(i) dt <- data.frame(C = rbinom(N, 50, 0.4), trt = rbinom(N, 1, 0.5)) dt$$Y <- 1 X <- model.matrix(~ trt, dt) foo <- lFormula(eval(myFormula), dt) Z <- t(as.matrix(foo$$reTrms$Zt))

b <- rnorm(length(table(dt$C)), 0, 3) dt$Y <- X %*% betas + Z %*% b + rnorm(nrow(dt))

lmm0 <- lmerTest::lmer(myFormula, dt)
lm0 <-  lm(Y ~ trt, dt)

pe.lmm[i] <- summary(lmm0)$$coefficients[2,1] # point estimate for mixed model pv.lmm[i] <- summary(lmm0)$$coefficients[2,5]  # p value for mixed model
pe.lm[i] <- summary(lm0)$$coefficients[2,1] # point estimate for linear model pv.lm[i] <- summary(lm0)$$coefficients[2,4]    # p value for linear model

}

hist(pe.lmm )


mean(pe.lmm)
## [1] 0.0522676

hist(pe.lm)


mean(pe.lm)
## [1] 0.05386004


So we can see that the point estimates are unbiased in both cases. However, when it comes to the p-values:

hist(pv.lmm )


mean(pv.lmm) ; median(pv.lmm)
## [1] 0.06440553
## [1] 0.007470051

hist(pv.lmm)


mean(pv.lm) ; median(pv.lm)
## [1] 0.406798
## [1] 0.3322566


Thus we see that the mixed model successfully returns low probabilities of obtaining these or larger point estimates if there was actually no effect of trt, the bast majority of the time. On the othe hand the linear model fails to do so. In particular, with the often-used 0.05 threshold we find:

mean(pv.lmm <= 0.05) ; mean(pv.lm <= 0.05)
## [1] 0.76
## [1] 0.175


That is, the mixed model finds a "significant" result at the 0.05 level 76% of the time, whereas the linear model, only does so 17.5% of the time.