# Covariance Estimation and its Application in Large-Scale Online Controlled
Experiments

Tao Xiong, Yihan Bao, Penglei Zhao, and Yong Wang

###### Abstract

During the last few decades, online controlled experiments (also known as A/B
tests) have been adopted as a golden standard for measuring business
improvements in industry. In our company, there are more than a billion users
participating in thousands of experiments simultaneously, and with statistical
inference and estimations conducted to thousands of online metrics in those
experiments routinely, computational costs would become a large concern. In
this paper we propose a novel algorithm for estimating the covariance of
online metrics, which introduces more flexibility to the trade-off between
computational costs and precision in covariance estimation. This covariance
estimation method reduces computational cost of metric calculation in large-
scale setting, which facilitates further application in both online controlled
experiments and adaptive experiments scenarios like variance reduction,
continuous monitoring, Bayesian optimization, etc., and it can be easily
implemented in engineering practice.

###### Index Terms:

Online Controlled Experiments, Covariance Estimation, Large-Scale Data
Processing, Variance Reduction, Continuous Monitoring, Bayesian Optimization.

##  I Introduction

Over the past few decades, online controlled experiments are proven to embody
the best scientific design for establishing the casual relationship between
treatment effect and users’ observable behaviour [1]. Although the statistical
theory underlying online controlled experiments is well developed [2], there
are difficulties in applying these methods to large-scale systems. In 2010,
Google [3] first proposed an overlapping experiment infrastructure to tackle
these problems, and later the idea of overlapping infrastructure is adopted by
many giant companies, e.g., Facebook [4], LinkedIn [5], Twitter [6], and also
our company.

Given the intrinsic complexity of overlapping infrastructures, there are two
significant difficulties in applying statistical inferences, one is how to
rigorously define “metric” and “treatment effect” in this setting; the other,
regarded as a top challenge for many companies, is how to compute the metrics
at scale reliably and efficiently [7]. We first introduce the basic concepts
of overlapping experiment infrastructure and the classical Rubin casual model
[8, 9]. Under the framework of Rubin causal model, we define the metrics in
online controlled experiments mathematically in order to better adapt to
complex real world situations and facilitate further analysis. Then, if we
want to evaluate multiple metrics, assessing the correlation within them is a
major concern. However, calculating the covariance among large-sample groups
could be computationally unacceptable [10]. In order to conquer this problem,
we propose a covariance estimation method, which helps us to reduce
computational costs, so that it is applicable in large-scale settings.
Theoretical induction and numerical simulation are both carried out to
evaluate our covariance estimation method in terms of accuracy and
computational efficiency.

Moreover, our covariance estimation method can also be applied to other tasks
in online metric analytics which rely on the estimation of covariance and
variance (variance can be considered as a special form of covariance). We
proposed 3 applications: variance reduction [11, 12, 13], continuous
monitoring [14] and Bayesian optimization [15], and we evaluate our covariance
estimation with numerical examples accordingly. Besides, the benchmarks of our
proposed method deployed in our experimentation system are also illustrated.

The main contributions of this paper are summarized as follows:

  * •

We give a rigorous definition of online metrics under the framework of Rubin
causal model, which can integrate classical causal inference theories with
modern large-scale experimental design.

  * •

We propose a novel and computational efficient method for estimating
covariance with more flexibility to trade off between computational costs and
precision in online metric analytics.

  * •

We demonstrate that our covariance estimation method can also be applied to
other tasks in A/B testing platform. To our knowledge, no other companies have
adopted such an ameliorated covariance estimation algorithm in engineering
practice.

For reproducibility of all the simulations in the paper, we provide all
related code under a public Github repository at
https://github.com/xt2357/covariancesimulation

##  II Background

###  II-A Overlapping Experiment Infrastructure

We first introduce the overlapping infrastructure in order to better
illustrate this particular application, although the proposed method itself
does not rely on the infrastructure.

####  II-A1 Basic Concepts

In the context of A/B Testing platform, for every experiment we have live
traffic as incoming data, and traffic is segmented to different experiment
groups. Typically, user is the randomization unit, where users are distributed
to groups with different treatments to study their user-level metrics in
randomized experiments [16].

To further explain the overlapping infrastructure, we introduce 3 key concepts
below:

  * •

A Domain is a segmentation of incoming live traffic.

  * •

A Layer is a partition of system parameters, where each subset of parameters
is carefully determined to ensure that the parameters in different layers have
little interaction. Experiments can be carried out within a layer.

  * •

An Experiment is a segmentation of traffic where zero or more system
parameters in the current layer can be assigned alternate values.

In overlapping infrastructure, layers and domains can be nested. Domains
contain layers, layers contain experiments, but can also contain domains; this
nesting framework introduces more flexibility to partition the system
parameters.

####  II-A2 Traffic Diversion Algorithm

Domains, layers and experiments are organized into a tree-like structure, as
illustrated in Figure 1.

Figure 1: A typical tree structure of domains and layers in overlapping
infrastructure.

With this tree structure, the incoming traffic is diverted by the
randomization unit, for example, user is the randomization unit in the our A/B
Testing platform, and every request is diverted by its unique user-id. In the
example of Figure 1, every user is randomly diverted to only one domain, but
those users in the overlapping domain are in these 3 layers simultaneously. A
randomized reallocation mechanism of users is applied to every layer of the
overlapping domain to ensure that users are independently diverted, that is to
say, every user in the overlapping domain is independently re-randomized into
the experiments within each layer, and is re-randomized into groups within the
allocated experiment.

The randomization mechanism described above can be achieved using a hash
function $f$ with a seed: for example, given the user-id of a user, we
calculate $f$(user-id, seed) modulo 10000, then all users with the same mod
will be grouped together for traffic diversion. Assume we have a set of users
with mod equals to 345, then they will all be diverted into the overlapping
domain if mod 345 is configured to be included in this domain as a stream of
traffic. With a unique seed for each randomization, we can assume that all the
randomization processes within domains, layers and experiments are
independent. For example, the seed used in the segmentation of the non-
overlapping domain and the overlapping domain will not equal to the seed used
in the randomized relocation of the UI layer. In fact, the diversion algorithm
described above defines a deterministic assignment of all potential users into
domains, layers and experiments.

###  II-B Rubin Causal Model

Our covariance estimation method is built on the framework of Rubin casual
model [8, 9, 17]. A key notion underlying the Rubin causal model is potential
outcomes. For example, a particular unit $i$ would have an outcome of
$Y_{1}(i)$ if it is exposed the treatment ($T=1$), and would have an outcome
of $Y_{0}(i)$ if it is exposed the treatment ($T=0$); potential outcome is
independent of whether it actually receives treatment or not. To measure unit-
level casual effect, we should calculate the difference of corresponding
outcomes with or without treatment, $Y_{1}(i)-Y_{0}(i)$; however, it is
impossible to observe both the potential outcomes with or without treatment at
the same time, and this dilemma is the “fundamental problem of causal
inference” [9, 18].

Although we cannot directly measure the unit-level treatment effect due to
this dilemma, we can still estimate the average treatment effect (ATE) over
the entire population in randomized experiments by calculating the difference
of means between the observations from treatment groups and the observations
from control groups, which is an unbiased estimator of ATE. [9].

Noted that there is an important assumption throughout the above procedure,
called “stable unit treatment value assumption” (SUTVA) [19], which requires
that “the observation of potential outcome on one unit should not be affected
by the particular assignment of treatments to the other units”, and this is
why we can assume all units in the treatment group are only observed with
$Y_{1}^{\prime}s$.

##  III Methodology

Under the framework of Rubin causal model, we first give the rigorous
definition of online metrics, then introduce the covariance estimation problem
of online metrics.

###  III-A Definition of Online Metrics

Suppose for a randomized experiment, there are experiment groups $1,2,...G$,
users $1,2,...,U$, and metrics $1,2,...,M$, where $G$, $U$ and $M$ are
predetermined nonrandom quantity. Then for every pair of experiment group $g$,
user $u$ and metric $m$, we define $Y(g,u,m)$ to be the corresponding
potential outcome, and $I_{g}(u)$ to be the indicator of whether user $u$ is
actually assigned to the group $g$, with $I_{g}(u)=1$ meaning user $u$ is
actually in group $g$. Another indicator $Z(g,u,m)$ denotes whether this
potential outcome $Y(g,u,m)$ is successfully observed in actual, with value 1
indicating success. Then the sum of all observed outcome for metric $m$ within
the group $g$ is defined as
$S(g,m)=\sum_{u\in\\{1,2,...U\\}}{I_{g}(u)Y(g,u,m)Z(g,u,m)}$, and the average
value of metric $m$ for the group $g$ is
$A(g,m)=\sum_{u\in\\{1,2,...U\\}}{I_{g}(u)Y(g,u,m)Z(g,u,m)}/$
$\sum_{u\in\\{1,2,...U\\}}{I_{g}(u)Z(g,u,m)}$

More generally, the calculation of most complex metrics can be derived by some
function $f$ mapping from the average metrics. We denote a complex metric for
a given group $g$ as $V_{g}$, where $V_{g}=f(A(g,1),A(g,2),...,A(g,M))$.
Moreover, in the setting of online controlled experiments, for a given group
$g$, $I_{g}(u)$ of every user $u\in\\{1,2,...U\\}$ are considered to be i.i.d.
random variables. Thus, for any group $g$, the average metrics in vector form
$\vec{A}(g)=(A(g,1),A(g,2),...,A(g,M))$ would be asymptotically (multivariate)
normal distributed in the ”big data” setting [20, 21]:

$(\vec{A}(g)-\vec{\alpha}_{g})\xrightarrow{D}\mathcal{N}(\vec{0},\Sigma_{g})$

where $\vec{\alpha}_{g}$ is the true value $\vec{A}(g)$ converges to, and
$\Sigma_{g}$ is a symmetric positive semi-definite covariance matrix for
average metrics of group $g$. Then by Delta method [21], the complex metric
$V_{g}$ is also asymptotically normal distributed when $f$ is differentiable:

$f(\vec{A}(g))\approx f(\vec{\alpha}_{g})+\nabla
f(\vec{\alpha}_{g})^{T}\cdot(\vec{A}(g)-\vec{\alpha}_{g})$

which implies that

$\left(f(\vec{A}(g))-f(\vec{\alpha}_{g})\right)\,{\xrightarrow{D}}\,\mathcal{N}\left(0,\nabla
f(\vec{\alpha}_{g})^{T}\cdot\Sigma\cdot\nabla f(\vec{\alpha}_{g})\right)$

For simplicity, we will use the average metric for illustration in the rest of
the paper, but the conclusions can be easily generalized to the other complex
metrics using differentiable functions.

###  III-B The Covariance Estimation Problem

Consider we are doing an online experiment, and we are interested to discover
the relationships between some key metrics over some specific time period. For
example, in an experiment of social media data analysis, we would like to
check if the average number of exposures to a video leads to a higher average
number of comments for this video in the last week, and we can formalize these
kind of experiments as following:

In an experiment, the average of metric $m$ for group $g$ during period $t$ is
defined as

$A_{t}(g,m)=\sum_{u}{I_{g}(u)Y_{t}(g,u,m)Z_{t}(g,u,m)}/$
$\sum_{u}{I_{g}(u)Z_{t}(g,u,m)},$

and we denote the numerator as $S_{t}(g,m)$, and denominator as $N_{t}(g,m)$.

After this rigorous definition based on potential outcome, we know that the
value of actual observation of metric $m$ for user $u$ in group $g$ over
period $t$, denoted as $X_{t}(g,u,m)$, would simply be equivalent to
$Y_{t}(g,u,m)$, given that $I_{g}(u)=1$ and $Z_{t}(g,u,m)=1$.

Now the problem is, in order to measure the relationship of metric $m_{1}$ and
$m_{2}$, how to measure their covariance efficiently?

####  III-B1 Existing Methods

In large-scale experimentation platform, observations from different users are
often considered i.i.d. samples since the population is very large and the
size of the observations is small compared to the size of the population. In
the case of estimating the covariance of two metrics in the same period $t$,
suppose we have observed two metrics value $m_{1}$ and $m_{2}$ for each user
$u$, then estimating the covariance of two sample averages is very similar to
estimating the sample covariance, which is trivial.

For estimating the covariance of two periods $t$ and $t^{\prime}$ of the same
metric, if for every user we observe both values in the two periods $t$ and
$t^{\prime}$ successfully, the covariance can also be estimated using the
sample covariance formula; but in most real-world cases, not all users use our
product in both periods $t$ and $t^{\prime}$ , and it will cause the problem
of missing data.

The naive approach is to calculate the sample covariance only with the
observations of users who appear in both periods, but it only works when the
data is missing completely at random [21], which is not a reasonable
assumption in reality, because in fact inactive users would systematically
suffer from more missing data than active users.

There are many existing works on handling missing data, one approach is based
on modelling the missing data mechanism, which requires rather strong
assumptions on moments of squares and crossproducts of elements [22]. Another
approach adopts the Delta method after data augmentation, which requires
calculating the user-level sample covariance with the augmented data [21].

We notice that all these existing methods of covariance estimation are based
on user-level data aggregation and manipulation. However, in the context of
large-scale online controlled experiments, where we have thousands of groups,
metrics, and tens of millions of users involved in different experiments with
periods of several weeks, the user-level data processing are computationally
expensive. Imaging a 30-days experiment that contains ten million users using
our company’s product every day, and continuously we collect their daily
activity data under different metrics, then even a single covariance
estimation process for only one metric between two days needs a user-level
data join of $10,000,000$ observations; and this process will repeat $30\times
30/2=450$ times for estimating pair-wise covariance among 30 days, which is a
huge requirement for computational resources in real world applications.

####  III-B2 Our Bucket-Based Method

We present a novel idea for estimating the covariance with much lower
computational costs and more flexibility, and solving the data missing problem
without the need of modelling the data missing mechanism. The method we
present is based on a preprocessing procedure, which is called “bucketing”.

In the realization of “bucketing”, we introduce an extra deterministic
randomization process for assigning all users into different buckets
$1,2,...B$. Noted that this randomization process should be independent with
the traffic diversion process described in section 2 by choosing a totally
different hash seed: given any user-id, we can calculate the hash value of the
user-id, and get the remainder $b$ of its hash value modulo $B$, then users
with the same remainder $b$ are assigned into the bucket $b\in\\{1,2,...,B\\}$
respectively. Formally, we use indicator $I_{b}(u)$ to represent the bucketing
result of the user $u$, $I_{b}(u)=1$ if and only if the user $u$ is assigned
to bucket $b$, otherwise 0.

With this bucketing algorithm, we can reduce the observations into a number of
$B$ potential outcome sums, for example, given $n$ observations of the metric
$m$ for a group $g$ in period $t$, we sum up the observations within each
bucket respectively, and the potential outcome sum for bucket $b$ with respect
to group $g$ and metric $m$ within period $t$ is denoted as
$S_{t}(g,m,b)=\sum_{u}{I_{g}(u)I_{b}(u)Y_{t}(g,u,m)Z_{t}(g,u,m)}$, it also
equals to the summation of those successful observations from users who are
assigned to bucket $b$. Given assumptions from Rubin causal model and the
independence property of the hash function used in traffic diversion and
bucketing, we prove in the following Theorem 1 that this estimator is
applicable in calculating covariance between any two average metrics, and for
every average metric it can be calculated with only $B$ potential outcome
sums, which significantly reduces computational costs when an appropriate $B$
is chosen.

For any potential outcome sums of metrics $m$ and $m^{\prime}$ across all
buckets for group $g$ at different period $t$ and $t^{\prime}$ respectively,
denoted as $S_{t}(g,m)$ and $S_{t^{\prime}}(g,m^{\prime})$, we first define
the bucket-level sample covariance as

$K(S_{t}(g,m),S_{t^{\prime}}(g,m^{\prime}))$

which equals to

$\displaystyle\frac{1}{B-1}\sum^{B}_{b=1}$
$\displaystyle{(S_{t}(g,m,b)-\bar{S}_{t}(g,m,b))}\times$
$\displaystyle{(S_{t^{\prime}}(g,m^{\prime},b)-\bar{S}_{t^{\prime}}(g,m^{\prime},b))},$

where $\bar{S}_{t}(g,m,b)$ is the average of corresponding values across $B$
buckets.

Theorem 1.
For any period of $t$ and $t^{\prime}$ for metrics $m$ and $m^{\prime}$ of
group $g$, the corresponding average metric is
$A_{t}(g,m)={S_{t}(g,m)}/{N_{t}(g,m)}$, and
$A_{t^{\prime}}(g,m^{\prime})={S_{t^{\prime}}(g,m^{\prime})}/{N_{t^{\prime}}(g,m^{\prime})}$,
and for simplicity we denote them as $A_{t}$, $A_{t^{\prime}}$, $S_{t}$,
$S_{t^{\prime}}$, $N_{t}$, $N_{t^{\prime}}$, then the estimation for its
covariance $\text{cov}[A_{t}(g,m),A_{t^{\prime}}(g,m^{\prime})]$:

$\displaystyle\text{cov}[A_{t},A_{t^{\prime}}]=$ $\displaystyle
CB[\frac{K(S_{t},S_{t^{\prime}})}{N_{t}N_{t^{\prime}}}+\frac{S_{t}S_{t^{\prime}}K(N_{t},N_{t^{\prime}})}{N_{t}^{2}N_{t^{\prime}}^{2}}$
$\displaystyle-\frac{S_{t^{\prime}}K(S_{t},N_{t^{\prime}})}{N_{t}N_{t^{\prime}}^{2}}-\frac{S_{t}K(S_{t^{\prime}},N_{t})}{N_{t}^{2}N_{t^{\prime}}}]$

where $B$ is the number of buckets and $C$ is the correction term satisfying

$C=1-E[I_{g}(u)].$

$E[I_{g}(u)]$ is the expectation of $I_{g}(u)$ across any user $u$.
Noted it does not matter which user to choose here, because $I_{g}(u)$ of
different users are i.i.d. random variables. In fact, $E[I_{g}(u)]$ is the
probability of assigning users to group $g$, and the correction term can be
approximating to 1 when the size of observations is very small compared to the
size of population. The proof of Theorem 1 is in the appendix.

###  III-C Numerical Examples and Interpretation

####  III-C1 Settings

First, we simulate $N=10,000$ users as a population using the following steps:

  * •

Generate $Y_{t}(g,u,m)$ and $Y_{t^{\prime}}(g,u,m^{\prime})$ from
$\mathcal{N}(\vec{\mu},\Sigma)$ for each user, and the covariance matrix
$\Sigma$ is constructed such that $Y_{t}(g,u,m)$ and
$Y_{t^{\prime}}(g,u,m^{\prime})$ are correlated.

  * •

For a metric $m$ in observation period $t$, define the activeness of user $u$
as the percentage of its observation value among all users in group $g$.
Formally, for any $Y_{t}(g,u,m)$, we define corresponding activeness
$e_{t}(g,u,m)=\sum_{u^{\prime}}\mathbbm{1}\\{u^{\prime}|Y_{t}(g,u^{\prime},m)<Y_{t}(g,u,m)\\}/N$,
and the same for $Y_{t^{\prime}}(g,u,m^{\prime})$.

  * •

Generate $Z_{t}(g,u,m)$ from user activeness for probabilistic modeling of
missing data, specifically, for each user $u$, there is a probability of
$\text{min}(0.5,1-e_{t}(g,u,m))$ for the corresponding observation of
$Y_{t}(g,u,m)$ to be missing, thus we can generate
$Z_{t}(g,u,m)\sim\textit{Bern}(1.0-\text{min}(0.5,1-e_{t}(g,u,m))),$ and the
same for $Z_{t^{\prime}}(g,u,m^{\prime})$.

Then we simulate the following procedure 100,000 times: each time we generate
$I_{g}(u)\sim\textit{Bern}(0.1)$ for each user in the population, and get the
corresponding observations by $Y_{t}(g,u,m)$,
$Y_{t^{\prime}}(g,u,m^{\prime})$, $Z_{t}(g,u,m)$ and
$Z_{t^{\prime}}(g,u,m^{\prime})$ for users with $I_{g}(u)=1$. A bucketing
process is then carried out, assigning each user to one of the buckets
$1,2,...B$ independently, and with the same probability of $1/B$.

For comparison, we run the naive approach, the data augmentation approach and
our bucket-based approach for estimating covariance; the true covariance can
be estimated by the sample covariance calculated from $A_{t}(g,u,m)$ and
$A_{t^{\prime}}(g,u,m^{\prime})$ among 100,000 trials, and the mean and
standard deviation of the 100,000 estimations by each approach is compared to
the true covariance for evaluating the estimation.

####  III-C2 Results and Interpretation

TABLE I: Simulation Results: Ground truth, average, standard deviation and the time for calculation of covariance in millisecond. Method | Ground Truth | Avg. | SD. | Time
---|---|---|---|---
Naive | 8.825 | 21.085 | 1.708 | 135
Data Augmentation | 8.825 | 9.790 | 0.685 | 1219
Bucketing of $B=100$ | 8.825 | 8.792 | 1.684 | 155
Bucketing of $B=200$ | 8.825 | 8.807 | 1.269 | 298
Bucketing of $B=500$ | 8.825 | 8.801 | 0.931 | 733
Bucketing of $B=1000$ | 8.825 | 8.802 | 0.790 | 1458

We can observe from Table I that with our bucket-based covariance estimator,
the average is closer to the ground truth of the covariance compared to other
methods; also the fluctuation decreases as the number of buckets increases.
However, the computational time also increases as the number of buckets
increases, thus we need to trade off between computational costs and accuracy
by choosing an appropriate number of buckets.

Moreover, we find that the data augmentation approach that assumes i.i.d.
observations is biased upward, which can be explained by the fact that we had
sampled 10% of the population distributed to different groups without
replacement, but the observations can be treated as i.i.d. only when the
sample size is small compared to the population. To better illustrate this
problem, we adopt the same setting as the previous section, and simulate the
data augmentation method with different sampling ratios, the results are
summarized in Table II.

TABLE II: Simulation Results: Data augmentation method with different sampling ratios. The results indicate that the upward bias decreases when the ratio decreases. Ratio | Ground Truth | Avg. | SD.
---|---|---|---
0.2 | 3.862 | 4.798 | 0.222
0.1 | 8.639 | 9.572 | 0.689
0.05 | 18.081 | 19.065 | 1.977
0.01 | 94.121 | 94.164 | 22.596

In real A/B testing platform, we often encounter the situation of large ratio
(more than 5%) experiment groups in real world, in which case the assumption
of independence is no longer applicable. This poses a difficulty to evaluate
the metrics accurately, but our bucket-based method is not built upon the
i.i.d. assumption of observations, thus can still work well in such
situations.

##  IV Deployment and Simulation Results

In this section, we introduce three applications of our covariance estimation
method: variance reduction, continuous monitoring, and Bayesian optimization.

###  IV-A Variance Reduction

####  IV-A1 Background

One of the major purposes for carrying out online controlled experiments is to
measure whether a treatment or strategy given to users would significantly
influence their behaviours, and it can be reflected by the change of key
performance metrics. This change is defined as “average treatment effect”
(ATE) in A/B Testing platform. To illustrate, for a given metric $m$ from
control group $g$ and experiment group $g^{\prime}$ during a specific period
$t$, the average metrics $A_{t}(g,m)$ and $A_{t}(g^{\prime},m)$ are unbiased
estimators of their corresponding true averages of potential outcomes from
population, and the estimator of ATE can be written in such form:
$\Delta_{t}(m)=A_{t}(g^{\prime},m)-A_{t}(g,m)$, which is an unbiased estimator
of $\delta=\mathbf{E}(\Delta)$.

As discussed before, the challenge with measuring ATE is the ability to detect
it when it indeed exists, usually referenced as “sensitivity”. One way to
improve sensitivity is to carry out variance reduction [11, 12, 13].

####  IV-A2 CUPED

For the purpose of variance reduction, Deng et al. [11] proposed a method
called CUPED (Controlled-experiment Using Pre-Experiment Data), in which they
introduce a control variate $Y$ with a known expectation $\theta$. Define:

$\hat{\Delta}_{t}(m)=\Delta_{t}(m)-\beta(\bar{Y}-\theta),$

where $\beta$ can be any constant. It can be easily shown that
$\hat{\Delta}_{t}(m)$ is still an unbiased estimator of $\Delta$, regardless
of the value of $\beta$. Moreover, the variance of $\hat{\Delta}_{t}(m)$ is

$\displaystyle\text{var}(\hat{\Delta}_{t}(m))$
$\displaystyle=\text{var}(\Delta_{t}(m)-\beta\bar{Y})$
$\displaystyle=\text{var}(\Delta_{t}(m))+\beta^{2}\text{var}(\bar{Y})-2\beta\text{cov}(\Delta_{t}(m),\bar{Y}),$

$\text{var}(\hat{\Delta}_{t}(m))$ reaches its minimum when
$\beta=\text{cov}(\Delta_{t}(m),\bar{Y})$ $/\text{var}(\bar{Y})$. With this
optimal $\beta$ we have:

$\text{var}(\hat{\Delta}_{t}(m))=\text{var}(\Delta_{t}(m))(1-\rho^{2}),$

where $\rho$ is widely known as the correlation coefficient. This estimator
$\hat{\Delta}_{t}(m)$ with reduced variance while keeping the unbiased
property improves the sensitivity of the online metric $m$, and allows for a
more precise assessment of the metric value. For example, in a simulation Deng
et al. carried out in the same paper [11], for metric “queries-per-user”,
using “queries-per-user in the 1-week pre-experiment period” as control
control variate $Y$, variance reduction rate can reach more than 45 percent.

The major problem here is that we need to estimate the optimal value of
$\beta$ when conducting the variance reduction procedure, and it is important
to point out that estimating $\text{cov}(\Delta_{t}(m),\bar{Y})$ can be
expensive when group size is gigantic. By applying our covariance estimation
method, the cost of estimating the covariance would reduce from joining two
data sets on user-level (in most cases larger than ten million) to bucket-
level, where typically we choose the number of bucket size less than 1000.

####  IV-A3 Numerical Examples

To illustrate the effectiveness of our covariance estimation method in
variance reduction case, we simulate the control variate $Y$ with zero mean,
and being correlated with the observations. Then, we want to evaluate the
precision of the optimal $\beta$ estimated with the proposed method under
different bucket numbers and correlation coefficients.

Let sample size $n=10,000$, and for this fixed sample size, we consider the
correlation between $\Delta_{t}(m)$ and $Y$ while $\rho$ = 0.3, 0.5, 0.6, 0.8;
and bucket number $B$ = 50, 100, 200, 500, 1000 with a repetition of 1000
times.

TABLE III: Simulation Results: Relative error of the optimal $\beta$ estimated using our method compared to the optimal value of $\beta$ in theory. Relative Error | B=50 | B=100 | B=200 | B=500 | B=1000
---|---|---|---|---|---
$\rho$=0.3 | 0.1708 | 0.1197 | 0.0817 | 0.0534 | 0.0410
$\rho$=0.5 | 0.0866 | 0.0574 | 0.0413 | 0.0278 | 0.0201
$\rho$=0.6 | 0.0630 | 0.0437 | 0.0301 | 0.0189 | 0.0142
$\rho$=0.8 | 0.0252 | 0.0174 | 0.0127 | 0.0080 | 0.0059

The results are summarized in Table III, the relative error of the estimated
optimal $\beta$ decreases as the number of buckets increases. It is also
interesting to notice that the relative error decreases rapidly when the
correlation increases, which indicates that the precision of our method would
increase when a stronger correlated control variate is chosen. (The optimal
$\beta$ approaches to zero when the correlation approaches to zero, in which
case even a small random error of the estimation will lead to a large relative
error).

###  IV-B Continuous Monitoring

####  IV-B1 Background

Hypothesis testing is a powerful tool for conducting statistical inference in
online controlled experiments, and the null-hypothesis statistical testing
(e.g., t-test or z-test) is widely used in A/B tests, where they summarize the
test result with a p-value, and reject the null hypothesis $H_{0}$ when the
p-value is less than the significant level $\alpha$. It is guaranteed that the
probability under the null hypothesis of making a Type-I error (reject $H_{0}$
when $H_{0}$ is true) is less than that pre-determined $\alpha$.

In a t-test we cannot continuously monitor the result and then early stop the
test once a significant signal is detected (p-value less than $\alpha$),
otherwise the probability of making Type I error will not be bounded by
$\alpha$ anymore. Instead, before conducting a t-test for an online-controlled
experiment, we need to determine the sufficient sample size for the t-test in
advance, and after the experiment starts, only make the conclusion when the
sample size is reached. This procedure is rather inflexible when carrying out
online controlled experiments; for example, the users might need to decide a
large sample size in advance while the true effect is easy to detect with a
much smaller sample, and the opportunity cost of waiting for extra samples can
be large in such cases. A test which allows continuous monitoring can help
users detect the true effects as quickly as possible, and adjust the sample
size dynamically with more flexibility; but continuous monitoring in the
t-test will cause a severe inflation of Type-I error, making the test an
invalid inference [23].

There are many methods proposed to solve the problem of continuous monitoring,
one approach is to construct an always-valid p-value despite continuous
monitoring from the users [23], but this requires rather strong assumptions
about the collected data. Another approach is based on Bayesian testing [14],
where given two prior probabilities $P(H_{0})$ and $P(H_{1})$ for $H_{0}$ and
$H_{1}$ to be true respectively, we update the posterior odds by

$\frac{P(H_{1}|Data)}{P(H_{0}|Data)}=\frac{P(H_{1})}{P(H_{0})}\times\frac{P(Data|H_{1})}{P(Data|H_{0})}$

The last term above is the likelihood ratio of the observed data, which is
also known as “Bayes factor”. Deng et al. proved that continuous monitoring
based on specific stopping criterion about the posterior odds can control the
false discovery rate (FDR, the proportion of false discoveries among the
discoveries) at a predetermined level. To be specific, stopping the monitoring
and rejecting $H_{0}$ when the posterior odds is greater than $K$ can
guarantee a FDR upper bound of $1/(K+1)$.

We adapted the Bayes factor method into our experimentation system by applying
the bucked-based covariance estimation described in this paper:

Observed data is distributed over many periods, i.e., each day from Day $1$ to
Day $d$ $\\{D_{1},D_{2},...,D_{d}\\}$ can be considered as different period
$t$ of the experiment, then in every period we observe
$X_{D_{1}}(g,u,m),X_{D_{2}}(g,u,m),...,X_{D_{d}}(g,u,m)$ for every user $u$ in
group $g$ of metric $m$, and the average value of observations across for
every period $t$ is $A_{t}(g,m)=\sum_{u}{X_{t}(g,u,m)}/{N_{t}(g,m)}$ for
$t\in\\{D_{1},D_{2},...,D_{d}\\}$, where $N_{t}(g,m)$ is the number of
observations in period $t$ for group $g$ of metric $m$. Then the Bayes factor
based on average metrics $A_{D_{1}}(g,m),A_{D_{2}}(g,m),...,A_{D_{d}}(g,m)$ is
calculated and users can monitor the Bayes factor continuously by each period,
and stop the experiment once the stopping criterion is satisfied.

The difficulty is that there are overlapping users among periods, the
observations of metrics from the same user from different periods are
inevitably correlated, and we must calculate the likelihood ratio of the whole
path of the observations, otherwise we will get an inaccurate likelihood and
be unable to control the FDR. Fortunately, $A_{t}(g,m)$ of different periods
can be considered (multivariate) normal distributed by applying central limit
theorem, and the covariance matrix of
$A_{D_{1}}(g,m),A_{D_{2}}(g,m),...,A_{D_{d}}(g,m)$ can be estimated
efficiently by our bucket-based covariance estimator, and the likelihood ratio
can be calculated using multivariate normal density formula.

####  IV-B2 Numerical Examples

Suppose for group $g$ we have periods of $d=30$ days and $N=4000$ users, and
for every user the average metric from Day 1 to Day $d$ is an i.i.d. random
vector with a length of $d$ from a multivariate normal distribution where the
samples of different periods are correlated, the mean of this normal
distribution is $\vec{\mu}=0$ under the null hypothesis $H_{0}$, and
$\vec{\mu}=0.3$ under $H_{1}$. We simulate 10000 runs under $H_{0}$ and other
10000 runs under $H_{1}$, with 20000 runs of the Bayes factor continuous
monitoring are simulated in total, and the prior odds of $H_{1}$ over $H_{0}$
is 1.

We can add observations from all users together to get the average metric over
all periods for group $g$, $A_{D_{1}}(g,m),$
$A_{D_{2}}(g,m),...,A_{D_{d}}(g,m)$ under both $H_{0}$ and $H_{1}$
respectively; then calculate Bayes factor for each period $t$, and reject
$H_{0}$ when the observed Bayes factor is larger than $9$, guaranteeing a FDR
bound of $1/(9+1)=0.1$.

Here we implement three methods to derive the likelihood $P(Data|H_{0})$ and
$P(Data|H_{1})$ for calculating the Bayes factor over every period $t$: the
non-covariance method is to treat
$A_{D_{1}}(g,m),A_{D_{2}}(g,m),...,A_{D_{d}}(g,m)$ as independent random
variables, where the density of every $A_{t}(g,m)$ is calculated respectively
by normal density formula (using the true variance of each $A_{t}(g,m)$) and
the likelihood of the whole path of
$A_{D_{1}}(g,m),A_{D_{2}}(g,m),...,A_{D_{d}}(g,m)$ is the multiplication of
all these densities; the second method is that we calculate the likelihood by
the multivariate normal density formula using the true covariance matrix of
$A_{D_{1}}(g,m),A_{D_{2}}(g,m),...,A_{D_{d}}(g,m)$; and in the third method,
the likelihood are calculated in the same way as the second method, except
that the covariance matrix is estimated by the our bucket-based covariance
estimator.

TABLE IV: Simulation Results: FDR is the false discovery rate which should be bounded by FDR $\leq$ 0.1 in our setting, and Power is the ratio of true rejections to the number of tests where $H_{1}$ is true. Method | FDR | Power
---|---|---
Non-Covariance | 0.179 | 0.827
True Corvariance | 0.079 | 0.693
Estimated Corvariance ($B=300$) | 0.103 | 0.742
Estimated Corvariance ($B=200$) | 0.115 | 0.753

As summarized in TABLE IV, the non-covariance method fails to control FDR
although it demonstrates the highest power, while the other methods that take
covariance into consideration show better FDR results. It is worth notice that
FDR is inflated slightly when the number of buckets reduced from 300 to 200,
which may be caused by the precision of the covariance matrix estimation, in
fact, we should choose a higher number of bucket $B$ when the number of
periods in experiments increase, because the number of parameters in the
covariance matrix grows as the periods increase. In reality, most experiments
will not be continuing more than 30 days, and a bucket number of 300 is
sufficient enough in practice, thus a huge improvement of performance in
estimating the covariance matrix is achieved when we have tens of million of
observations within each period.

###  IV-C Bayesian Optimization

####  IV-C1 Background

Bayesian optimization is a powerful method for optimizing expensive black-box
functions [24]. This method commonly starts with evaluating the objective
functions at a few randomly-selected points, and then fitting a surrogate
model to the collected data. The posterior surrogate model can provide the
estimation of the function value as well as its uncertainty at each point,
then we can construct an acquisition function that can balance between
exploitation and exploration from this posterior distribution to determine the
next query points. The optimization process proceeds sequentially several
rounds, fitting the surrogate model by all the data collected at each
iteration. In online controlled experiments, Bayesian optimization can solve
the problem of searching the optimal values of continuous parameters with
online metrics as the objective; online metrics we observe can be noisy,
taking online video playtime per-user as an example: let metric $m$ be
everyday playtime, and we have observations $X_{t}(g,u,m)$ for users
$u\in\\{1,2,...,U\\}$ in group $g$ from Day 1 to Day $d$
$\\{D_{1},D_{2},...,D_{d}\\}$, then we can define an average metric by
$A_{t}(g,m)$, and the average metric we observe in each iteration is a normal
distributed random variable.

In real-world scenarios, the experimenters would usually want a complex
objective composed of multiple online metrics, for example, we may use an
objective with the form $af_{1}+bf_{2}$, where $a$ and $b$ are weights,
$f_{1}$ and $f_{2}$ are two different online metrics, and the two metrics are
usually correlated with each other. For example, $f_{1}$ is the watch time per
user and $f_{2}$ is the number of comments per user. There is a correlation
between them because users who watch for a long time are more likely to
comment. It is expensive to estimate covariance directly when the sample size
reaches tens of millions, so the experimenters would sometimes ignore the
effect of correlation between metrics in practice. With the bucket-based
covariance estimation method, we can estimate the covariance
$\text{cov}(f_{1},f_{2})$ cost-effectively, thereby improving Bayesian
optimization’s accuracy and efficiency by taking the correlations of metrics
into consideration.

####  IV-C2 Numerical Examples

We compared the optimization performance with covariance estimation and
without covariance estimation. The objective function is
$g(\mathbf{x})=2f_{1}(\mathbf{x})+f_{2}(\mathbf{x})$ where $f_{1}(\mathbf{x})$
and $f_{2}(\mathbf{x})$ are the test function Hartmann6 with six dimensions.
Hartmann6 has six local minima and one global minimum -9.96711 in
$\mathcal{X}=\left\\{x_{i}\in(0,1),i=1,...,6\right\\}$. Our goal is to find
the minimum function value at $\mathcal{X}$ which we call best objective here.
For each point $\mathbf{x}\in\mathcal{X}$, the true function value is
$2f_{1}(\mathbf{x})+f_{2}(\mathbf{x})$, and the observation are computed by
$n=10000$ samples from a bivariate normal distribution with mean
$\mathbf{\mu}=[2f_{1}(\mathbf{x}),f_{2}(\mathbf{x})]$ and covariance matrix
$\Sigma=[[\sigma_{1}^{2},\text{cov}(f_{2}(\mathbf{x}),f_{1}(\mathbf{x}))]^{T},$
$[\text{cov}(f_{1}(\mathbf{x}),f_{2}(\mathbf{x})),\sigma_{2}^{2}]^{T}]$.
Figure 2 below shows that covariance estimation can significantly improve
Bayesian optimization’s convergence speed in both positive and negative
correlation.

In the simulation, the covariance matrix for the upper plot is
[[0.1375,0.10825318],[0.10825318,0.1125]] which represents $f_{1}$ and $f_{2}$
have a positive correlation; and the covariance matrix for lower plot is
[[0.084375,-0.11095398],[-0.11095398,0.1700625]] which represents $f_{1}$ and
$f_{2}$ have a negative correlation. We can see from the results that the
procedure with covariance estimation outperformed the procedure without
covariance estimation.

Figure 2: Value of the best objective by each iteration of optimization, with
horizontal line indicates the global optimum for the problem, and the lower
line is approaching the global optimum faster.

##  V Benchmarks in Practice

We now demonstrate the performance benchmarks for our bucket-based covariance
estimator in a subset of our company’s running experiments, which contains 345
experiments from one business scenery of our company. For every metric in our
system, there are approximately 10 billion samples from these experiments
every day, and for every metric we have a corresponding daily routine which
divert these samples into a number of $B$ (typically less than 1000) buckets
for each experiment group, then the covariance estimation of any two days can
be calculated through web services interactively and simultaneously using
these buckets reduced from samples. However, for data-join method we have to
set up an offline task which can join the two days of data together by users
and then estimate the covariance for each experiment group by the data
augmentation method described previously.

Figure 3 displays the benchmarks of estimating the pairwise covariance of one
metric within 2,3,4 days respectively. For our bucket-based method we need to
preprocess samples from every day respectively, thus the complexity grows
linearly as the number of days $n$ increases; while for the data-join method
we need to launch an offline task for each pair of days, which lead to a
complexity of $O(n^{2})$, and it will be difficult to compute when we have
thousands of running experiments, with thousands of metrics being calculated
and one billion active users participating in our system.

Figure 3: Running time for pairwise covariance estimation of one metric for
345 running experiments among different days, a corresponding spark program
which has 300 executors and 4G memory for each executor is employed in each
method.

##  VI Conclusion

Starting from the classical framework of Rubin casual model and the
independence assumptions of the hash functions in overlapping online
controlled experiment infrastructure, we propose a bucket-based covariance
estimation method and theoretical-rigorously prove that this estimator can be
generalized to calculate any covariance among different metrics in different
periods. Moreover, this method can be easily deployed in numerous applications
in any online controlled experiments platform, and we have illustrated the
effectiveness and advantages of introducing our covariance estimation method
in these applications by both simulation and real online data. The bucket-
based method can be flexible to trade off between computational costs and
precision in large-scale applications, and improve the performance of numerous
applications in online controlled experiments in aspects of efficiency,
accuracy and effectiveness. Although our covariance estimation can help tackle
large data problems, the experimental design still need to be carefully
carried out, for example, to increase the sensitivity of metrics[13, 25].
These problems are left for further research.

[Proof of the Theorem 1]

Step 1: For metrics $m$ and $m^{\prime}$ of group $g$ at $t$ and $t^{\prime}$
respectively, let $O_{t}(g,u,m):=Y_{t}(g,u,m)Z_{t}(g,u,m)$ and
$O_{t^{\prime}}(g,u,m^{\prime}):=Y_{t^{\prime}}(g,u,m^{\prime})Z_{t^{\prime}}(g,u,m^{\prime})$,
in fact, $O_{t}(g,u,m)$ is a new potential outcome which is the product of two
potential outcomes $Y_{t}(g,u,m)$ and $Z_{t}(g,u,m)$,
$O_{t^{\prime}}(g,u,m^{\prime})$ as well. We first prove that for any
$S_{t}(g,m)=\sum_{u}{I_{g}(u)O_{t}(g,u,m)}$ and
$S_{t^{\prime}}(g,m^{\prime})=\sum_{u}{I_{g}(u)O_{t^{\prime}}(g,u,m^{\prime})}$,
we have

$\begin{split}&E[B(1-E[I_{g}(u)])K(S_{t}(g,m),S_{t^{\prime}}(g,m^{\prime}))]\\\
&=\text{cov}[S_{t}(g,m),S_{t^{\prime}}(g,m^{\prime})]\end{split}$ (1)

Proof of (1). By definition we have

$S_{t}(g,m,b)=\sum_{u}{I_{g}(u)I_{b}(u)O_{t}(g,u,m)}$

where $S_{t}(g,m,b)$ is the potential outcome sum of bucket $b$ and $I_{b}(u)$
is the indicator representing if the user $u$ is assigned to bucket $b$. Thus
we have

$\begin{split}&E[K(S_{t}(g,m),S_{t^{\prime}}(g,m^{\prime}))]=E[\frac{B}{B-1}[\frac{\sum_{b}{S_{t}(g,m,b)S_{t^{\prime}}(g,m^{\prime},b)}}{B}\\\
&-\frac{\sum_{b}{S_{t}(g,m,b)}}{B}\frac{\sum_{b^{\prime}}{S_{t^{\prime}}(g,m^{\prime},b^{\prime})}}{B}]]\\\
&=E[\frac{\sum_{b}{S_{t}(g,m,b)S_{t^{\prime}}(g,m^{\prime},b)}}{B-1}]\\\
&-E[\frac{1}{B(B-1)}\sum_{b}\sum_{b^{\prime}}{S_{t}(g,m,b)S_{t^{\prime}}(g,m^{\prime},b^{\prime})}]\\\
&=E[\frac{B\sum_{b}{S_{t}(g,m,b)S_{t^{\prime}}(g,m^{\prime},b)}}{B(B-1)}]\\\
&-E[\frac{1}{B(B-1)}\sum_{b}\sum_{b^{\prime}}{S_{t}(g,m,b)S_{t^{\prime}}(g,m^{\prime},b^{\prime})}]\\\
&=E[\frac{\sum_{b}{S_{t}(g,m,b)S_{t^{\prime}}(g,m^{\prime},b)}}{B}]\\\
&-E[\frac{1}{B(B-1)}\sum_{b\neq
b^{\prime}}{S_{t}(g,m,b)S_{t^{\prime}}(g,m^{\prime},b^{\prime})}]\end{split}$
(2)

and

$\displaystyle E[S_{t}(g,m,b)S_{t^{\prime}}(g,m^{\prime},b^{\prime})]$
$\displaystyle=E[\sum_{u}{I_{g}(u)I_{b}(u)O_{t}(g,u,m)}\sum_{u}{I_{g}(u)I_{b^{\prime}}(u)O_{t^{\prime}}(g,u,m^{\prime})}]$
$\displaystyle=E[\sum_{u}{I^{2}_{g}(u)I_{b}(u)I_{b^{\prime}}(u)O_{t}(g,u,m)O_{t^{\prime}}(g,u,m^{\prime})}$
$\displaystyle+\sum_{u\neq
u^{\prime}}{I_{g}(u)I_{g}(u^{\prime})I_{b}(u)I_{b^{\prime}}(u^{\prime})O_{t}(g,u,m)O_{t^{\prime}}(g,u^{\prime},m^{\prime})}]$
$\displaystyle=E[I^{2}_{g}(u)I_{b}(u)I_{b^{\prime}}(u)]\sum_{u}{O_{t}(g,u,m)O_{t^{\prime}}(g,u,m^{\prime})}$
$\displaystyle+E[I_{g}(u)I_{g}(u^{\prime})I_{b}(u)I_{b^{\prime}}(u^{\prime})]\sum_{u\neq
u^{\prime}}{O_{t}(g,u,m)O_{t^{\prime}}(g,u^{\prime},m^{\prime})}$

Notice that $I_{b}(u)I_{b^{\prime}}(u)$ always equals to 0 when $b\neq
b^{\prime}$, and from the assumption that the hash function with different
seeds generate random results independently, and a hash function with a
specific seed generate i.i.d. results for different users, we have:

$\begin{split}&E[S_{t}(g,m,b)S_{t^{\prime}}(g,m^{\prime},b^{\prime})]\\\
&=\frac{E^{2}[I_{g}(u)]}{B^{2}}\sum_{u\neq
u^{\prime}}{O_{t}(g,u,m)O_{t^{\prime}}(g,u^{\prime},m^{\prime})}\end{split}$
(3)

if $b\neq b^{\prime}$, and

$\begin{split}&E[S_{t}(g,m,b)S_{t^{\prime}}(g,m^{\prime},b^{\prime})]\\\
&=\frac{E[I_{g}(u)]}{B}\sum_{u}{O_{t}(g,u,m)O_{t^{\prime}}(g,u,m^{\prime})}\\\
&+\frac{E^{2}[I_{g}(u)]}{B^{2}}\sum_{u\neq
u^{\prime}}{O_{t}(g,u,m)O_{t^{\prime}}(g,u^{\prime},m^{\prime})}\end{split}$
(4)

if $b=b^{\prime}$.

Apply (3) and (4) to (2) we have

$\begin{split}&E[K(S_{t}(g,m),S_{t^{\prime}}(g,m^{\prime}))]\\\
&=\frac{E[I_{g}(u)]}{B}\sum_{u}{O_{t}(g,u,m)O_{t^{\prime}}(g,u,m^{\prime})}\end{split}$
(5)

By definition we have

$\displaystyle\text{cov}[S_{t}(g,m),S_{t^{\prime}}(g,m^{\prime})]$
$\displaystyle=\text{cov}[\sum_{u}{I_{g}(u)O_{t}(g,u,m)},\sum_{u}{I_{g}(u)O_{t^{\prime}}(g,u,m^{\prime})}$
$\displaystyle=\sum_{u}\sum_{u^{\prime}}\text{cov}[I_{g}(u)O_{t}(g,u,m),I_{g}(u)O_{t^{\prime}}(g,u^{\prime},m^{\prime})]$
$\displaystyle=\sum_{u}{I_{g}(u)O_{t}(g,u,m),I_{g}(u)O_{t^{\prime}}(g,u,m^{\prime})]}$
$\displaystyle=\sum_{u}{O_{t}(g,u,m)O_{t^{\prime}}(g,u,m^{\prime})cov[I_{g}(u),I_{g}(u)]}$
$\displaystyle=E[I_{g}(u)](1-E[I_{g}(u)])\sum_{u}{O_{t}(g,u,m)O_{t^{\prime}}(g,u,m^{\prime})}$

Then compare the result above to (5), the proof is complete.

Step 2: Let $N$ be the total number of potential users in the population, we
have

$(\frac{S_{t}(g,m)}{N},\frac{N_{t}(g,m)}{N},\frac{S_{t^{\prime}}(g,m^{\prime})}{N},\frac{N_{t^{\prime}}(g,m^{\prime})}{N}):=(a,b,c,d)$

which are asymptotically multivariate normal distributed, according to the
alternative form of Delta method, we have

$\text{cov}[\frac{a}{b},\frac{c}{d}]\approx\frac{\text{cov}[a,c]}{bd}+\frac{ac\cdot\text{cov}[b,d]}{b^{2}d^{2}}-\frac{c\cdot\text{cov}[a,d]}{bd^{2}}-\frac{a\cdot\text{cov}[b,c]}{b^{2}d}$

Proof Done.

## References

  * [1] R. Kohavi, R. Longbotham, D. Sommerfield, and R. M. Henne, “Controlled experiments on the web: survey and practical guide,” _Data mining and knowledge discovery_ , vol. 18, no. 1, pp. 140–181, 2009.
  * [2] J. Manzi, _Uncontrolled: The surprising payoff of trial-and-error for business, politics, and society_. Basic Books (AZ), 2012.
  * [3] D. Tang, A. Agarwal, D. O’Brien, and M. Meyer, “Overlapping experiment infrastructure: More, better, faster experimentation,” in _Proceedings of the 16th ACM SIGKDD international conference on Knowledge discovery and data mining_ , 2010, pp. 17–26.
  * [4] E. Bakshy, D. Eckles, and M. S. Bernstein, “Designing and deploying online field experiments,” in _Proceedings of the 23rd international conference on World wide web_ , 2014, pp. 283–292.
  * [5] Y. Xu, N. Chen, A. Fernandez, O. Sinno, and A. Bhasin, “From infrastructure to culture: A/b testing challenges in large scale social networks,” in _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2015, pp. 2227–2236.
  * [6] H. S. Al-Khalifa and R. M. Al-Eidan, “An experimental system for measuring the credibility of news content in twitter,” _International Journal of Web Information Systems_ , 2011.
  * [7] S. Gupta, R. Kohavi, D. Tang, Y. Xu, R. Andersen, E. Bakshy, N. Cardin, S. Chandran, N. Chen, D. Coey _et al._ , “Top challenges from the first practical online controlled experiments summit,” _ACM SIGKDD Explorations Newsletter_ , vol. 21, no. 1, pp. 20–35, 2019.
  * [8] J. D. Angrist, G. W. Imbens, and D. B. Rubin, “Identification of causal effects using instrumental variables,” _Journal of the American statistical Association_ , vol. 91, no. 434, pp. 444–455, 1996.
  * [9] D. Rubin, “Estimating causal effects of treatments in randomized and nonrandomized studies,” _Journal of Educational Psychology_ , vol. 66, no. 5, pp. 688–701, 1974.
  * [10] J. Fan, Y. Liao, and H. Liu, “An overview of the estimation of large covariance and precision matrices,” _The Econometrics Journal_ , vol. 19, no. 1, pp. C1–C32, 2016.
  * [11] A. Deng, Y. Xu, R. Kohavi, and T. Walker, “Improving the sensitivity of online controlled experiments by utilizing pre-experiment data,” in _Proceedings of the sixth ACM international conference on Web search and data mining_ , 2013, pp. 123–132.
  * [12] Y. Guo and A. Deng, “Flexible online repeated measures experiment,” _arXiv preprint arXiv:1501.00450_ , 2015.
  * [13] H. Xie and J. Aurisset, “Improving the sensitivity of online controlled experiments: Case studies at netflix,” in _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2016, pp. 645–654.
  * [14] A. Deng, J. Lu, and S. Chen, “Continuous monitoring of a/b tests without pain: Optional stopping in bayesian testing,” in _2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA)_. IEEE, 2016, pp. 243–252.
  * [15] B. Letham, B. Karrer, G. Ottoni, E. Bakshy _et al._ , “Constrained bayesian optimization with noisy experiments,” _Bayesian Analysis_ , vol. 14, no. 2, pp. 495–519, 2019.
  * [16] S. Deng, R. Longbotham, T. Walker, and Y. Xu, “Choice of the randomization unit in online controlled experiment,” in _JSM proceedings_ , 2011.
  * [17] J. S. Sekhon, “The Neyman-Rubin Model of Causal Inference and Estimation via Matching Methods,” _Oxford handbook of political methodology_ , pp. 271–, 2008.
  * [18] P. W. Holland, “Statistics and causal inference,” _Journal of the American statistical Association_ , vol. 81, no. 396, pp. 945–960, 1986.
  * [19] D. B. Rubin, “Comment: Which ifs have causal answers,” _Journal of the American statistical association_ , vol. 81, no. 396, pp. 961–962, 1986.
  * [20] A. W. Van der Vaart, _Asymptotic statistics_. Cambridge university press, 2000, vol. 3.
  * [21] A. Deng, U. Knoblich, and J. Lu, “Applying the delta method in metric analytics: A practical guide with novel ideas,” in _Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2018, pp. 233–242.
  * [22] S. Leon, A. A. Tsiatis, and M. Davidian, “Semiparametric estimation of treatment effect in a pretest-posttest study,” _Biometrics_ , vol. 59, no. 4, pp. 1046–1055, 2003.
  * [23] R. Johari, P. Koomen, L. Pekelis, and D. Walsh, “Peeking at a/b tests: Why it matters, and what to do about it,” in _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2017, pp. 1517–1525.
  * [24] D. R. Jones, M. Schonlau, and W. J. Welch, “Efficient global optimization of expensive black-box functions,” _Journal of Global optimization_ , vol. 13, no. 4, pp. 455–492, 1998.
  * [25] R. Budylin, A. Drutsa, I. Katsev, and V. Tsoy, “Consistent transformation of ratio metrics for efficient online controlled experiments,” in _Proceedings of the Eleventh ACM International Conference on Web Search and Data Mining_ , 2018, pp. 55–63.