Some business metrics need to be normalized. For example, a video streaming service might ask how many rebuffers a user experienced. But if the user played many videos, it will be better to normalize this metric by number of videos played. Another example is when we sum total revenue for a supplier, but we also need to normalize it for the number of product lines.

In previous posts, we have talked about how to compute relative effects, which is another type of normalization. In this case we divided the treatment mean by the control mean, both derived from a single model like OLS. For a ratio metric we divide 2 completely different random variables, derived from different models.

Bivariate Model in OLS

Say the numerator metric is \(N\), and the denominator metric is \(D\). When we first studied how to divide a treatment mean by a control mean, we said that these two quantities have no covariance. It makes sense when the underlying OLS model is simply \(y = \beta_0 + T \beta_1 + \varepsilon\) and \(T\) is randomized. However, when thinking about two distinct metrics, which might affect one another, we will need to jointly model them.

Consider the bivariate OLS model that jointly models

\[[N, D] = \beta_0 + T \beta_1 + X \beta_2 + \varepsilon.\]

In this case, the outcome \(Y = [N, D]\) is a \(1 \times 2\) row vector, so \(\beta_0 = [\beta_{0,N}, \beta_{0,D}]\) is not a constant, it’s a vector. Likewise \(\beta_1 = [\beta_{1,N}, \beta_{1,D}]\) is a vector, and when \(X\) has \(p\) columns \(\beta_2 = [\beta_{2,N}, \beta_{2,D}]\) is a \(p \times 2\) matrix. The noise \(\varepsilon = [\varepsilon_N, \varepsilon_D]\) is also a vector, with covariance matrix \(\Sigma_\varepsilon\).

Under this model, we define the mean of each metric in each arm by averaging the model over the covariates of all users

\[\begin{align} \mu_{N,T} &= E_X\big[E[N | T = 1, X]\big] = \beta_{0,N} + \beta_{1,N} + E[X] \beta_{2,N} \\ \mu_{D,T} &= E_X\big[E[D | T = 1, X]\big] = \beta_{0,D} + \beta_{1,D} + E[X] \beta_{2,D} \\ \mu_{N,C} &= E_X\big[E[N | T = 0, X]\big] = \beta_{0,N} + E[X] \beta_{2,N} \\ \mu_{D,C} &= E_X\big[E[D | T = 0, X]\big] = \beta_{0,D} + E[X] \beta_{2,D} \end{align}\]

The treatment effect on the ratio metric is then

\[\begin{align} \tau &= \frac{\mu_{N,T}}{\mu_{D,T}} - \frac{\mu_{N,C}}{\mu_{D,C}} \end{align}\]

The point estimate of the effect is derived from the coefficients of the models, replacing \(E[X]\) with the sample average \(\bar{X}\)

\[\begin{align} \frac{\hat{\mu}_{N,T}}{\hat{\mu}_{D,T}} &= \frac{\hat{\beta}_{0,N} + \hat{\beta}_{1,N} + \bar{X} \hat{\beta}_{2,N}}{\hat{\beta}_{0,D} + \hat{\beta}_{1,D} + \bar{X} \hat{\beta}_{2,D}} \\ \frac{\hat{\mu}_{N,C}}{\hat{\mu}_{D,C}} &= \frac{\hat{\beta}_{0,N} + \bar{X} \hat{\beta}_{2,N}}{\hat{\beta}_{0,D} + \bar{X} \hat{\beta}_{2,D}} \end{align}\]

Note that \(\tau\) is built from a ratio of means: we first average the numerator and the denominator within each arm, then divide. This is different from averaging a per-user ratio, which has the same geometric-mean subtlety that we saw with the log linear model in the post on relative effects.

Covariance of a bivariate model

We are already familiar with how to model the covariance of the parameter estimates under a single outcome. In this formulation there are two error terms, and covariance between the two. The clever step here is if the covariates for \(N\) are the same as the covariates for \(D\), e.g. \(Z = [1, T, X]\) then we can express

\[\begin{align} \text{Cov}(\text{vec}(\hat{B})) &= \Sigma_\varepsilon \otimes (Z^\top Z)^{-1} = \begin{pmatrix} \sigma_N^2 (Z^\top Z)^{-1} & \sigma_{ND} (Z^\top Z)^{-1} \\ \sigma_{ND} (Z^\top Z)^{-1} & \sigma_D^2 (Z^\top Z)^{-1} \end{pmatrix} \\ \Sigma_\varepsilon &= \begin{pmatrix} \sigma_N^2 & \sigma_{ND} \\ \sigma_{ND} & \sigma_D^2 \end{pmatrix} \\ \hat{\sigma}_N^2 &= \frac{1}{n-k} \sum_{i=1}^n \hat{\varepsilon}_{N,i}^2 \\ \hat{\sigma}_D^2 &= \frac{1}{n-k} \sum_{i=1}^n \hat{\varepsilon}_{D,i}^2 \\ \hat{\sigma}_{ND} &= \frac{1}{n-k} \sum_{i=1}^n \hat{\varepsilon}_{N,i} \hat{\varepsilon}_{D,i} \end{align}\]

Variance of the Treatment Effect with the Delta Method

We will get the variance of \(\hat{\tau}\) in two steps. First we apply the delta method within each arm to get the variance of each arm’s ratio. Then we apply it again to the difference of the two ratios.

Step 1: The Ratio Within Each Arm

The estimated means in each arm are linear combinations of the coefficients. Define the baseline vectors \(b_T = (1, 1, \bar{X})\) and \(b_C = (1, 0, \bar{X})\). Then

\[\begin{align} \hat{\mu}_{N,T} &= b_T \hat{\beta}_N & \hat{\mu}_{D,T} &= b_T \hat{\beta}_D \\ \hat{\mu}_{N,C} &= b_C \hat{\beta}_N & \hat{\mu}_{D,C} &= b_C \hat{\beta}_D \end{align}\]

To get the covariance of these means, it helps to recall the case of a single outcome. For the regression \(y = Z\beta + \varepsilon\) with \(\text{Var}(\varepsilon) = \sigma^2\), the estimated mean of arm \(t \in \{T, C\}\) is simply \(\hat{\mu}_t = b_t \hat{\beta}\). Its variance can be compactly written in terms of a fixed variable from the observed data, and \(\sigma^2\).

\[\begin{align} c_t &= b_t (Z^\top Z)^{-1} b_t^\top \\ \text{Var}(\hat{\mu}_t) &= b_t \text{Cov}(\hat{\beta}) b_t^\top \\ &= b_t \sigma^2 (Z^\top Z)^{-1} b_t^\top \\ &= \sigma^2 c_t \end{align}\]

In the bivariate model, we need the full \(2 \times 2\) covariance of \((\hat{\mu}_{N,t}, \hat{\mu}_{D,t}) = (b_t \hat{\beta}_N, b_t \hat{\beta}_D)\). The elements are governed by

\[\begin{align} \text{Var}(\hat{\mu}_{N,t}) &= b_t \text{Cov}(\hat{\beta}_N) b_t^\top &= b_t \sigma_N^2 (Z^\top Z)^{-1} b_t^\top &= c_t \sigma_N^2 \\ \text{Var}(\hat{\mu}_{D,t}) &= b_t \text{Cov}(\hat{\beta}_D) b_t^\top &= b_t \sigma_D^2 (Z^\top Z)^{-1} b_t^\top &= c_t \sigma_D^2 \\ \text{Cov}(\hat{\mu}_{N,t}, \hat{\mu}_{D,t}) &= b_t \text{Cov}(\hat{\beta}_N, \hat{\beta}_D) b_t^\top &= b_t \sigma_{ND} (Z^\top Z)^{-1} b_t^\top &= c_t \sigma_{ND} \end{align}\]

Hence we can write \[\text{Cov}\begin{pmatrix} \hat{\mu}_{N,t} \\ \hat{\mu}_{D,t} \end{pmatrix} = c_t \, \Sigma_\varepsilon.\]

The ratio in arm \(t\) is \(R_t = g(\mu_{N,t}, \mu_{D,t}) = \frac{\mu_{N,t}}{\mu_{D,t}}\), with gradient

\[\nabla g_t = \left(\frac{1}{\mu_{D,t}},\ -\frac{\mu_{N,t}}{\mu_{D,t}^2}\right).\]

Applying the delta method, \(\text{Var}(\hat{R}_t) \approx \nabla g_t^\top (c_t \Sigma_\varepsilon) \nabla g_t\), it expands to

\[\boxed{\text{Var}(\hat{R}_t) \approx \frac{c_t}{\mu_{D,t}^2}\left(\sigma_N^2 - 2 R_t \sigma_{ND} + R_t^2 \sigma_D^2\right)}\]

The covariance term, \(-2 R_t \sigma_{ND}\), is critical. When the numerator and denominator are positively correlated, ignoring their covariance overstates the variance of the ratio, sometimes by a large factor.

Sanity check the base case

As a sanity check, consider the case with no covariates, so \(Z = [1, T]\). The baseline vectors are \(b_T = (1, 1)\) and \(b_C = (1, 0)\), and

\[Z^\top Z = \begin{pmatrix} n & n_T \\ n_T & n_T \end{pmatrix}, \qquad (Z^\top Z)^{-1} = \frac{1}{n_T n_C} \begin{pmatrix} n_T & -n_T \\ -n_T & n \end{pmatrix}.\]

Then \(c_T = \frac{n_T - 2 n_T + n}{n_T n_C} = \frac{1}{n_T}\) and \(c_C = \frac{n_T}{n_T n_C} = \frac{1}{n_C}\). \(c_t\) now plays the role of the very familiar \(1/n_t\) multiplier in the variance of a mean. The formula for \(\text{Var}(\hat{R}_T)\) reduces to the familiar variance of a ratio of means,

\[\text{Var}(\hat{R}_T) \approx \frac{1}{n_T \mu_{D,T}^2}\left(\sigma_N^2 - 2 R_T \sigma_{ND} + R_T^2 \sigma_D^2\right).\]

When there are covariates, \(c_T\) is no longer exactly \(1/n_T\), but when \(T\) is randomized it remains very close to it. We explain this in the appendix.

Step 2: The Difference of Ratios

The treatment effect is \(\tau = R_T - R_C\), so

\[\text{Var}(\hat{\tau}) = \text{Var}(\hat{R}_T) + \text{Var}(\hat{R}_C) - 2 \text{Cov}(\hat{R}_T, \hat{R}_C).\]

The two arms share \(\beta_2\), so in general their ratios are not independent. The covariance between the estimated ratios of the two arms is

\[\text{Cov}\left(\begin{pmatrix} \hat{\mu}_{N,T} \\ \hat{\mu}_{D,T} \end{pmatrix}, \begin{pmatrix} \hat{\mu}_{N,C} \\ \hat{\mu}_{D,C} \end{pmatrix}\right) = c_{TC} \, \Sigma_\varepsilon, \qquad c_{TC} = b_T (Z^\top Z)^{-1} b_C^\top,\]

yielding

\[\begin{align} \text{Cov}(\hat{R}_T, \hat{R}_C) \approx c_{TC} \nabla g_T^\top \Sigma_\varepsilon \nabla g_C \\ \boxed{\text{Var}(\hat{\tau}) \approx \text{Var}(\hat{R}_T) + \text{Var}(\hat{R}_C) - 2 c_{TC} \, \nabla g_T^\top \Sigma_\varepsilon \nabla g_C} \end{align}\]

In practice the cross term is tiny. When there are no covariates, \(c_{TC} = 0\) exactly, and the two arms’ ratios are uncorrelated. With covariates, the arms are only linked through the shared estimate \(\hat{\beta}_2\), and \(c_{TC}\) shrinks at the rate \(1/n^2\) while \(c_T\) and \(c_C\) shrink at the rate \(1/n\).

Simulation in R

We now execute the math above on simulated data. Continuing the video streaming example, the numerator \(N\) is the number of rebuffers and the denominator \(D\) is the number of videos played. Each user has \(p = 2\) pretreatment covariates, \(X_1\) = plays in the pre-period and \(X_2\) = rebuffers in the pre-period, and the treatment \(T\) is randomized with probability 0.5. The treatment increases plays, but it increases rebuffers by proportionally more, so rebuffers per play go up.

The noise terms \(\varepsilon_N\) and \(\varepsilon_D\) are positively correlated, with \(\sigma_{ND} = \rho \, \sigma_N \sigma_D\).

Code

simulate_experiment <- function(n) {
  # Pretreatment covariates
  X_1 <- rnorm(n, mean = 20, sd = 5)
  X_2 <- 0.1 * X_1 + rnorm(n, mean = 0, sd = 1)
  X <- cbind(X_1, X_2)

  # Randomized treatment. Note that this masks R's shorthand T for TRUE
  T <- rbinom(n, 1, 0.5)

  # Correlated noise
  eps <- mvrnorm(n, mu = c(0, 0), Sigma = Sigma_eps)

  N <- beta_0["N"] + T * beta_1["N"] + drop(X %*% beta_2[, "N"]) + eps[, 1]
  D <- beta_0["D"] + T * beta_1["D"] + drop(X %*% beta_2[, "D"]) + eps[, 2]
  list(N = N, D = D, T = T, X = X)
}

# True coefficients, one column per metric
beta_0 <- c(N = 0.5, D = 5)
beta_1 <- c(N = 0.3, D = 1)
beta_2 <- cbind(N = c(0.05, 0.6), D = c(0.8, 0))

# True noise covariance
sigma_N <- 1.5
sigma_D <- 4
rho <- 0.6
sigma_ND <- rho * sigma_N * sigma_D
Sigma_eps <- matrix(c(sigma_N^2, sigma_ND,
                      sigma_ND, sigma_D^2), 2, 2,
                    dimnames = list(c("N", "D"), c("N", "D")))

set.seed(100)
n <- 20000
sim <- simulate_experiment(n)
N <- sim$N
D <- sim$D
T <- sim$T
X <- sim$X

Since we know the true coefficients, we can compute the true \(\tau\) by plugging \(E[X] = (20, 2)\) into the definitions of \(\mu_{N,T}\), \(\mu_{D,T}\), \(\mu_{N,C}\) and \(\mu_{D,C}\).

Code

E_X <- c(20, 2)
mu_N_T <- unname(beta_0["N"] + beta_1["N"] + sum(E_X * beta_2[, "N"]))
mu_D_T <- unname(beta_0["D"] + beta_1["D"] + sum(E_X * beta_2[, "D"]))
mu_N_C <- unname(beta_0["N"] + sum(E_X * beta_2[, "N"]))
mu_D_C <- unname(beta_0["D"] + sum(E_X * beta_2[, "D"]))
tau <- mu_N_T / mu_D_T - mu_N_C / mu_D_C
c(R_T = mu_N_T / mu_D_T, R_C = mu_N_C / mu_D_C, tau = tau)

        R_T         R_C         tau 
0.136363636 0.128571429 0.007792208

Fitting the Bivariate Model

In R, passing cbind(N, D) as the outcome fits the bivariate OLS model. The coefficient matrix \(\hat{B}\) has one column per metric, \(\hat{\beta}_N\) and \(\hat{\beta}_D\).

Code

fit <- lm(cbind(N, D) ~ T + X)
B_hat <- coef(fit)
beta_N_hat <- B_hat[, "N"]
beta_D_hat <- B_hat[, "D"]
B_hat

                     N           D
(Intercept) 0.49823977  4.98934922
T           0.32121636  0.99843562
XX_1        0.04858011  0.80016834
XX_2        0.59871457 -0.01325133

Next we estimate \(\Sigma_\varepsilon\) from the residuals, using the formulas for \(\hat{\sigma}_N^2\), \(\hat{\sigma}_D^2\) and \(\hat{\sigma}_{ND}\). The matrix product \(\hat{\varepsilon}^\top \hat{\varepsilon}\) computes all three sums at once.

Code

Z <- model.matrix(fit)
k <- ncol(Z)
eps_hat <- resid(fit)
Sigma_eps_hat <- crossprod(eps_hat) / (n - k)
Sigma_eps_hat

         N         D
N 2.260786  3.604297
D 3.604297 15.966950

The estimate is close to the true \(\Sigma_\varepsilon\), which has \(\sigma_N^2 = 2.25\), \(\sigma_{ND} = 3.6\) and \(\sigma_D^2 = 16\). We can also verify that \(\text{Cov}(\text{vec}(\hat{B})) = \Sigma_\varepsilon \otimes (Z^\top Z)^{-1}\) agrees with the covariance matrix that R reports for the bivariate model.

Code

ZtZ_inv <- solve(crossprod(Z))
cov_vec_B_hat <- kronecker(Sigma_eps_hat, ZtZ_inv)
max(abs(cov_vec_B_hat - vcov(fit)))

[1] 2.494532e-15

Point Estimate

The estimated means use the baseline vectors \(b_T = (1, 1, \bar{X})\) and \(b_C = (1, 0, \bar{X})\), where both arms use the same overall average \(\bar{X}\).

Code

X_bar <- colMeans(X)
b_T <- c(1, 1, X_bar)
b_C <- c(1, 0, X_bar)

mu_N_T_hat <- sum(b_T * beta_N_hat)
mu_D_T_hat <- sum(b_T * beta_D_hat)
mu_N_C_hat <- sum(b_C * beta_N_hat)
mu_D_C_hat <- sum(b_C * beta_D_hat)

R_T_hat <- mu_N_T_hat / mu_D_T_hat
R_C_hat <- mu_N_C_hat / mu_D_C_hat
tau_hat <- R_T_hat - R_C_hat
c(R_T_hat = R_T_hat, R_C_hat = R_C_hat, tau_hat = tau_hat)

    R_T_hat     R_C_hat     tau_hat 
0.135977298 0.127125333 0.008851965

Step 1: Variance of the Ratio Within Each Arm

First we compute \(c_T = b_T (Z^\top Z)^{-1} b_T^\top\) and \(c_C = b_C (Z^\top Z)^{-1} b_C^\top\). Since \(T\) is randomized, they are very close to \(1/n_T\) and \(1/n_C\).

Code

n_T <- sum(T)
n_C <- n - n_T
c_T <- drop(t(b_T) %*% ZtZ_inv %*% b_T)
c_C <- drop(t(b_C) %*% ZtZ_inv %*% b_C)
c(c_T = c_T, one_over_n_T = 1 / n_T, c_C = c_C, one_over_n_C = 1 / n_C)

         c_T one_over_n_T          c_C one_over_n_C 
1.001907e-04 1.001904e-04 9.981071e-05 9.981036e-05

Then we apply the delta method, \(\text{Var}(\hat{R}_t) \approx \nabla g_t^\top (c_t \Sigma_\varepsilon) \nabla g_t\), using the gradient \(\nabla g_t = \left(\frac{1}{\mu_{D,t}}, -\frac{\mu_{N,t}}{\mu_{D,t}^2}\right)\).

Code

grad_g_T <- c(1 / mu_D_T_hat, -mu_N_T_hat / mu_D_T_hat^2)
grad_g_C <- c(1 / mu_D_C_hat, -mu_N_C_hat / mu_D_C_hat^2)

var_R_T_hat <- drop(t(grad_g_T) %*% (c_T * Sigma_eps_hat) %*% grad_g_T)
var_R_C_hat <- drop(t(grad_g_C) %*% (c_C * Sigma_eps_hat) %*% grad_g_C)
c(var_R_T_hat = var_R_T_hat, var_R_C_hat = var_R_C_hat)

 var_R_T_hat  var_R_C_hat 
3.277250e-07 3.643962e-07

The expanded form of the boxed formula, \(\frac{c_t}{\mu_{D,t}^2}\left(\sigma_N^2 - 2 R_t \sigma_{ND} + R_t^2 \sigma_D^2\right)\), gives the same answer.

Code

sigma_N_sq_hat <- Sigma_eps_hat["N", "N"]
sigma_D_sq_hat <- Sigma_eps_hat["D", "D"]
sigma_ND_hat <- Sigma_eps_hat["N", "D"]

c_T / mu_D_T_hat^2 *
  (sigma_N_sq_hat - 2 * R_T_hat * sigma_ND_hat + R_T_hat^2 * sigma_D_sq_hat)

[1] 3.27725e-07

Step 2: Variance of the Difference of Ratios

Finally we compute the cross term using \(c_{TC} = b_T (Z^\top Z)^{-1} b_C^\top\), and combine everything into \(\text{Var}(\hat{\tau})\).

Code

c_TC <- drop(t(b_T) %*% ZtZ_inv %*% b_C)
cov_R_T_R_C_hat <- c_TC * drop(t(grad_g_T) %*% Sigma_eps_hat %*% grad_g_C)

var_tau_hat <- var_R_T_hat + var_R_C_hat - 2 * cov_R_T_R_C_hat
se_tau_hat <- sqrt(var_tau_hat)
c(c_TC = c_TC, cov_R_T_R_C_hat = cov_R_T_R_C_hat,
  var_tau_hat = var_tau_hat, se_tau_hat = se_tau_hat)

           c_TC cov_R_T_R_C_hat     var_tau_hat      se_tau_hat 
  -3.552744e-10   -1.227291e-12    6.921236e-07    8.319396e-04

As expected, \(c_{TC}\) is several orders of magnitude smaller than \(c_T\) and \(c_C\), so the cross term barely moves the variance. With the standard error in hand, we can report a confidence interval and p value.

Code

c(tau_hat = tau_hat,
  lower = tau_hat - qnorm(0.975) * se_tau_hat,
  upper = tau_hat + qnorm(0.975) * se_tau_hat,
  p_value = 2 * pnorm(-abs(tau_hat / se_tau_hat)))

     tau_hat        lower        upper      p_value 
8.851965e-03 7.221393e-03 1.048254e-02 1.938121e-26

We said that the two step calculation is identical to applying the delta method once to \(\tau\) as a function of the full coefficient vector \(\text{vec}(\hat{B})\). The gradient of \(\tau = \frac{b_T \beta_N}{b_T \beta_D} - \frac{b_C \beta_N}{b_C \beta_D}\) with respect to \((\beta_N, \beta_D)\) is \(\left(\frac{b_T}{\mu_{D,T}} - \frac{b_C}{\mu_{D,C}}, \ -\frac{\mu_{N,T}}{\mu_{D,T}^2} b_T + \frac{\mu_{N,C}}{\mu_{D,C}^2} b_C\right)\), and it gives the same variance.

Code

grad_tau <- c(b_T / mu_D_T_hat - b_C / mu_D_C_hat,
              -mu_N_T_hat / mu_D_T_hat^2 * b_T + mu_N_C_hat / mu_D_C_hat^2 * b_C)
c(one_shot = drop(t(grad_tau) %*% cov_vec_B_hat %*% grad_tau),
  two_step = var_tau_hat)

    one_shot     two_step 
6.921236e-07 6.921236e-07

Checking the Variance Over Many Experiments

A single experiment cannot tell us whether \(\widehat{\text{Var}}(\hat{\tau})\) is correct. To check it, we repeat the experiment 1,000 times and compare the standard deviation of \(\hat{\tau}\) across experiments to the average estimated standard error. We also compute a naive standard error that ignores the covariance between the numerator and the denominator, by setting \(\hat{\sigma}_{ND} = 0\), which is what we would get from two separate regressions.

Code

estimate_tau <- function(N, D, T, X) {
  fit <- lm(cbind(N, D) ~ T + X)
  B_hat <- coef(fit)
  Z <- model.matrix(fit)
  ZtZ_inv <- solve(crossprod(Z))
  Sigma_eps_hat <- crossprod(resid(fit)) / (nrow(Z) - ncol(Z))

  X_bar <- colMeans(X)
  b_T <- c(1, 1, X_bar)
  b_C <- c(1, 0, X_bar)
  mu_T_hat <- drop(b_T %*% B_hat)  # (mu_N_T_hat, mu_D_T_hat)
  mu_C_hat <- drop(b_C %*% B_hat)  # (mu_N_C_hat, mu_D_C_hat)

  grad_g_T <- c(1 / mu_T_hat["D"], -mu_T_hat["N"] / mu_T_hat["D"]^2)
  grad_g_C <- c(1 / mu_C_hat["D"], -mu_C_hat["N"] / mu_C_hat["D"]^2)
  c_T <- drop(t(b_T) %*% ZtZ_inv %*% b_T)
  c_C <- drop(t(b_C) %*% ZtZ_inv %*% b_C)
  c_TC <- drop(t(b_T) %*% ZtZ_inv %*% b_C)

  var_tau <- function(Sigma) {
    c_T * drop(t(grad_g_T) %*% Sigma %*% grad_g_T) +
      c_C * drop(t(grad_g_C) %*% Sigma %*% grad_g_C) -
      2 * c_TC * drop(t(grad_g_T) %*% Sigma %*% grad_g_C)
  }

  c(tau_hat = unname(mu_T_hat["N"] / mu_T_hat["D"] - mu_C_hat["N"] / mu_C_hat["D"]),
    se_tau_hat = sqrt(var_tau(Sigma_eps_hat)),
    se_tau_hat_naive = sqrt(var_tau(diag(diag(Sigma_eps_hat)))))
}

sims <- t(replicate(1000, {
  sim <- simulate_experiment(n)
  estimate_tau(sim$N, sim$D, sim$T, sim$X)
}))

z <- qnorm(0.975)
c(tau = tau,
  mean_tau_hat = mean(sims[, "tau_hat"]),
  sd_tau_hat = sd(sims[, "tau_hat"]),
  mean_se_tau_hat = mean(sims[, "se_tau_hat"]),
  mean_se_tau_hat_naive = mean(sims[, "se_tau_hat_naive"]),
  coverage = mean(abs(sims[, "tau_hat"] - tau) <= z * sims[, "se_tau_hat"]),
  coverage_naive = mean(abs(sims[, "tau_hat"] - tau) <= z * sims[, "se_tau_hat_naive"]))

                  tau          mean_tau_hat            sd_tau_hat 
          0.007792208           0.007762030           0.000809931 
      mean_se_tau_hat mean_se_tau_hat_naive              coverage 
          0.000827188           0.001047469           0.951000000 
       coverage_naive 
          0.985000000

The average of \(\hat{\tau}\) is close to the true \(\tau\), and the average standard error from the delta method matches the actual standard deviation of \(\hat{\tau}\) across experiments, so the 95% confidence interval covers the true \(\tau\) about 95% of the time. The naive standard error, which ignores \(\sigma_{ND}\), is too large, and its confidence interval is overly conservative. When the numerator and denominator are positively correlated, jointly modeling them gives a more precise, and still correct, measure of uncertainty.

Appendix

How Covariates Change \(c_t\)

To see what covariates do to \(c_T\), we return to a single outcome \(y = \beta_0 + T \beta_1 + X \beta_2 + \varepsilon\) with \(\text{Var}(\varepsilon) = \sigma^2\). We will find a different way to express \(c_T\), as a sum of \(1/n_T\) and another term that is a function of covariate imbalance.

Because the model contains an intercept and \(T\), the residuals sum to zero within each arm, so the fitted model passes through each arm’s own averages. For the treatment arm, the raw average of \(y\) equals the fitted value at the arm’s own covariate average, \(\bar{y}_T = \hat{\beta}_0 + \hat{\beta}_1 + \bar{X}_T \hat{\beta}_2\), where \(\bar{y}_T\) and \(\bar{X}_T\) are the averages of \(y\) and \(X\) over treated users. Subtracting this from \(\hat{\mu}_T = b_T \hat{\beta} = \hat{\beta}_0 + \hat{\beta}_1 + \bar{X} \hat{\beta}_2\) gives

\[\hat{\mu}_T = \bar{y}_T + \Delta_T \hat{\beta}_2, \qquad \Delta_T = \bar{X} - \bar{X}_T.\]

The estimated treatment mean is the raw treatment average, adjusted for how far the treatment arm’s covariates are from the overall average \(\bar{X}\). Note that under randomization \(\bar{X}_T\) is close to \(\bar{X}\), and the gap shrinks as the sample grows, but in any finite sample they differ.

Since \(\hat{\mu}_T\) is a sum of two random quantities, \(\bar{y}_T\) and \(\hat{\beta}_2\), its variance in general has three terms

\[\text{Var}(\hat{\mu}_T) = \text{Var}(\bar{y}_T) + \Delta_T \, \text{Var}(\hat{\beta}_2) \, \Delta_T^\top + 2 \Delta_T \, \text{Cov}(\hat{\beta}_2, \bar{y}_T).\]

The first two terms are familiar. The average of \(n_T\) independent outcomes has \(\text{Var}(\bar{y}_T) = \sigma^2 / n_T\). For the second term, the Frisch-Waugh-Lovell theorem says that \(\hat{\beta}_2\) is the coefficient from regressing \(y\) on \(\tilde{X}\), the covariates demeaned within each arm, so \(\hat{\beta}_2 = (\tilde{X}^\top \tilde{X})^{-1} \tilde{X}^\top y\) and \(\text{Var}(\hat{\beta}_2) = \sigma^2 (\tilde{X}^\top \tilde{X})^{-1}\).

The third term, the cross term, is exactly zero. Let \(1_T\) be the vector that is 1 for treated users and 0 otherwise, so that \(\bar{y}_T = \frac{1}{n_T} 1_T^\top y\). Both \(\bar{y}_T\) and \(\hat{\beta}_2\) are linear in \(y\), and \(\text{Cov}(y) = \sigma^2 I\), so

\[\text{Cov}(\hat{\beta}_2, \bar{y}_T) = (\tilde{X}^\top \tilde{X})^{-1} \tilde{X}^\top \, \sigma^2 I \, \frac{1_T}{n_T} = \frac{\sigma^2}{n_T} (\tilde{X}^\top \tilde{X})^{-1} \tilde{X}^\top 1_T.\]

The product \(\tilde{X}^\top 1_T\) is the sum of the demeaned covariates over the treated users. Since each treated user’s covariates were demeaned by the treatment arm’s own average, \(\tilde{X}^\top 1_T = \sum_{i \in T} (X_i - \bar{X}_T)^\top = 0\). This is why the demeaning must be done within each arm. If we demeaned by the overall average \(\bar{X}\) instead, the sum would be \(n_T (\bar{X}_T - \bar{X})^\top\), which is not zero.

With the cross term gone, the variance is the sum of the first two terms

\[\text{Var}(\hat{\mu}_T) = \frac{\sigma^2}{n_T} + \Delta_T \, \sigma^2 (\tilde{X}^\top \tilde{X})^{-1} \Delta_T^\top.\]

Earlier we showed that \(\text{Var}(\hat{\mu}_T) = \sigma^2 c_T\). Both expressions are the variance of the same estimator, so setting them equal and canceling \(\sigma^2\) gives an alternative formula for \(c_T\),

\[\boxed{c_T = \frac{1}{n_T} + \Delta_T (\tilde{X}^\top \tilde{X})^{-1} \Delta_T^\top}\]

The second term is a quadratic form, so it is never negative: \(c_T \geq 1/n_T\), with equality only when the treatment arm’s covariates average exactly to \(\bar{X}\). It is a penalty for covariate imbalance, and it measures how far the model has to extrapolate the treatment mean from \(\bar{X}_T\) to \(\bar{X}\). The same argument applies to \(c_C\) with \(\Delta_C = \bar{X} - \bar{X}_C\).

Under randomization the penalty is very small, for two reasons:

\(\Delta_T = \frac{n_C}{n}(\bar{X}_C - \bar{X}_T)\) is a difference between two sample means that have the same expected value. It is centered at zero, and its magnitude shrinks at the rate \(1/\sqrt{n}\).
\(\tilde{X}^\top \tilde{X}\) is a sum over all \(n\) users, roughly \(n \, \text{Var}(X)\), so its inverse shrinks at the rate \(1/n\).

Together the penalty shrinks at the rate \(1/n^2\), while \(1/n_T\) shrinks at the rate \(1/n\). Taking expectations, the relative excess \(n_T c_T - 1\) is roughly \(p \, n_C / n^2\) when \(X\) has \(p\) columns, which is negligible for any reasonably sized experiment.

This breaks down without randomization. If treated users have systematically different covariates, then \(\Delta_T\) does not shrink, the penalty is on the same order as \(1/n_T\), and the variance of \(\hat{\mu}_T\) can be much larger than \(\sigma^2 / n_T\).

Finally, this does not mean covariates are harmful. The variance of each mean is \(c_t\) times the noise variance. Adding covariates increases \(c_t\) by a negligible amount, but can greatly reduce the residual variances \(\sigma_N^2\) and \(\sigma_D^2\), and with them the variance of the ratio. This is the same tradeoff that makes CUPED effective.