CUPED (variance reduction with pre-experiment covariates)

Deep dive

How it works

In: per unit an arm, an outcome Y and an optional pre-period covariate X · Out: effect, standard error, theta · State: none between runs; five sums per fit (n, sum x, sum y, sum xy, sum xx) · Cost: O(n), two passes

1. Pick X from before assignment, usually the same metric over preceding weeks; freeze the definition before looking at results.
2. Fit theta = cov(Y, X) / var(X) on all units of both arms that have X.
3. Compute the common mean mean X over the same units.
4. For each unit with X: Y' = Y - theta (X - mean X). Units without X keep Y' = Y.
5. Compare arms as in a plain A/B test: difference of means of Y', Welch standard error on Y'.
6. Report the adjusted effect and CI; the raw estimate is a check, since both target the same number.

Why variance falls. var(Y - theta X) = var Y + theta^2 var X - 2 theta cov(Y, X), minimal at theta = cov / var X, where it equals var Y (1 - rho^2) (Deng et al. 2013). The interval narrows by sqrt(1 - rho^2): 29 % at rho = 0.7.

Why it stays unbiased. Assignment is random and X was fixed earlier, so E[X] is equal in both arms and the subtracted term cancels in the difference. Estimating theta from the same data adds a small-sample bias that is minor in large samples (Lin 2013).

Pooled or per-arm theta. One pooled theta is the CUPED default. Lin (2013) fits a slope per arm (treatment-by-covariate interaction, centered covariates); with a full set of interactions it cannot hurt asymptotic precision, while the additive version can (Freedman's critique, answered in Lin 2013). EMA keeps interactions out of the primary model; pre-specify one choice.

Several covariates. With a vector X, theta = Sigma_xx^-1 Sigma_xy and the factor is 1 - R^2: ordinary least squares ANCOVA. Each slope costs about one degree of freedom: negligible at 10^5 units, real at 50 (EMA asks for few, pre-specified covariates).

Missing pre-period. New users have no X. Fit on units that have it and leave the rest unadjusted (step 4), which equals setting X = mean X. Missingness is fixed before assignment, so bias stays 0; the factor becomes 1 - f rho^2 for the share f with a pre-period.

Relation to other designs. Stratified randomization balances X at assignment; post-stratification is the same regression with a categorical X; EMA calls the two complementary. A model's prediction of Y can serve as X (Poyarkov et al. 2016, boosted trees), trained on pre-experiment data only.

Code

Simulated trial (true effect 0.1, 4000 units, rho = 0.7), pooled theta, 4000 repetitions for three cases. Asserts unbiasedness, the variance ratio against 1 - (1 - miss) rho^2, a smaller standard error and CI coverage. Compiled and run in release mode.

struct Rng(u64); // xorshift64
impl Rng {
    fn next(&mut self) -> u64 { let mut x = self.0; x ^= x << 13; x ^= x >> 7; x ^= x << 17; self.0 = x; x }
    fn unif(&mut self) -> f64 { ((self.next() >> 11) as f64 + 0.5) / (1u64 << 53) as f64 }
    fn normal(&mut self) -> f64 { // Box-Muller
        (-2.0 * self.unif().ln()).sqrt() * (6.283_185_307_179_586 * self.unif()).cos()
    }
}
struct Unit { arm: usize, x: Option<f64>, y: f64 } // x = pre-period metric; None for new users

/// Randomized trial: corr(x, y) = rho in each arm, true effect tau, a share `miss` has no pre-period data.
fn simulate(rng: &mut Rng, n: usize, tau: f64, rho: f64, miss: f64) -> Vec<Unit> {
    (0..n).map(|_| {
        let (z1, z2) = (rng.normal(), rng.normal());
        let arm = (rng.next() & 1) as usize;
        let x = if rng.unif() < miss { None } else { Some(3.0 + 1.5 * z1) };
        let y = 5.0 + 2.0 * (rho * z1 + (1.0 - rho * rho).sqrt() * z2) + tau * arm as f64;
        Unit { arm, x, y }
    }).collect()
}
/// Pooled over BOTH arms and only units that have x: theta = cov(x, y) / var(x), plus the common mean of x.
fn fit_theta(units: &[Unit]) -> (f64, f64) {
    let obs: Vec<(f64, f64)> = units.iter().filter_map(|u| u.x.map(|x| (x, u.y))).collect();
    let n = obs.len() as f64;
    let (mx, my) = (obs.iter().map(|p| p.0).sum::<f64>() / n, obs.iter().map(|p| p.1).sum::<f64>() / n);
    let cov: f64 = obs.iter().map(|p| (p.0 - mx) * (p.1 - my)).sum();
    let var: f64 = obs.iter().map(|p| (p.0 - mx).powi(2)).sum();
    (cov / var, mx)
}
/// Difference of arm means and its Welch standard error; `adjust` = Some((theta, mean_x)) applies CUPED.
fn effect(units: &[Unit], adjust: Option<(f64, f64)>) -> (f64, f64) {
    let (mut n, mut s, mut ss) = ([0.0; 2], [0.0; 2], [0.0; 2]);
    for u in units {
        let v = match (adjust, u.x) { (Some((th, mx)), Some(x)) => u.y - th * (x - mx), _ => u.y };
        n[u.arm] += 1.0; s[u.arm] += v; ss[u.arm] += v * v;
    }
    let var = |a: usize| (ss[a] - s[a] * s[a] / n[a]) / (n[a] - 1.0);
    (s[1] / n[1] - s[0] / n[0], (var(0) / n[0] + var(1) / n[1]).sqrt())
}
/// Repeat the experiment `reps` times; check unbiasedness, the variance ratio and the standard error.
fn run(rng: &mut Rng, rho: f64, miss: f64) {
    let (tau, reps, n) = (0.1, 4000, 4000);
    let (mut raw, mut adj, mut se_raw, mut se_adj, mut cover) = (vec![], vec![], 0.0, 0.0, 0);
    for _ in 0..reps {
        let units = simulate(rng, n, tau, rho, miss);
        let (th, mx) = fit_theta(&units);
        let (d0, s0) = effect(&units, None);
        let (d1, s1) = effect(&units, Some((th, mx)));
        (raw.push(d0), adj.push(d1));
        se_raw += s0 / reps as f64;
        se_adj += s1 / reps as f64;
        cover += ((d1 - tau).abs() < 1.96 * s1) as u32;
    }
    let stats = |v: &[f64]| { let m = v.iter().sum::<f64>() / v.len() as f64;
        (m, v.iter().map(|a| (a - m).powi(2)).sum::<f64>() / (v.len() as f64 - 1.0)) };
    let ((m0, v0), (m1, v1)) = (stats(&raw), stats(&adj));
    let predicted = 1.0 - (1.0 - miss) * rho * rho; // only units with a pre-period shrink
    println!("rho {rho}, missing {miss}: mean effect raw {m0:.4}, cuped {m1:.4} (true {tau})");
    println!("  variance ratio {:.3}, predicted {predicted:.3}, 1 - rho^2 = {:.3}", v1 / v0, 1.0 - rho * rho);
    println!("  mean SE raw {se_raw:.4}, cuped {se_adj:.4}; coverage {:.3}", cover as f64 / reps as f64);
    let mc = |v: f64| 4.0 * (v / reps as f64).sqrt(); // 4 Monte Carlo standard errors
    assert!((m0 - tau).abs() < mc(v0) && (m1 - tau).abs() < mc(v1)); // both unbiased
    assert!((v1 / v0 - predicted).abs() < 0.03);
    assert!(se_adj < se_raw && (se_adj / se_raw).powi(2) < predicted + 0.03);
    assert!(cover as f64 / reps as f64 > 0.935);
}
fn main() {
    let mut rng = Rng(0x9E37_79B9_7F4A_7C15);
    run(&mut rng, 0.7, 0.0);
    run(&mut rng, 0.7, 0.2); // new users without a pre-period
    run(&mut rng, 0.0, 0.0); // useless covariate: no gain, no harm
}

Measured on this Mac (Apple M4 Pro), release build, about 2 s for all three cases. rho = 0.7, none missing: variance ratio 0.501 (theory 0.510), mean standard error 0.0633 to 0.0452. With 20 % missing: 0.592 (theory 0.608), 0.0493. With rho = 0: 0.999. Mean estimates 0.0993 raw and 0.0995 adjusted for a true 0.1; coverage 0.950 to 0.953.

Complexity and constants

• O(n), five accumulators (40 bytes) per metric; p covariates cost O(n p^2) plus O(p^3). Needed units: n' = n (1 - rho^2).
• Cost drivers: joining the pre-period metric to each unit, window length, and metric count (one theta each).

Efficiency

• Cost: a streaming pass of five sums per unit and metric, cheap next to the join that supplies X; the saving is units: factor 1 - rho^2 (0.51 at rho = 0.7, arithmetic).
• Levers:
  • covariate with higher correlation (same metric, longer window) — gain rho^2 — measure rho on past data first
  • several pre-period features — factor 1 - R^2 — enough units per slope, few pre-specified features
  • stratify on the same variable at assignment — removes chance imbalance (EMA guideline) — few strata
  • precompute the per-unit pre-period metric daily — no join at analysis time — many experiments share it
  • winsorized or robust X and Y — rho not driven by outliers — heavy-tailed metrics
• Naive → standard → state of the art:
  • raw difference of means · n units per arm
  • CUPED, one pre-period metric, pooled theta · n (1 - rho^2)
  • covariate vector or model prediction, per-arm slopes, robust errors · n (1 - R^2), bounded by how predictable Y is

Implementation notes

• Compute mean X once over both arms. Centering each arm on its own mean makes the adjustment vanish when nobody is missing: you silently get the raw estimate.
• The covariate window ends at assignment, not at first exposure.
• Run an A/A test on history: adjusted effect near 0, coverage near 95 %.
• Beyond 10^7 units use Welford or pairwise sums; the textbook sum xy - n mean_x mean_y loses digits.
• Ratio metrics: adjust inside the delta method (A/B testing).

Pitfalls

• A covariate touched by the treatment (spend in the experiment week, a window that moves with each unit's first exposure) biases the estimate. EMA: no covariates measured after randomization.
• Choosing covariate or window after seeing results invalidates the p-value; pre-specify. A large gap between adjusted and raw estimates means imbalance, a bug or a leaking covariate.
• Strong curvature or outliers: rho overstates the gain; add a bin or spline of X.
• Unequal arms with one additive slope can be worse than raw (Lin 2013); use per-arm slopes.

Variants and alternatives

In production

• GrowthBook documents CUPED; Spotify's confidence library fits a linear model on pre-exposure data.
• Clinical trials: baseline-adjusted ANCOVA for continuous endpoints (EMA guideline).
• Field trials: covariance analysis on an earlier measurement (from memory).

Implementations

Sources

• A. Deng, Y. Xu, R. Kohavi, T. Walker, Improving the Sensitivity of Online Controlled Experiments by Utilizing Pre-Experiment Data, WSDM 2013. paper
• W. Lin, Agnostic notes on regression adjustments to experimental data: Reexamining Freedman's critique, Annals of Applied Statistics 7(1), 2013, 295-318. paper
• EMA CHMP, Guideline on adjustment for baseline covariates in clinical trials, EMA/CHMP/295050/2013, adopted 26 February 2015. pdf
• A. Poyarkov, A. Drutsa, A. Khalyavin, G. Gusev, P. Serdyukov, Boosted Decision Tree Regression Adjustment for Variance Reduction in Online Controlled Experiments, KDD 2016. paper