# SSRM: Sequential Sample Ratio Mistmatch Test
## Michael Lindon (michael.lindon@optimizely.com)

## What is an SRM?

In almost all online controlled experiments (OCEs) that use random designs, the experimenter must specify a traffic allocation for each of the variations in the experiment. If you have never done this, chances are you have been using a default traffic allocation that assigns equal amounts of traffic to all the experiment variations.

A common cause for concern is when the number of visitors assigned to each variation appears somewhat different from what one might expect under the intended traffic allocation. Indeed, this is a signal of an underlying data quality issue or engineering bug in the implementation, which has the potential to invalidate any causal conclusions drawn from the experiment. Loosely speaking, a **Sample Ratio Mistmatch** (SRM) is said to be presented when the observed sample ratios are significantly different from the expected sample ratios. It is widely regarded as a good experimentation practice to check experiments for SRMs. We can formulate this problem statistically to give a precise mathematical definition of an SRM.

## Statistical Model
The simplest example of an A/B test involves two variations, a control and a treatment, and equal traffic allocation to each. To be precise, a new visitor is assigned to the treatment with probability 0.5 upon entering the experiment. Statistically speaking, this assignment is an example of a $\text{Binomial}(1,p)$ random variable. The intended traffic allocation corresponds to the *null hypothesis:* $p = 0.5$.  Detection of a sample ratio mismatch can be formulated as a hypothesis test which tests precisely this null hypothesis.

A more general example is provided when there are multiple variations in an experiment. In this case a user assignment is an example of a $\text{Multinomial}(1,p)$ random variable $x$, which is simply a higher dimensional generalization of the $\text{Binomial}(1,p)$. If there are $d$ variations in the experiment, then $x$ is a vector of length $d$ containing $d-1$ zeros and a single one to indicate the outcome. As an example $x=(0,0,1)$ indicates that the user was assigned to the third variation. Naturally, the parameter $p$ is also a $d$ dimensional vector, with the constraint that the elements are in $[0,1]$ and sum to 1, so that these values correspond to probabilities. The intended traffic allocation to each variation corresponds to a null hypothesis about the parameter $p$. If we expected equal traffic to each variations, then the null and alternative hypotheses are

$$\begin{align}
H_0: x_1, x_2, \dots \sim \text{Multinomial}(1,p)\hspace{1cm}p=p_0,\\
H_1: x_1, x_2, \dots \sim \text{Multinomial}(1,p)\hspace{1cm}p\neq p_0,\\
\end{align}$$
with $p_0=[1/3, 1/3, 1/3]$.


There are Bayesian and frequentist approaches to testing this, each with sequential and non-sequential tests. The ssrm package contains a test developed at Optimizely which is sequential and has both Bayesian and frequentist interpretations. 

## SSRM
The wonderful thing about a sequential version of a sample ratio mismatch test is that you don't have to wait for the experiment to finish before it can be executed. The SSRM test detects in real-time, from the very beginning of an experiment, if an SRM is present. This allows bugs to be rapidly detected in the early stages of an experiment and corrected before exposing new visitors to the variations. See the accompanying [white paper](https://arxiv.org/abs/2011.03567) for the underlying mathematical details. Let's start by generating some data to work with.


```python
from scipy.stats import multinomial
import numpy as np

# Set the seed of our random number generator for reproducibility. Don't worry about this
np.random.seed(0)

# Our intended allocation probabilities
p_0 = [0.1, 0.5, 0.4]

# The actual allocation probabilities
p = [0.1, 0.49, 0.41]

# Specify number of visitors
n = 100000

# Generate allocations
data = multinomial.rvs(1, p, size=n)
```

You have just simulated `n` allocations under probabilities that differ from the intended. The null hypothesis is therefore not correct for this dataset, and it would be correct to reject the null hypothesis.

The SSRM API supports both sequential Bayesian and sequential frequentist analyses. There is a tuning parameter for both, but the above example supplies a good default that works for most use-cases.

This tutorial uses`pandas` and `matplotlib` to transform and visualize data. If you do not already have these libraries installed, uncomment the code below to install them now.


```python
import sys
!{sys.executable} -m pip install matplotlib
```

## Frequentist Analysis
To obtain sequential p-values, simply import the `sequential_p_values` function. The probability under the null hypothesis that the sequential p-value falls below $\alpha$ is less than $\alpha$. In other words, if you reject the null and declare an SRM to be present when the sequential p-value is less than $\alpha$, then the probability of reporting a **false positive** is at most $\alpha$. 


```python
from ssrm_test import ssrm_test
import matplotlib.pyplot as plt

pvals = ssrm_test.sequential_p_values(data, p_0)
```


```python
fig = plt.figure()
plt.plot(pvals)
fig.suptitle("Sequential p-value for Testing SRM", fontsize=20)
plt.xlabel("Visitor Count", fontsize=18)
plt.ylabel("p-value", fontsize=16)
plt.ylim(0, 1)
plt.show()
```

## Bayesian Analysis
One can obtain posterior probabilities (of an SRM) in the following way


```python
posterior_probabilities = ssrm_test.sequential_posterior_probabilities(data, p_0)
```


```python
fig = plt.figure()
plt.plot(posterior_probabilities)
fig.suptitle("Posterior Probability of SRM", fontsize=20)
plt.xlabel("Visitor Count", fontsize=18)
plt.ylabel("Probability", fontsize=16)
plt.ylim(0, 1)
plt.show()
```

## The Advantages of a Sequential Test
Consider the following tweet in which the user identifies an SRM using a non sequential frequentist test.

</img>

We don't have access to this exact dataset, but we can create one just like it. In total there were 1637070 visitors in the experiment and the empirical sample ratio was 50.2%. With this information we can create a similar dataset as follows.


```python
# Set the seed of our random number generator for reproducibility. Don't worry about this
np.random.seed(1)

# Our intended allocation probabilities
p_0 = [0.5, 0.5]

# The actual allocation probabilities
p = [0.502, 0.498]

# Specify number of visitors
n = 1637070

# Generate allocations
data = multinomial.rvs(1, p, size=n)

x = np.sum(data, axis=0)[0]
print(
    "Empirical sample ratios are {}% and {}%".format(
        round(100 * x / n, 1), round(100 * (n - x) / n, 1)
    )
)
```

The null hypothesis is that $p=0.5$. Let's compute the p-value from a two sided likelihood ratio test.


```python
from scipy.stats import binom

print("{} visitors assigned to the control".format(x))
print("The non sequential p-value is {}".format(2 * binom.sf(x, n, 0.5)))
```

The user has been able to identify an SRM at the 0.05-alpha level test by using a standard likelihood ratio test at the end of the experiment. Note that this happened at the **end of the experiment**. Unfortunately, this dataset is now invalid for testing and this large experiment of $1637070$ visitors has gone to waste. One needs to fix the underlying issue and start a new experiment, collecting new data.

Ideally, you wouldn't have to wait until the end of an experiment to test if there is an SRM. This is exactly the main benefit a *sequential* SRM test provides: it can test for an SRM while your experiment is running, after **every single visitor**, without inflating the false positive rate. Let's run the SSRM test on the same dataset.


```python
print("Computing... (this may take some time)")
pvals = ssrm_test.sequential_p_values(data, p_0)
fig = plt.figure()
plt.plot(pvals)
fig.suptitle("Sequential p-value for Testing SRM", fontsize=20)
plt.xlabel("Visitor Count", fontsize=16)
plt.ylabel("p-value", fontsize=16)
plt.ylim(0, 1)
plt.show()
```


```python
sig_reached = len([p for p in pvals if p > 0.05])
print(
    "SSRM rejected the null (no SRM) at the alpha=0.05 level after {} visitors.".format(
        sig_reached
    )
)
print(
    "That is only {}% of the way through the experiment.".format(
        round(100 * sig_reached / n), 3
    )
)
```

As illustrated with the above example, SSRM allowed us to detect an SRM much sooner than the end of the experiment. Larger SRM's, and hence larger errors, are detected even sooner.

## Time Aggregated Data
Sometimes you may only have time aggregated data, instead of visitor-level data. That is, you might only know the total number of visitors assigned to each variation within an interval (or bucket) of time. For 3 variations, such time aggregated data can be represented as $[[100,120,99], [123,112,111], \dots$, where there is one inner array for each time bucket. The remaining steps are identical.


```python
# Set the seed of our random number generator for reproducibility. Don't worry about this
np.random.seed(0)

# Our intended allocation probabilities
p_0 = [0.1, 0.5, 0.4]

# The actual allocation probabilities
p = [0.1, 0.49, 0.41]

# Generate allocations
data = multinomial.rvs(1000, p, size=100)

print("Now the data looks like: ")
print(data[1:10])
```


```python
pvals = ssrm_test.sequential_p_values(data, p_0)
fig = plt.figure()
plt.plot(pvals)
fig.suptitle("Sequential p-value for Testing SRM", fontsize=20)
plt.xlabel("Time Bucket Index", fontsize=16)
plt.ylabel("p-value", fontsize=16)
plt.ylim(0, 1)
plt.show()
```

## Working With Pandas DataFrames
Suppose your data is inside a `pandas` DataFrame. The following cells will show you how to manipulate them into the format expected by the SSRM API. First, if you haven't done so already, install `pandas` by uncommenting out the code below, and let's create a sample DataFrame to work with.


```python
import sys
!{sys.executable} -m pip install pandas

import pandas as pd

df = pd.DataFrame(data)
df.columns = ["variation_a", "variation_b", "variation_c"]
df.head()
```

The dataframe df gives us the total number of visitors assigned to each variation within an interval of time. It's very easy to extract these values using `df.values`. This gives an array of arrays, which can be fed into the ssrm api.


```python
pvals = ssrm_test.sequential_p_values(df.values, p_0)
fig = plt.figure()
plt.plot(pvals)
fig.suptitle("Sequential p-value for Testing SRM", fontsize=20)
plt.xlabel("Time Bucket Index", fontsize=16)
plt.ylabel("p-value", fontsize=16)
plt.ylim(0, 1)
plt.show()
```

83ce31922f5fc8dc1d82feca479b878f55588996
ipynb
Jupyter Notebook
notebooks/introduction.ipynb
michaellindon/ssrm
c44f5382471a2c2ddea4bf597072304745af8578
Apache-2.0
notebooks/introduction.ipynb
michaellindon/ssrm
c44f5382471a2c2ddea4bf597072304745af8578
Apache-2.0
notebooks/introduction.ipynb
michaellindon/ssrm
c44f5382471a2c2ddea4bf597072304745af8578
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES