A/B Testing — Statistical Hypothesis Testing

import pandas as pd
import numpy as np
from scipy import stats

# ── Full experiment dataset — 6 weeks, 2 variants, 2 devices, 2 sources ──────
rows = []
data = [
    # week, variant, users, conv, device, source, revenue
    (1,"control",  1621,186,"desktop","organic",9300),
    (1,"control",  1597,184,"mobile", "organic",9200),
    (1,"treatment",1628,212,"desktop","organic",10600),
    (1,"treatment",1613,209,"mobile", "organic",10450),
    (2,"control",  1594,181,"desktop","paid",   9050),
    (2,"control",  1590,181,"mobile", "paid",   9050),
    (2,"treatment",1608,215,"desktop","paid",   10750),
    (2,"treatment",1598,213,"mobile", "paid",   10650),
    (3,"control",  1658,191,"desktop","organic",9550),
    (3,"control",  1643,190,"mobile", "organic",9500),
    (3,"treatment",1671,221,"desktop","organic",11050),
    (3,"treatment",1659,218,"mobile", "organic",10900),
    (4,"control",  1612,184,"desktop","paid",   9200),
    (4,"control",  1604,183,"mobile", "paid",   9150),
    (4,"treatment",1624,218,"desktop","paid",   10900),
    (4,"treatment",1618,215,"mobile", "paid",   10750),
    (5,"control",  1588,181,"desktop","organic",9050),
    (5,"control",  1576,179,"mobile", "organic",8950),
    (5,"treatment",1601,213,"desktop","organic",10650),
    (5,"treatment",1593,210,"mobile", "organic",10500),
    (6,"control",  1634,186,"desktop","paid",   9300),
    (6,"control",  1618,184,"mobile", "paid",   9200),
    (6,"treatment",1647,221,"desktop","paid",   11050),
    (6,"treatment",1638,218,"mobile", "paid",   10900),
]

df = pd.DataFrame(data, columns=["week","variant","users","conversions",
                                  "device","source","revenue_gbp"])

# ── Overall totals per variant ────────────────────────────────────────────────
overall = (df.groupby("variant")
             .agg(total_users      =("users",       "sum"),
                  total_conversions=("conversions",  "sum"),
                  total_revenue    =("revenue_gbp",  "sum"))
             .reset_index())

overall["conv_rate_pct"] = (overall["total_conversions"]
                            / overall["total_users"] * 100).round(4)
overall["avg_order_gbp"] = (overall["total_revenue"]
                            / overall["total_conversions"]).round(2)

print("Overall experiment summary:")
print(overall.to_string(index=False))

# ── Extract scalars for testing ───────────────────────────────────────────────
ctrl  = overall[overall["variant"]=="control"].iloc[0]
treat = overall[overall["variant"]=="treatment"].iloc[0]

n_c, x_c = int(ctrl["total_users"]),  int(ctrl["total_conversions"])
n_t, x_t = int(treat["total_users"]), int(treat["total_conversions"])
p_c = x_c / n_c
p_t = x_t / n_t

print(f"\nControl:   n={n_c:,}  conversions={x_c:,}  rate={p_c*100:.4f}%")
print(f"Treatment: n={n_t:,}  conversions={x_t:,}  rate={p_t*100:.4f}%")
print(f"Absolute difference: {(p_t - p_c)*100:+.4f} pp")
print(f"Relative uplift:     {(p_t - p_c)/p_c*100:+.2f}%")

INTENDED_SPLIT = 0.50   # intended fraction in each group

total_users    = n_c + n_t
expected_c     = total_users * INTENDED_SPLIT
expected_t     = total_users * INTENDED_SPLIT
actual_split_c = n_c / total_users
actual_split_t = n_t / total_users

# ── Chi-square goodness-of-fit for SRM ───────────────────────────────────────
# Observed counts vs expected counts under 50/50 split
observed  = np.array([n_c, n_t])
expected  = np.array([expected_c, expected_t])

srm_chi2, srm_p = stats.chisquare(f_obs=observed, f_exp=expected)

print("Sample Ratio Mismatch (SRM) Check:")
print(f"  Total users:          {total_users:,}")
print(f"  Control:   {n_c:,}  ({actual_split_c*100:.2f}%  expected 50.00%)")
print(f"  Treatment: {n_t:,}  ({actual_split_t*100:.2f}%  expected 50.00%)")
print(f"\n  Chi-square statistic: {srm_chi2:.4f}")
print(f"  p-value:              {srm_p:.4f}")

if srm_p < 0.01:
    verdict = "FAIL — SRM detected. Experiment result is INVALID."
elif srm_p < 0.05:
    verdict = "WARNING — borderline SRM. Investigate assignment pipeline."
else:
    verdict = "PASS — no significant imbalance. Randomisation looks clean."

print(f"\n  SRM verdict: {verdict}")
print(f"\n  Imbalance: {abs(n_t - n_c)} more users in treatment "
      f"({abs(actual_split_t - actual_split_c)*100:.3f}pp off target)")

ALPHA = 0.05    # significance threshold

# ── Two-proportion Z-test ─────────────────────────────────────────────────────
# Pooled proportion under null hypothesis (no difference)
p_pool = (x_c + x_t) / (n_c + n_t)

# Standard error of the difference under H0
se_pool = np.sqrt(p_pool * (1 - p_pool) * (1/n_c + 1/n_t))

# Z-statistic
z_stat = (p_t - p_c) / se_pool

# Two-tailed p-value
p_val_z = 2 * (1 - stats.norm.cdf(abs(z_stat)))

print("Two-proportion Z-Test:")
print(f"  Pooled proportion (H0): {p_pool*100:.4f}%")
print(f"  Standard error:         {se_pool*100:.5f} pp")
print(f"  Z-statistic:            {z_stat:.4f}")
print(f"  p-value (two-tailed):   {p_val_z:.6f}")
print(f"  Significant at α={ALPHA}: {'YES' if p_val_z < ALPHA else 'NO'}")

# ── Chi-square test of independence on 2×2 contingency table ──────────────────
# Rows: control / treatment  |  Cols: converted / not converted
contingency = np.array([
    [x_c, n_c - x_c],    # control:   converted, not converted
    [x_t, n_t - x_t],    # treatment: converted, not converted
])

chi2, p_val_chi2, dof, expected_freq = stats.chi2_contingency(contingency)

print(f"\nChi-Square Test of Independence:")
print(f"  Contingency table:")
print(f"              Converted   Not Converted")
print(f"  Control:    {x_c:>9,}   {n_c-x_c:>13,}")
print(f"  Treatment:  {x_t:>9,}   {n_t-x_t:>13,}")
print(f"\n  Chi2 statistic: {chi2:.4f}")
print(f"  Degrees of freedom: {dof}")
print(f"  p-value:            {p_val_chi2:.6f}")
print(f"  Significant at α={ALPHA}: {'YES' if p_val_chi2 < ALPHA else 'NO'}")

print(f"\nBoth tests agree: {'YES' if (p_val_z < ALPHA) == (p_val_chi2 < ALPHA) else 'NO — INVESTIGATE'}")

CONFIDENCE = 0.95
WEEKLY_USERS = 40000          # total platform weekly active users at 100% rollout
AVG_ORDER_GBP = 50.00         # average order value
WEEKS_PER_YEAR = 52

# ── Unpooled SE for confidence interval (not the pooled SE used in the test) ──
se_unpooled = np.sqrt(p_c*(1-p_c)/n_c + p_t*(1-p_t)/n_t)

# ── Z critical value for desired confidence level ─────────────────────────────
z_crit = stats.norm.ppf(1 - (1 - CONFIDENCE) / 2)

point_est  = p_t - p_c
ci_lower   = point_est - z_crit * se_unpooled
ci_upper   = point_est + z_crit * se_unpooled

print(f"95% Confidence Interval for Absolute Uplift:")
print(f"  Point estimate:  {point_est*100:+.4f} pp")
print(f"  CI lower bound:  {ci_lower*100:+.4f} pp")
print(f"  CI upper bound:  {ci_upper*100:+.4f} pp")
print(f"  Interval:        [{ci_lower*100:+.2f} pp,  {ci_upper*100:+.2f} pp]")
print(f"\n  Relative uplift: {point_est/p_c*100:+.2f}%  "
      f"(CI: [{ci_lower/p_c*100:+.1f}%,  {ci_upper/p_c*100:+.1f}%])")

# ── Revenue impact at 100% rollout ────────────────────────────────────────────
def annual_revenue_uplift(uplift_pp):
    extra_conv_weekly = WEEKLY_USERS * uplift_pp
    return extra_conv_weekly * AVG_ORDER_GBP * WEEKS_PER_YEAR

rev_point     = annual_revenue_uplift(point_est)
rev_lower     = annual_revenue_uplift(ci_lower)
rev_upper     = annual_revenue_uplift(ci_upper)

print(f"\nAnnualised Revenue Uplift at 100% Rollout ({WEEKLY_USERS:,} users/week, £{AVG_ORDER_GBP} AOV):")
print(f"  Conservative (CI lower): £{rev_lower:,.0f}/year")
print(f"  Point estimate:          £{rev_point:,.0f}/year")
print(f"  Optimistic (CI upper):   £{rev_upper:,.0f}/year")
print(f"\n  Recommendation: use £{rev_lower:,.0f} (lower bound) for business case — "
      f"never the point estimate.")

ALPHA   = 0.05    # type I error rate
POWER   = 0.80    # desired power (1 - type II error rate)
EFFECT  = p_t - p_c   # observed effect — also the MDE we'd want to detect

# ── Z critical values ─────────────────────────────────────────────────────────
z_alpha = stats.norm.ppf(1 - ALPHA / 2)   # two-tailed: 1.96
z_beta  = stats.norm.ppf(POWER)            # 80% power: 0.842

# ── Minimum sample size per group ────────────────────────────────────────────
# Standard formula: n = (z_alpha + z_beta)^2 * 2 * p_bar * (1 - p_bar) / delta^2
p_bar = (p_c + p_t) / 2   # average proportion for sample size formula
delta = EFFECT             # minimum detectable effect = observed effect

n_required = ((z_alpha + z_beta)**2 * 2 * p_bar * (1 - p_bar)) / (delta**2)
n_required_ceil = int(np.ceil(n_required))

print(f"Power Analysis:")
print(f"  Baseline conversion rate (control):  {p_c*100:.4f}%")
print(f"  Target conversion rate (treatment):  {p_t*100:.4f}%")
print(f"  Minimum detectable effect (MDE):     {delta*100:.4f} pp")
print(f"  Desired power:                        {POWER*100:.0f}%")
print(f"  α (two-tailed):                       {ALPHA}")
print(f"\n  Z_α/2 (critical value):  {z_alpha:.4f}")
print(f"  Z_β   (power):           {z_beta:.4f}")
print(f"\n  Minimum sample per group: {n_required_ceil:,}")
print(f"  Minimum sample total:     {n_required_ceil*2:,}")
print(f"\n  Actual sample per group:  ~{(n_c+n_t)//2:,}")
print(f"  Actual sample total:      {n_c+n_t:,}")
print(f"  Met requirement:          {'YES' if min(n_c,n_t) >= n_required_ceil else 'NO'}")

# ── Achieved power with actual sample sizes ───────────────────────────────────
se_actual = np.sqrt(p_bar * (1 - p_bar) * 2 / ((n_c + n_t) / 2))
z_achieved = delta / se_actual - z_alpha
achieved_power = stats.norm.cdf(z_achieved)

print(f"\n  Achieved power with actual n: {achieved_power*100:.1f}%")

# ── How small an effect could we detect at 80% power with actual n? ──────────
mde_actual = (z_alpha + z_beta) * np.sqrt(2 * p_bar * (1-p_bar) / ((n_c+n_t)/2))
print(f"  Actual MDE at 80% power:      {mde_actual*100:.4f} pp "
      f"({mde_actual/p_c*100:.1f}% relative)")

# ── Segment breakdown: device and source ─────────────────────────────────────
for seg_col in ["device", "source"]:
    print(f"\nSegment analysis by {seg_col}:")
    print(f"{'Segment':<12} {'Variant':<12} {'Users':>7} {'Conv':>6} {'Rate':>8} {'Uplift':>9} {'p-val':>8}")
    print("─" * 62)

    for seg_val in df[seg_col].unique():
        seg_df = df[df[seg_col]==seg_val]
        agg = (seg_df.groupby("variant")
                     .agg(users=("users","sum"), conv=("conversions","sum"))
                     .reset_index())

        sc = agg[agg["variant"]=="control"].iloc[0]
        st = agg[agg["variant"]=="treatment"].iloc[0]

        pc_s = sc["conv"] / sc["users"]
        pt_s = st["conv"] / st["users"]
        uplift_s = (pt_s - pc_s) * 100

        # Z-test for this segment
        pp_s  = (sc["conv"]+st["conv"]) / (sc["users"]+st["users"])
        se_s  = np.sqrt(pp_s*(1-pp_s)*(1/sc["users"]+1/st["users"]))
        z_s   = (pt_s - pc_s) / se_s
        pv_s  = 2*(1-stats.norm.cdf(abs(z_s)))

        for row, v in [(sc,"control"),(st,"treatment")]:
            rate = row["conv"]/row["users"]*100
            upl  = f"{uplift_s:+.2f}pp" if v=="treatment" else "—"
            pv   = f"{pv_s:.4f}" if v=="treatment" else "—"
            print(f"  {seg_val:<10} {v:<12} {int(row['users']):>7,} "
                  f"{int(row['conv']):>6,} {rate:>7.2f}% {upl:>9} {pv:>8}")

# ── Novelty effect check: early vs late weeks ─────────────────────────────────
print(f"\nNovelty effect check — weekly treatment uplift:")
weekly = (df.groupby(["week","variant"])
            .agg(users=("users","sum"), conv=("conversions","sum"))
            .reset_index())

print(f"{'Week':<6} {'Control Rate':>14} {'Treatment Rate':>15} {'Uplift':>8}")
print("─" * 46)
for wk in range(1, 7):
    wdf = weekly[weekly["week"]==wk]
    wc  = wdf[wdf["variant"]=="control"].iloc[0]
    wt  = wdf[wdf["variant"]=="treatment"].iloc[0]
    pc_w = wc["conv"]/wc["users"]*100
    pt_w = wt["conv"]/wt["users"]*100
    print(f"  {wk:<4}   {pc_w:>12.2f}%   {pt_w:>13.2f}%   {pt_w-pc_w:>+7.2f}pp")