Analysing Experiments with Regression

24.1 A difference in means is a regression (40 min)

Concept. Regressing the outcome on a 0/1 treatment dummy, revenue ~ treated, gives:

• intercept = mean of the control group,
• coefficient on treated = treated mean − control mean.

So the regression is the difference in means. Its value comes from what it adds: standard errors, covariates, interactions, clustering, all in one familiar framework.

Create day24/ols_equals_diff.py:

# day24/ols_equals_diff.py   (python -m day24.ols_equals_diff)
"""The workhorse: a difference in means IS a regression on a treatment dummy."""
import numpy as np
import statsmodels.formula.api as smf

from causalkit.experiment_data import email_v2_results
from causalkit.lab import truth

df = email_v2_results().to_pandas()

diff = df.loc[df.treated == 1, "revenue"].mean() - df.loc[df.treated == 0, "revenue"].mean()
classical = smf.ols("revenue ~ treated", df).fit()
robust = smf.ols("revenue ~ treated", df).fit(cov_type="HC2")

print(f"Difference in means          {diff:+.3f}")
print(f"OLS coefficient on treated   {classical.params['treated']:+.3f}")
print(f"Classical SE {classical.bse['treated']:.3f}   Robust (HC2) SE {robust.bse['treated']:.3f}")
lo, hi = robust.conf_int().loc["treated"]
print(f"95% CI [{lo:+.3f}, {hi:+.3f}]   p = {robust.pvalues['treated']:.4f}")

# Relative lift and its CI (delta method for a ratio of two means)
t = df.loc[df.treated == 1, "revenue"]; c = df.loc[df.treated == 0, "revenue"]
ratio = t.mean() / c.mean()
se_log = np.sqrt(t.var() / (len(t) * t.mean() ** 2) + c.var() / (len(c) * c.mean() ** 2))
print(f"Relative lift {ratio - 1:+.2%}   95% CI [{ratio * np.exp(-1.96 * se_log) - 1:+.2%}, "
      f"{ratio * np.exp(1.96 * se_log) - 1:+.2%}]")
print(f"Truth: {truth('email_v2')['revenue']:+.3f} per customer, {truth('email_v2')['revenue_relative']:+.1%}")

python -m day24.ols_equals_diff

Difference in means          +2.446
OLS coefficient on treated   +2.446
Classical SE 1.011   Robust (HC2) SE 1.011
95% CI [+0.464, +4.428]   p = 0.0156
Relative lift +3.49%   95% CI [+0.65%, +6.41%]
Truth: +3.480 per customer, +5.0%

This one experiment estimates +$2.45, below the true +$3.48 but well inside its confidence interval [0.46, 4.43]. That's what a single experiment looks like: right on average, off by about one standard error in any given run. Report the interval, not just the point.

24.2 Robust standard errors and relative lifts (50 min)

Robust (HC2) standard errors don't assume both groups have the same variance. With equal group sizes they're almost identical to classical ones (as above). With unequal groups and a treatment that changes the variance, classical standard errors can be badly wrong. Use robust ones by default; it costs nothing.

Relative lift ("+3.5%") is how results are usually communicated. Its confidence interval needs care, because it's a ratio of two noisy means. The script uses the delta method on the log of the ratio, which gives an asymmetric interval [+0.65%, +6.41%].

Try it. Treat only 10% of customers instead of half. Compare classical and robust standard errors.

# day24/unequal_split.py   (python -m day24.unequal_split)
import numpy as np
import statsmodels.formula.api as smf

from causalkit.lab import population, run

pop = population()
rng = np.random.default_rng(241)
z = (rng.random(pop.height) < 0.10).astype(np.int8)          # only 10% treated
df = run(pop, z, "email_v2", seed=241).to_pandas()
for cov in ["nonrobust", "HC2"]:
    fit = smf.ols("revenue ~ treated", df).fit(cov_type=cov)
    print(f"{cov:<10} effect {fit.params['treated']:+.3f}  SE {fit.bse['treated']:.3f}")

nonrobust  effect +6.806  SE 1.635
HC2        effect +6.806  SE 1.847

24.3 Covariates: same target, less noise (60 min)

Concept. Adding pre-treatment covariates to the regression doesn't change what is estimated in a randomised experiment (randomisation already removed bias). It removes the part of the outcome's variation that the covariates explain, so the estimate is less noisy. Two common specifications:

• ANCOVA: revenue ~ treated + pre_revenue + loyalty.
• Lin (2013): interact treatment with centred covariates, revenue ~ treated * (pre_c + loy_c). Lin proved this version never makes precision worse, even if the effect differs across customers.

Create day24/covariate_adjustment.py. Each replay draws a fresh complete randomisation (not stratified), so covariates have chance imbalances to correct:

# day24/covariate_adjustment.py   (python -m day24.covariate_adjustment)
"""
Adding pre-treatment covariates: same target, less noise. Three specifications, 300 replays,
each with a fresh complete randomisation (no stratification, so covariates have imbalance to fix):
  plain     revenue ~ treated
  ANCOVA    revenue ~ treated + pre_revenue + loyalty
  Lin       revenue ~ treated * (centred covariates)   (Lin 2013: never hurts precision)
"""
import numpy as np
import polars as pl
import statsmodels.formula.api as smf

from causalkit.lab import population, run, truth

pop = population()
rng = np.random.default_rng(24)
cov = pop.select("customer_id", "loyalty", "pre_revenue").with_columns(
    pre_c=pl.col("pre_revenue") - pl.col("pre_revenue").mean(),
    loy_c=pl.col("loyalty") - pl.col("loyalty").mean())

SPECS = {"plain": "revenue ~ treated",
         "ANCOVA": "revenue ~ treated + pre_revenue + loyalty",
         "Lin": "revenue ~ treated * (pre_c + loy_c)"}
est = {k: [] for k in SPECS}
covered = {k: 0 for k in SPECS}
target = truth("email_v2")["revenue"]
for r in range(300):
    z = rng.permutation(np.repeat([0, 1], pop.height // 2))
    df = run(pop, z, "email_v2", seed=1000 + r).join(cov, on="customer_id").to_pandas()
    for name, f in SPECS.items():
        fit = smf.ols(f, df).fit(cov_type="HC2")
        est[name].append(fit.params["treated"])
        lo, hi = fit.conf_int().loc["treated"]
        covered[name] += lo <= target <= hi

rows = [{"specification": k, "mean": round(np.mean(v), 3), "SD across replays": round(np.std(v), 3),
         "95% CI coverage": f"{covered[k] / 300:.0%}"} for k, v in est.items()]
print(pl.DataFrame(rows))
print(f"True effect {target:+.3f}")

python -m day24.covariate_adjustment

│ specification ┆ mean  ┆ SD across replays ┆ 95% CI coverage │
│ plain         ┆ 3.481 ┆ 1.086             ┆ 94%             │
│ ANCOVA        ┆ 3.487 ┆ 0.942             ┆ 92%             │
│ Lin           ┆ 3.487 ┆ 0.942             ┆ 92%             │
True effect +3.480

All three are unbiased. Covariates cut the standard deviation by 13%. ANCOVA and Lin agree here because the effect is a constant 5% for everyone. Coverage of 92–94% over 300 replays is within simulation error of 95% (±2.5 points).

Only pre-treatment covariates. A covariate measured after treatment is a bad control (Day 17), even in a perfect experiment.

24.4 Binary outcomes (35 min)

For a 0/1 outcome like conversion, the same regression is called a linear probability model (LPM). Its coefficient is the effect in percentage points, directly. A logit model gives the same answer once converted to an average marginal effect, but its raw coefficient (log-odds) is hard to explain to a manager.

Create day24/binary_outcome.py:

# day24/binary_outcome.py   (python -m day24.binary_outcome)
"""Conversion (0/1): the linear probability model gives the effect in percentage points directly."""
import numpy as np
import statsmodels.formula.api as smf

from causalkit.experiment_data import email_v2_results
from causalkit.lab import truth

df = email_v2_results().to_pandas()

lpm = smf.ols("converted ~ treated", df).fit(cov_type="HC2")
logit = smf.logit("converted ~ treated", df).fit(disp=False)
ame = logit.get_margeff().summary_frame().loc["treated"]          # average marginal effect

print(f"Linear probability model: {lpm.params['treated'] * 100:+.2f} points  (SE {lpm.bse['treated'] * 100:.2f})")
print(f"Logit, average marginal effect: {ame['dy/dx'] * 100:+.2f} points  (SE {ame['Std. Err.'] * 100:.2f})")
print(f"Logit coefficient (log-odds, harder to explain): {logit.params['treated']:+.4f}")
print(f"Truth: {truth('email_v2')['converted'] * 100:+.2f} points")

Linear probability model: +2.42 points  (SE 0.49)
Logit, average marginal effect: +2.42 points  (SE 0.49)
Logit coefficient (log-odds, harder to explain): +0.1019
Truth: +1.51 points

In a randomised experiment, the LPM with robust standard errors is the standard choice. (This run overshoots the truth by 1.9 standard errors: unusual, but it happens about 6% of the time. Don't adjust your method to the truth after seeing it.)

24.5 Day 24 mini-project: a standard readout (70 min)

Goal. Package the analysis into a function every experiment can use: one row per metric with control mean, effect, robust SE, CI, relative lift and p-value.

Create causalkit/analysis.py:

"""
A standard experiment readout (Day 24 mini-project; CUPED added on Day 25).

    from causalkit.analysis import readout
    readout(df, metrics=["revenue", "orders", "converted"], covariates=["pre_revenue", "loyalty"])
"""
import numpy as np
import polars as pl
import statsmodels.formula.api as smf


def readout(df: pl.DataFrame, metrics: list[str], covariates: list[str] | None = None,
            treatment: str = "treated") -> pl.DataFrame:
    """One row per metric: control mean, absolute effect, robust SE, 95% CI, relative lift, p-value."""
    pdf = df.to_pandas()
    rows = []
    for m in metrics:
        rhs = " + ".join([treatment] + (covariates or []))
        fit = smf.ols(f"{m} ~ {rhs}", pdf).fit(cov_type="HC2")
        control_mean = pdf.loc[pdf[treatment] == 0, m].mean()
        eff, se = fit.params[treatment], fit.bse[treatment]
        rows.append({
            "metric": m,
            "control_mean": round(control_mean, 3),
            "effect": round(eff, 3),
            "se": round(se, 3),
            "ci_low": round(eff - 1.96 * se, 3),
            "ci_high": round(eff + 1.96 * se, 3),
            "lift": f"{eff / control_mean:+.2%}",
            "p_value": f"{fit.pvalues[treatment]:.4f}" if fit.pvalues[treatment] >= 0.0001 else "<0.0001",
        })
    return pl.DataFrame(rows)

Then day24/readout.py:

# day24/readout.py   (python -m day24.readout)
"""Mini-project: the email_v2 readout, without and with covariate adjustment, against the truth."""
import polars as pl

from causalkit.analysis import readout
from causalkit.experiment_data import email_v2_results
from causalkit.lab import truth

df = email_v2_results()
METRICS = ["revenue", "orders", "converted", "margin"]

with pl.Config(tbl_width_chars=120):
    print("Unadjusted")
    print(readout(df, METRICS))
    print("Adjusted for pre_revenue, pre_orders and strata (strata already include loyalty)")
    print(readout(df, METRICS, covariates=["pre_revenue", "pre_orders", "C(stratum)"]))
t = truth("email_v2")
print("Truth:", {m: round(t[m], 3) for m in METRICS})

python -m day24.readout

Unadjusted
│ metric    ┆ control_mean ┆ effect ┆ se    ┆ ci_low ┆ ci_high ┆ lift   ┆ p_value │
│ revenue   ┆ 70.095       ┆ 2.446  ┆ 1.011 ┆ 0.464  ┆ 4.428   ┆ +3.49% ┆ 0.0156  │
│ orders    ┆ 1.143        ┆ 0.057  ┆ 0.014 ┆ 0.029  ┆ 0.084   ┆ +4.95% ┆ <0.0001 │
│ converted ┆ 0.598        ┆ 0.024  ┆ 0.005 ┆ 0.015  ┆ 0.034   ┆ +4.05% ┆ <0.0001 │
│ margin    ┆ 21.029       ┆ 0.734  ┆ 0.303 ┆ 0.139  ┆ 1.328   ┆ +3.49% ┆ 0.0156  │
Adjusted for pre_revenue, pre_orders and strata (strata already include loyalty)
│ revenue   ┆ 70.095       ┆ 2.634  ┆ 0.873 ┆ 0.923  ┆ 4.345   ┆ +3.76% ┆ 0.0026  │
│ orders    ┆ 1.143        ┆ 0.053  ┆ 0.012 ┆ 0.029  ┆ 0.077   ┆ +4.61% ┆ <0.0001 │
│ converted ┆ 0.598        ┆ 0.024  ┆ 0.005 ┆ 0.014  ┆ 0.033   ┆ +3.94% ┆ <0.0001 │
│ margin    ┆ 21.029       ┆ 0.79   ┆ 0.262 ┆ 0.277  ┆ 1.304   ┆ +3.76% ┆ 0.0026  │
Truth: {'revenue': 3.48, 'orders': 0.057, 'converted': 0.015, 'margin': 1.044}

Adjustment shrinks every revenue and margin SE by about 14%. Don't add loyalty next to C(stratum): the strata already contain loyalty, so the two are perfectly collinear and the regression can't separate them.

Write a Results section in docs/EXPERIMENT_PLAN_email_v2.md using the adjusted table and the plan's decision rule.

24.6 Save your work with Git (15 min)

git add .
git commit -m "Day 24: regression analysis, robust SEs, covariates, readout function"
git push

Day 24 checklist

• I can show that a difference in means equals a regression on a treatment dummy
• I use robust standard errors by default and can report a relative lift with its CI
• I can explain why pre-treatment covariates reduce noise without changing the target
• I can analyse a binary outcome with a linear probability model
• I have a reusable readout function for every experiment