CONQ: CONtinuous Quantile Treatment Effects for Large-Scale Online Controlled Experiments

Although in many tech firms, the vast majority of metrics are measuring user engagement activities, another important category of metrics aim to gauge app or web performance, i.e. performance metrics. For such performance metrics, not only do we care about the average users’ experience, we are also interested in the cohort of users that are suffering from the worst app or web performances. Bearing these concerns in mind, it’s important to evaluate these metrics at at least multiple quantile locations in A/B tests, especially in the tail area. At Snap Inc., P50 and P90 are the typical objectives various teams optimize for in A/B experiments. (Liu et al., 2019) proposed a statistically valid and scalable method that can calculate QTE with $p$-values and confidence intervals by appealing to quantile estimator’s asymptotic normality (Cramér, 1999), however the varying bandwidths choice mentioned is rather ad-hoc and lacks theoretically guarantee. (Deng et al., 2018) introduced an approach that appeals to delta-method and the outer CI, circumventing the density estimation which still involves user-level correlation calculation for the variance term. Furthermore, these approaches still focus on treatment effects at a number of pre-determined quantile locations (.e.g., P50 or P90), where generalizations to a continuous range of quantile locations cannot be easily obtained without point-wise recalculation.

QTE at various quantile locations has been receiving increased attention these days in the industry, as it provides a fuller picture of distributional changes between control and treatment, and better captures heterogeneity (Lux, [n.d.]). At Uber, the approach they adopted is quantile regression, citing the advantage of relying on existing abundant literature (Koenker and Hallock, 2001; Hunter and Lange, 2000). However, the methodology proposed still requires an optimization algorithm, which is not suitable for the myriad of performance metrics and experiments at Snap Inc. at scale.

In the Statistics literature, the relationship between quantiles and CDF (cumulative distribution function) has been investigated extensively. Most notably, for $X_{1},X_{2},\cdots,X_{n}$ that are i.i.d. random variables from an unknown non-parametric $F$ with density $f$, for $0<p<1$, denote by $\xi_{p}=\inf\{x:F(x)\geq p\}$ the $p$-th quantile of $F$, $\xi_{n,p}$ the $p$-th sample quantile of $F_{n}$, (Bahadur, 1966) coined the famed Bahadur representation:

To put it simply, it stated that an asymptotic relationship between quantile and CDF can be established through its density function $f$. Refinements of Bahadur’s result in the i.i.d. setting were further proposed by Kiefer (Kiefer, 1967, 1970a, 1970b). (Wu et al., 2005) further extended such representation for a wide range of dependent cases. (Woodruff, 1952) proposed the Woodruff confidence interval, which essentially established a one to one relationship between point-wise confidence intervals of the CDF and quantiles.

Another line of research focused on a $\mathcal{L}_{\infty}$-norm deviation of the empirical CDF from the truth (uniform bound), instead of at a specific quantile location (point-wise). The celebrated Dvoretsky-Kiefer-Wolfowitz-Massart (DKWM) inequality (Dvoretzky et al., 1956; Massart, 1990) provides a tight bound on the CDF for i.i.d. samples:

, which by its expression, is indifferent to quantile locations. (Howard and Ramdas, 2019) further proposed two methods for point-wise quantile confidence intervals and two methods for $\mathcal{L}_{\infty}$-norm type confidence bounds, which are proven to be adaptive in a sequential setting, resonating with the always-valid $p$-values line of research (Johari et al., 2015, 2017). However, coming from a QTE perspective, we are still most interested in point-wise statistical inference, and (Howard and Ramdas, 2019)’s method requires either parameter tuning or numeric root finding, and generalization to the dependent case is non-trivial.

In this paper, we focus on the end-horizon time setting which is the common-case at Snap Inc., instead of a sequential test setting as considered by (Howard and Ramdas, 2019). We propose an easy to interpret and implement data-driven method for a range of quantile level treatment effects for percentage changes.

At Snap Inc., hundreds of metrics are evaluated in A/B experiments simultaneously, aiming to not only gauge user telemetry data, but also on app-performance. As a social network company with a camera focus, myriad facets of in-app performances are the direct goal of optimization for many teams. There’s almost always a dual performance metric for an engagement metric. For example, as we measure the number of app opens, we also measure app open latency.

This brings about another layer of complication in A/B tests. As a typical practice, randomization units are users; for engagement type metrics, each user contributes exactly one point, and all data-points can be considered i.i.d. observations, hence statistical inference can be by and large accomplished by appealing to the Central Limit Theorem. However, the dual performance metrics are not as straight-forward, as each user contributes multiple data-points (exactly the engagement count number of data-points), therefore we have a clustered set of observations which are no longer independent. Using the same notation in 2.2, for i.i.d. samples $X_{1},\cdots,X_{n}$, the sample $p$-th quantile $\xi_{n,p}$ is asymptotically normal:

With performance metrics, typically we need delta method to approximate the numerator in the variance term (Deng et al., 2018; Liu et al., 2019). However, this is only for a fixed quantile location $p$, and it’s hard to generalize to other locations without recalculation. Furthermore, kernel density estimation is typically required for $\hat{f}(\xi_{p})$, in which its bandwidth is hard to choose. Although some ad-hoc rules could work well in practice (Liu et al., 2019), small deviations in the bandwidth affect the result greatly.

To demonstrate such issue, here is an example on a typical performance metric measuring the latency on starting the Discover section inside the Snapchat app. The study improved this latency at P50 from $911.991$ to $908.173$, with a $p$-value of $0.044$ using the method mentioned in (Liu et al., 2019) with the recommended bandwidth choice. However, if we were to vary the bandwidth using the below choices (with normal reference rule (Silverman, 1986) chosen as $h_{n}=1.06\min\{s,\text{IQR}/1.34\}/n^{1/5}$ (IQR iinter-quantile range, $s$ is sample standard deviation and $n$ is sample size), we get drastically different $p$-values demonstrated in Table 1.

Using $0.05$ as a threshold, we can see that $p$-values vary from being very significant to not significant at all as the bandwidth varies. However, we lack a robust and data-driven approach on choosing such bandwidth which can guarantee good estimation quality for inference across all quantiles.

As noted in sub-Section 2.2, the famed Bahadur representation establishes the asymptotic relationship between sample quantiles and the empirical CDF, especially in various dependent cases (Wu et al., 2005).

Such Bahadur representation has a nice geometric interpretation as illustrated by 2. At the neighborhood of point $(\xi_{p},p)$, to convert from the confidence interval of empirical CDF to quantiles, we can simply divide by the tangent of the curve $F(\xi_{p})$, which is $F^{\prime}(\xi_{p})=f(\xi_{p})$. Therefore, for i.i.d. samples, we have

which no longer involves the density estimation. For clustered sample’s case, simply replacing the variance term $p(1-p)/n$ with $\nu$ would suffice, which $\nu$ can be empirically estimated using the Bag of Little Bootstrap (BLB) technique in the next section. Take one step further, we have the result that:

Woodruff (Woodruff, 1952) first proposed this interval empirically, which ingeniously circumvents the density estimation altogether, while still provides valid statistical inference results.

Below we discuss how we extend such CI to a continuous setting and clustered sample’s case, suitable for the large-scale A/B testing need at Snap Inc.

In the previous section, we discuss the relationship between empirical CDF and quantiles through the Bahadur representation. In this section, we discuss how the variance term $\nu$ can be efficiently calculated for at a continuous range of quantile locations using one set of bootstrap results, and how Woodruff type confidence intervals can be applied.

As mentioned in (Deng et al., 2018), ’ntiles’ operation (query quantile from set of observations) is expensive and should be avoided if possible. Furthermore, bearing in mind the one to one relationship between CI for quantiles and CI for CDF, in order to simultaneously get quantile level inference for all possible locations, we can first tackle the variance of CDF, and then appeal to Woodruff CI for transforming back to quantiles. In this paper, we use the Bag of Little Bootstrap in (Kleiner et al., 2014) with the balanced method (Gleason, 1988) for resampling to further reduce Monte Carlo variance, achieving higher bootstrap efficiency.

Suppose there are $N$ users, for user $i=1,\cdots,N$, there are total of $M_{i}$ events, denoted as $X_{i,1},\cdots,X_{i,M_{i}}$, while $\sum_{i=1}^{N}M_{i}=n$. We first take the logarithm of these metrics, then preserve a pre-defined digits of precision (we chose 2): $\text{round}(\log{X_{i,M_{i}}},2),\;i=1,\cdots,N$. Note this step is to reduce the total number of unique log-scaled values for efficiency. Furthermore, due to quantile’s invariance to monotone transformation (Berger, 2013), we have

This ensures we can get back the QTE on the original metric $X$ by simply taking the exponent. Furthermore, by delta method, we have

where the first term in the last expression can be bootstrapped on log-transformed data separately for control and treatment. This way of applying the delta-method is empirically more accurate than on the original scale metric, especially when percentage change is large, or when $\xi^{C}_{f(x)}(p)$ is small.

Then we split $N$ users together with their corresponding events into $s$ subsets ($s=100$ for Snap) for control and treatment separately, all events for any given user shall be preserved in the same subset. For each subset, generate a data sketch with unique rounded log-scaled metric values (denoted by $\tilde{X}_{j},\;j=1,\cdots,K$, $K$ here is the total number of unique values in the original dataset) with its corresponding count $c_{ij},\;i=1,\cdots,s,\;j=1,\cdots,K$:

Let $B$ be the number of bootstraps, to achieve balanced bootstrapping, i.e. where all samples appear exactly $B$ times, randomly permute the long vector $\boldsymbol{V}=\left\{\underbrace{1,\cdots,1}_{\text{B repetition}},\underbrace{2,\cdots,2}_{\text{B repetition}},\cdots\cdots,\underbrace{s,\cdots,s}_{\text{B repetition}}\right\}$ and then split $\boldsymbol{V}$ into $B$ vectors of length $s$, and treat each vector $\boldsymbol{v}_{i},\;i=1,\cdots,B$ as one bootstrap sample on the buckets. For each bootstrap sample vector $\boldsymbol{v}_{i}$, we can get a counter of how many times $1,\cdots,s$ each appears, and weight corresponding $c_{i,j}$ by how many times $i$-th bucket appears in $\boldsymbol{v}_{i}$. For each bootstrap sample, we can get the empirical cdf at all unique log-scaled values $\tilde{X}_{j},\;j=1,\cdots,K$. Using all bootstrap samples, we can get an estimate of the standard deviation of empirical cdf at all $\tilde{X}_{j}$ (denoted as $\mathbb{F}_{n,i}(\tilde{X}_{j}),\;i=1,\cdots,B$) as well, which would approximate $v$ at $\mathbb{F}_{n}(\tilde{X}_{j})$:

Bootstrapping on the empirical cdf avoids the ’ntiles’ operations, and enables us to use one set of bootstrap to get all potentially interested QTE by interpolating the estimated variance via a general continuity assumption. Below is a typical example of the interpolated bootstrapped standard error across all quantiles:

Since the original Woodruff method does not always produce a symmetric CI, empirically, we propose the following conservative estimate for each unique log-scaled metric value $\tilde{X}_{j}$ for control and treatment separately:

Furthermore, using these results joined with the original observations’ empirical CDF, we can then linearly interpolate these results on a fixed grid of quantiles (ensuring common quantile locations between control and treatment for QTE evaluation on such grid). At Snap, we chose a range of P20 to P99. The same example referenced in Figure 3 has the interpolated standard error, shown in Figure 4.

Once we have such common quantile locations and estimated log-scaled metrics’ SE for control and treatment, we have all the ingredients for traditional point-wise $t$-test based statistical inference, including percentage change $p$-values, confidence intervals, etc.

The following section provides a detailed summary of the aforementioned steps constituting the algorithm named CONQ (Continuous QTE for large-scaled experiments).

Here we summarize the end-to-end CONQ algorithm in Algorithm 1, to adopt earlier notations, superscripts $C$ and $T$ denote data associated with the control and the treatment group respectively, $n$ being the total number of events, $N$ being the total number of users. As can be seen, it only requires several hyper-parameters including $d$ for log-scaling precision, $s$ for number of buckets, $B$ for balanced bootstrap iterations, a pre-determined grid on quantile locations for evaluation, and $1-\alpha$ for confidence level. As for the bag of little bootstrap part, it can be simply construed as re-weighting observations in buckets based on samples on bucket indices. This greatly reduces computational cost in comparison to the method in (Deng et al., 2018) as the number of unique observations reduce from $n$ to $K$, and $K$ here can be further adjusted by the log-scaling precision $d$.

One key step that enables us to get a range of continuous QTE evaluation is the evaluation of $\mathbb{F}^{*,id}_{n^{id},i}$ on all $\tilde{X}^{id}_{j}$ for each bootstrap sample $i$ and $id\in\{C,T\}$. This is straightforward after we do a simple sort on $\tilde{X}^{id}_{j}$ and a weighted cumulative sum on the re-weighted count $c^{id}_{i,j}$. This avoids ’ntiles’ operation and improves efficiency.

We implemented the CONQ algorithm at Snap’s A/B experimentation pipeline, the flowchart illustrated in Figure 5 demonstrates the overall workflow orchestrated by Spark (Zaharia et al., 2016) and Airflow.

Recall our earlier example in Figure 1 where we see significant but directionally different results at P50 and P90 QTE, here using CONQ, in Figure 6, we can see that the regression mainly happens at P20 to P80, i.e. lower quantiles, while improvements mainly happen after P90, thus indicating presence of heterogeneity in treatment effects, especially for devices with varying performance. Further breakdown results by varying device performance clusters (2 to 7, with lower number indicating devices with general poorer performance, while higher number indicating higher performance) in Figure 7, the results corroborate with our speculation that low regressions mainly occur for high end devices, whereas improvements occur for low end devices.

For real experiments comparison, we use method in (Deng et al., 2018) (denoted by DELTA) with kernel density bandwidth choice $h$ chosen by the ad-hoc rule proposed in (Liu et al., 2019), and compare $\Delta_{\%}$ $p$-value results at P50 and P90 with CONQ. A total of 15 performance metrics on 696 treatment and control pairs are evaluated (total of 10440 $p$-values for P50 and P90). The performance metrics included measure a wide array of latencies and crashes for the Snapchat app, which itself should be representative enough for industry A/B experiments.

We compare the resulting $p$-values from DELTA and CONQ first using a scatter-plot. As can be seen in Figure 8, these two methods align quite well and mostly adhere a linear pattern as expected, indicating stable performance of CONQ across various settings and $p$-value regions. Then we compare the proportions of significance discoveries with a static threshold on $p$-values ranging from 0.0001 to 0.2, DELTA and CONQ also perform quite similarly at both P50 and P90 levels, with CONQ resulting in slightly more discoveries across all thresholds.

In this section, we evaluate CONQ on a subset of known A/A tests at Snap where no QTE is expected across all quantile locations. Then for any given quantile location, the set of $p$-values associated QTE across different treatment control pairs and various performance metrics should have no significance after multiple testing adjustment using the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995). The reason for using BH here instead of testing for a uniform distributed $p$-value, is the positive dependency relationships between many performance metrics, and the robustness of BH under such setting. Here we choose the nominal FDR level as $\alpha=0.05,0.1,0.2$. We evaluate a total of 15 performance metrics across 54 treatment and control pair combinations. The percentile grid for evaluation is chosen as P20 to P95 with $5\%$ increment. BH is applied conditioning on a fixed percentile. Since we know that all scenarios considered are A/A tests, any discoveries by BH are false discoveries. Table 3 summarizes the percentile and their corresponding count of false discoveries and its percentage (total of 810 $p$-value per quantile). As can be seen, CONQ performs stable across various quantile locations under a realistic FDR nominal thresholds of either $0.05$ or $0.1$. With a larger threshold at $0.2$, CONQ still performs robust with slightly elevated number of false positives in the tail area.