import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
# statsmodels = the R/SAS-style stats library (pre-installed on Colab): it does the standard errors,
# test statistics, intervals, and post-hoc comparisons for us, so we write far less by-hand code.
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.proportion import proportions_ztest, confint_proportions_2indep, proportion_confint
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.weightstats import DescrStatsW, CompareMeans
from scipy.stats.contingency import association   # Cramer's V in one call
CY="#0891b2"; DEEP="#0e7490"; LIGHT="#67e8f9"; INK="#1a2138"; GRID="#e6e9f2"; GREEN="#059669"; RED="#ef4444"; AMBER="#d97706"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
   "axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
   "axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
pd.set_option("display.width",120)
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
rng = np.random.default_rng(84)

try:    e = pd.read_excel("../../data/case-study-a-web-a-b-test--web_abtest.xlsx", sheet_name="Sessions")
except FileNotFoundError: e = pd.read_excel(BASE+"case-study-a-web-a-b-test--web_abtest.xlsx", sheet_name="Sessions")
print("shape:", e.shape, "| missing values:", e.isna().sum().sum())
display(e.describe(include="all").T[["count","unique","top","mean"]])

shape: (3600, 7) | missing values: 0

# conversion and time-on-page by variant
print(e.groupby("variant").agg(sessions=("converted","size"), conv_rate=("converted","mean"),
                               avg_seconds=("seconds_on_page","mean")).round(3))
# randomization check: device mix should be similar across A and B
print("\ndevice mix by variant (should be balanced):")
print(pd.crosstab(e.variant, e.device, normalize="index").round(3))

         sessions  conv_rate  avg_seconds
variant                                  
A            1830      0.099       54.410
B            1770      0.134       54.048

device mix by variant (should be balanced):
device   desktop  mobile  tablet
variant                         
A          0.309   0.610   0.081
B          0.328   0.596   0.076

print("Decision: proportion outcome + 2 independent groups -> TWO-PROPORTION z-TEST")
print("H0: pB = pA   (the redesign does not change conversion)")
print("H1: pB > pA   (the redesign converts better) -- one-sided")
A=e[e.variant=="A"]; B=e[e.variant=="B"]
cA,nA=A.converted.sum(),len(A); cB,nB=B.converted.sum(),len(B)
print(f"\nsuccess-failure check: A {cA}/{nA-cA}, B {cB}/{nB-cB} (all >> 10 -> z is valid)")

Decision: proportion outcome + 2 independent groups -> TWO-PROPORTION z-TEST
H0: pB = pA   (the redesign does not change conversion)
H1: pB > pA   (the redesign converts better) -- one-sided

success-failure check: A 181/1649, B 237/1533 (all >> 10 -> z is valid)

# the test: counts and sample sizes for [B, A], one-sided "larger"
z, p_one = proportions_ztest([cB, cA], [nB, nA], alternative="larger")
# the 95% confidence interval for the lift (pB - pA)
lo, hi = confint_proportions_2indep(cB, nB, cA, nA, method="wald")
pA, pB = cA/nA, cB/nB
print(f"A = {pA:.2%}, B = {pB:.2%}, absolute lift = {(pB-pA)*100:+.2f} pts, relative lift = {(pB/pA-1)*100:+.1f}%")
print(f"two-proportion z-test: z = {z:.2f}, one-sided p = {p_one:.4f}")
print(f"95% CI for the lift: [{lo*100:+.2f}, {hi*100:+.2f}] pts")
print("REJECT H0: the redesign converts significantly better" if p_one<0.05 else "fail to reject H0")

A = 9.89%, B = 13.39%, absolute lift = +3.50 pts, relative lift = +35.4%
two-proportion z-test: z = 3.28, one-sided p = 0.0005
95% CI for the lift: [+1.40, +5.59] pts
REJECT H0: the redesign converts significantly better

# was the test adequately powered? Judge it at the design stage, from a 10% baseline and lifts set in advance
# (power at the observed lift only restates the p-value)
from statsmodels.stats.proportion import proportion_effectsize
from statsmodels.stats.power import NormalIndPower
pw, base = NormalIndPower(), 0.10
for lift in [0.02, 0.03]:
    power = pw.power(effect_size=proportion_effectsize(base+lift, base), nobs1=nA, alpha=0.05, ratio=nB/nA, alternative="larger")
    print(f"design power at ~{nA} per arm, {base:.0%} baseline, +{lift*100:.0f} pts: {power:.2f}")
h80 = pw.solve_power(nobs1=nA, alpha=0.05, power=0.80, ratio=nB/nA, alternative="larger")
mde = np.sin(np.arcsin(np.sqrt(base)) + h80/2)**2 - base   # turn Cohen's h back into a lift
print(f"lift detectable with 80% power: about +{mde*100:.1f} pts")
fig,ax=plt.subplots(figsize=(6.6,3.4))
ax.bar(["A (control)","B (redesign)"],[pA,pB],color=[LIGHT,CY],width=0.55)
ax.errorbar([0,1],[pA,pB],yerr=[1.96*np.sqrt(pA*(1-pA)/nA),1.96*np.sqrt(pB*(1-pB)/nB)],fmt="none",ecolor=INK,capsize=7,lw=1.6)
for i,v in enumerate([pA,pB]): ax.text(i,v+0.006,f"{v:.1%}",ha="center",fontweight="bold")
ax.set_ylim(0,0.17); ax.set_ylabel("conversion rate"); ax.set_title("Conversion rate by variant (95% CI)")
plt.tight_layout(); plt.show()

design power at ~1830 per arm, 10% baseline, +2 pts: 0.61
design power at ~1830 per arm, 10% baseline, +3 pts: 0.88
lift detectable with 80% power: about +2.6 pts

	count	unique	top	mean
session_id	3600	3600	S403599	NaN
variant	3600	2	A	NaN
device	3600	3	mobile	NaN
landing_page	3600	2	control	NaN
seconds_on_page	3600.0	NaN	NaN	54.232194
converted	3600.0	NaN	NaN	0.116111
visit_date	3600	10	2026-05-04	NaN

Case Study: A Web A/B Test 🌐

⚙️ Setup¶