3. Statistical Foundations

This chapter covers how to summarize data honestly, how to reason from a sample to the population it came from, and how to decide whether an observed difference is real or noise. It ends with causal inference, which addresses the question every stakeholder actually asks and that correlation alone cannot answer. Chapter 2 gave the probability machinery; this chapter is where it meets data you actually have.

<a id="31-descriptive-statistics"></a>

3.1 Descriptive Statistics

Descriptive statistics compress a dataset into a few numbers, and the whole skill lies in knowing what each number hides as well as what it reveals.

Measures of Central Tendency

Faced with a thousand values, the first question is always "what is the typical one?" Three answers compete, and picking the wrong one is the most common way a summary misleads without technically lying.

The mean is the arithmetic average: the balance point of the data. It uses every value, which is its strength and its weakness. A single extreme value drags it arbitrarily far.

The median is the middle value when sorted, with half the data on each side. It ignores magnitude entirely and responds only to rank, which makes it immune to outliers. Formally, its breakdown point is 50%: you must corrupt half the data before you can move it arbitrarily.

The mode is the most frequent value. It is the only one of the three that works on nominal categories, and it is the right answer to "what is the most common plan?"

The practical rule: for symmetric data the three roughly coincide and the mean is most efficient. For skewed data they diverge, and the median usually describes a typical case better. Since Chapter 2 established that revenue, latency, and session length are typically log-normal, this is not an edge case but the default situation in data science.

  SYMMETRIC                        RIGHT-SKEWED (typical of money, time)

        _--_                         |\
       /    \                        | \___
      /      \                       |     -----____
     ---------                      -+-----------------
         ^                           ^   ^        ^
     all three                    mode median   mean
     coincide                                    ^
                                     the mean is dragged toward the tail
                                     by a few very large values

  RULE OF THUMB FOR SKEW DIRECTION
   mean > median  ->  right-skewed (a long tail of large values)
   mean < median  ->  left-skewed
   mean = median  ->  symmetric (or symmetric-ish)

  BREAKDOWN POINT: how much corruption before the statistic is garbage
   mean:   1/n  (one bad value is enough)
   median: 50%  (you must corrupt half the data)

import numpy as np
from scipy import stats

salaries = np.array([42, 45, 47, 48, 50, 52, 55, 58, 61, 850])  # thousands
                                                    # one founder included

print(f"mean   = {salaries.mean():.1f}")            # 130.8
print(f"median = {np.median(salaries):.1f}")        # 51.0
print(f"mode   = {stats.mode(salaries, keepdims=False).mode}")   # 42

# The mean describes nobody here: nine of ten people earn far less.
print(f"share below the mean = {(salaries < salaries.mean()).mean():.0%}")  # 90%

# Remove the outlier and the mean collapses toward the median.
print(f"mean without founder = {salaries[:-1].mean():.1f}")     # 50.9

# A trimmed mean is the middle ground: drop a percentage from each end,
# keeping some efficiency of the mean with much of the median's robustness.
print(f"10% trimmed mean = {stats.trim_mean(salaries, 0.1):.1f}")   # 52.0

# For categorical data only the mode is defined.
plans = np.array(["free", "pro", "free", "premium", "free", "pro"])
print(stats.mode(plans, keepdims=False).mode)       # 'free'
# np.mean(plans) would be a TypeError, which is the correct outcome.

Measure	Uses all data	Outlier-robust	Works on categories	Best when
Mean	Yes	No	No	Symmetric, no outliers
Median	No (rank only)	Yes	Ordinal only	Skewed or contaminated
Mode	No	Yes	Yes	Categorical, multimodal

Key Takeaways

• The mean uses every value and so moves arbitrarily far with a single outlier; the median needs half the data corrupted.
• For skewed data the mean and median diverge, and the median usually describes a typical case better.
• Comparing mean to median is a fast skew diagnostic requiring no plot.
• Only the mode is defined for nominal categories.

🧪 Practice

1. Compute all three measures for a small salary dataset, then add one extreme value and recompute. Note which moved.
2. Generate log-normal data and report the percentage of observations falling below the mean.
3. Interview: A company reports average employee salary of 95,000 and the union reports 62,000. Can both be honest? (Hint: ask which statistic each side chose and what the distribution looks like.)

Measures of Dispersion

Two datasets can share a mean and be nothing alike. One might cluster tightly and the other sprawl. Dispersion measures the spread, and it is frequently more decision-relevant than the centre: a delivery service averaging 30 minutes with a 5-minute spread is a different business from one averaging 30 minutes with a 40-minute spread.

The range is maximum minus minimum. It is trivially computed and nearly useless, since it depends entirely on the two most extreme points.

Variance is the mean squared deviation from the mean, and standard deviation is its square root, returning to the original units. Chapter 2 defined these for a distribution; here the concern is estimating them from a sample. The subtlety is the divisor. Dividing by n underestimates the population variance, because deviations are measured from the sample mean, which is itself pulled toward the data. Dividing by n - 1 corrects this and is called Bessel's correction.

The interquartile range is the spread of the middle 50%, robust in the same way the median is. The median absolute deviation is more robust still.

The coefficient of variation, standard deviation divided by mean, is unitless and lets you compare variability across quantities on different scales.

  SAME MEAN, DIFFERENT SPREAD

   A:        ||||||||              sd = 1.2   "reliable"
   B:  |  |  |  |  |  |  |  |      sd = 4.8   "unpredictable"
       +---------+---------+
                mean

  WHY n-1 AND NOT n

   The sample mean sits, by construction, at the centre of YOUR sample,
   not at the true population mean. Deviations measured from it are
   therefore systematically too small.

      population mean  *
                        \
      sample:      o  o  o  o          deviations from the SAMPLE mean
                      ^                are smaller than from the TRUE mean
                 sample mean

   Dividing by n-1 inflates the estimate just enough to remove the bias.
   The effect vanishes as n grows: at n=1000 it is a 0.1% correction.

import numpy as np
from scipy import stats

a = np.array([48, 49, 50, 51, 52])
b = np.array([20, 35, 50, 65, 80])

for name, x in [("A", a), ("B", b)]:
    print(f"{name}: mean={x.mean():.1f} sd={x.std(ddof=1):.2f} "
          f"IQR={stats.iqr(x):.1f} range={x.ptp()}")
# A: mean=50.0 sd=1.58 IQR=2.0  range=4
# B: mean=50.0 sd=23.72 IQR=30.0 range=60

# ddof matters. NumPy defaults to ddof=0 (population); pandas defaults
# to ddof=1 (sample). The same column can give two different answers
# depending on which library computed it.
print(f"ddof=0: {b.std(ddof=0):.3f}   ddof=1: {b.std(ddof=1):.3f}")
# 21.213 vs 23.717 -- a 12% difference at n=5

# Demonstrate the bias that n-1 corrects.
rng = np.random.default_rng(0)
true_var = 9.0                                    # population variance
biased, unbiased = [], []
for _ in range(20_000):
    s = rng.normal(0, 3, size=5)                  # small samples exaggerate it
    biased.append(s.var(ddof=0))
    unbiased.append(s.var(ddof=1))
print(f"E[var, ddof=0] = {np.mean(biased):.3f}  (true {true_var})")   # ~7.2
print(f"E[var, ddof=1] = {np.mean(unbiased):.3f}  (true {true_var})") # ~9.0
# ddof=0 underestimates by a factor of (n-1)/n = 4/5 = 0.8, exactly.

# Robust alternatives, for data with outliers.
contaminated = np.append(a, 500)
print(f"sd  = {contaminated.std(ddof=1):.1f}")                  # 183.7
print(f"IQR = {stats.iqr(contaminated):.1f}")                   # 2.5
print(f"MAD = {stats.median_abs_deviation(contaminated):.1f}")  # 1.5
# One bad value destroyed the standard deviation and left IQR untouched.

# Coefficient of variation compares spread across different scales.
revenue = np.array([1000, 1200, 950, 1100])       # dollars
latency = np.array([45, 52, 41, 48])              # milliseconds
for name, x in [("revenue", revenue), ("latency", latency)]:
    print(f"{name}: CV = {x.std(ddof=1) / x.mean():.3f}")
# revenue: 0.104   latency: 0.100
# Nearly identical relative variability, even though one is measured in
# dollars around 1000 and the other in milliseconds around 46.

Key Takeaways

• Spread is often more decision-relevant than the centre, especially for service-level promises.
• Use ddof=1 for a sample variance; NumPy defaults to 0 and pandas to 1, which silently produces different numbers.
• Standard deviation is destroyed by a single outlier, while the interquartile range and median absolute deviation are not.
• The coefficient of variation is unitless and enables comparison across quantities measured on different scales.

🧪 Practice

1. Compute variance with both ddof=0 and ddof=1 on a sample of five and explain the gap.
2. Show empirically that the interquartile range is unaffected by adding an extreme value that doubles the standard deviation.
3. Interview: Two delivery services both average 30 minutes. What single additional number would you ask for, and why? (Hint: think about what a customer actually experiences on a bad day.)

Skewness and Kurtosis

Mean and standard deviation describe location and spread, but two distributions sharing both can still look completely different. The third and fourth standardized moments capture shape: whether the distribution leans, and how heavy its tails are.

Skewness measures asymmetry. Positive (right) skew means a long tail of large values, which is the signature of money, durations, and counts. Negative (left) skew means a long tail of small values, seen in exam scores near a ceiling or in time-to-deadline data. Zero skew means symmetric.

Kurtosis measures tail heaviness relative to a normal distribution. The commonly reported excess kurtosis subtracts 3 so that a normal scores zero. Positive excess kurtosis means heavier tails and more extreme outliers than normal; negative means lighter tails.

Kurtosis matters more than its obscurity suggests. Financial returns and latency distributions are leptokurtic, meaning extreme events happen far more often than a normal model predicts. Risk models that assume normality systematically underestimate the frequency of disasters, which is not an academic point.

The practical use of both is as a diagnostic before choosing a method. High skew suggests a transformation or a robust statistic. High kurtosis warns that outlier-sensitive methods will behave badly.

  SKEWNESS                           KURTOSIS (excess)

   negative      zero      positive    heavy tails (>0)   light tails (<0)
     __/|        _--_      |\__             /\                 ____
    /   |       /    \     |   \           /  \               /    \
   /    |      /      \    |    \_        /    \             /      \
  -------     ---------   --------      _/      \_          /        \
  tail left   symmetric   tail right    sharp peak,         flat, bounded
                                        fat tails           (e.g. uniform)

  ROUGH INTERPRETATION GUIDE

   |skew| < 0.5     approximately symmetric
   0.5 to 1.0       moderately skewed
   > 1.0            highly skewed -- consider a transform

   excess kurtosis > 0   more outliers than a normal would produce
   excess kurtosis < 0   fewer outliers, flatter than normal

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

samples = {
    "normal":      rng.normal(0, 1, 100_000),
    "log-normal":  rng.lognormal(0, 1, 100_000),
    "uniform":     rng.uniform(-2, 2, 100_000),
    "exponential": rng.exponential(1, 100_000),
    "t (df=3)":    rng.standard_t(3, 100_000),
}

for name, x in samples.items():
    # bias=False applies the small-sample correction; fisher=True (default)
    # returns EXCESS kurtosis, so a normal scores ~0 rather than ~3.
    print(f"{name:<12} skew={stats.skew(x, bias=False):>7.3f}  "
          f"excess kurtosis={stats.kurtosis(x, bias=False):>8.3f}")
# normal        skew~ 0.00   excess kurtosis~  0.00
# log-normal    skew~ 6.2    excess kurtosis~ 90+     highly skewed, fat tail
# uniform       skew~ 0.00   excess kurtosis~ -1.20   flat, light tails
# exponential   skew~ 2.00   excess kurtosis~  6.00   theoretical values
# t (df=3)      skew~ 0.00   excess kurtosis~ large   symmetric but fat-tailed

# Why kurtosis matters for risk: count the 4-sigma events.
for name in ("normal", "t (df=3)"):
    x = samples[name]
    z = (x - x.mean()) / x.std(ddof=1)
    print(f"{name}: |z| > 4 occurs {(np.abs(z) > 4).mean() * 100:.3f}% of the time")
# normal:   ~0.006%  -- matches the theoretical 0.0063%
# t (df=3): ~0.3%    -- roughly fifty times more often
# A model assuming normality would call each of those a once-in-a-decade
# event while they arrive monthly.

# Skew as a decision trigger: does a log transform fix it?
rev = samples["log-normal"]
print(f"skew before log = {stats.skew(rev, bias=False):.3f}")     # ~6.2
print(f"skew after log  = {stats.skew(np.log(rev), bias=False):.3f}")  # ~0.0

Key Takeaways

• Skewness measures asymmetry and kurtosis measures tail heaviness, both invisible in the mean and standard deviation.
• Positive excess kurtosis means extreme events are far more frequent than a normal distribution predicts.
• Absolute skew above 1 is a signal to consider a transformation or a robust statistic.
• Assuming normality for fat-tailed data systematically understates the frequency of extreme outcomes.

🧪 Practice

1. Compute skewness and excess kurtosis for a normal, a uniform, and an exponential sample and match them to the shapes.
2. Take a right-skewed column, apply a log transform, and report skewness before and after.
3. Interview: Why do financial risk models that assume normality fail during crises? (Hint: compare the predicted and observed frequency of five-sigma moves.)

Percentiles and Quantiles

A single centre-and-spread summary assumes the distribution has a simple shape. Quantiles make no such assumption: they describe the distribution by reporting where it sits at each point along its range, which is why they dominate in performance engineering and in any setting where the tail is the point.

A percentile is the value below which a given percentage of observations fall. The 90th percentile of latency is the time that 90% of requests beat. Quantiles are the same idea on a 0-to-1 scale, and quartiles are the special case splitting data into four parts.

The reason this matters: averages hide tails, and users experience tails. If mean latency is 100ms but the 99th percentile is 4 seconds, then one request in a hundred is unbearable. For a page making 50 backend calls, most page loads will contain at least one such call. Reporting the mean would describe a system nobody is using.

A wrinkle that trips people up: there is no single definition of a percentile. When the desired position falls between two data points, methods differ in how they interpolate. NumPy offers nine. The differences are negligible on large samples and can be material on small ones.

  THE BOX PLOT IS A QUANTILE SUMMARY

        Q1     median    Q3
         |       |        |
    |----[=======|========]----|        o     o
    |    |                |    |
   min  25%      50%      75%  max      outliers beyond the whiskers
        |<------ IQR ---->|
        |<-1.5*IQR->|          |<-1.5*IQR->|
        whisker limits: the furthest data point still inside the fence

  WHY THE TAIL IS THE PRODUCT

   latency distribution
    |####
    |#######
    |##########
    |#############_
    |                ----____________________
    +-----------------------------------------> ms
     50   100                            4000
      ^                                    ^
    mean = 100ms                     p99 = 4000ms
    "the system is fast"          "1 in 100 users waits 4 seconds"

   With 50 backend calls per page, P(at least one slow call)
     = 1 - 0.99^50 = 39.5% of page loads.

import numpy as np

rng = np.random.default_rng(0)
# Realistic latency: mostly fast, with a slow tail.
latency = np.concatenate([
    rng.lognormal(mean=3.5, sigma=0.4, size=9800),    # the healthy majority
    rng.lognormal(mean=7.0, sigma=0.8, size=200),     # the slow 2%
])

for q in (50, 90, 95, 99, 99.9):
    print(f"p{q:<5} = {np.percentile(latency, q):8.1f} ms")
print(f"mean  = {latency.mean():8.1f} ms")
# p50   ~    33 ms
# p90   ~    58 ms
# p95   ~    70 ms
# p99   ~  1100 ms     <- more than thirty times the median
# p99.9 ~  4000 ms
# mean  ~    65 ms     <- describes neither the typical nor the bad case

# The compounding problem, computed rather than asserted.
p_slow = 0.01
for calls in (1, 10, 50):
    print(f"{calls:>2} calls -> P(at least one p99 hit) = "
          f"{1 - (1 - p_slow) ** calls:.1%}")
#  1 calls ->  1.0%
# 10 calls ->  9.6%
# 50 calls -> 39.5%    tail latency is a page-level problem, not a rare one

# Interpolation method matters on small samples.
small = np.array([1, 2, 3, 4])
for method in ("linear", "lower", "higher", "midpoint", "nearest"):
    print(f"{method:<9} p75 = {np.percentile(small, 75, method=method)}")
# linear 3.25 | lower 3 | higher 4 | midpoint 3.5 | nearest 3
# All are defensible. Fix one method and use it consistently, especially
# when comparing numbers computed by different tools.

# A critical operational point: percentiles do NOT average.
shard_p99 = [1400, 1350, 1500]
print(f"mean of shard p99s = {np.mean(shard_p99):.0f} ms")   # 1417 -- WRONG
# The true global p99 requires the pooled distribution, or a mergeable
# sketch such as t-digest or HDR histogram. Averaging percentiles across
# shards, hosts, or time buckets is one of the most common metric bugs.

Key Takeaways

• Quantiles describe a distribution without assuming its shape, which is why tail-sensitive systems are specified with them.
• Users experience tail latency, and averages conceal it entirely.
• Tail events compound: a 1% slow call rate makes a 50-call page slow 40% of the time.
• Percentiles cannot be averaged across shards or time windows; pool the data or use a mergeable sketch.

🧪 Practice

1. Generate a skewed latency sample and report p50, p95, p99 alongside the mean. Describe what each tells an on-call engineer.
2. Compute p75 of a four-element array under three interpolation methods and explain the differences.
3. Interview: A dashboard averages hourly p99 latency into a daily number. What is wrong with it? (Hint: ask what the average of two percentiles is a percentile of.)

Correlation Coefficients

Chapter 2 introduced correlation as rescaled covariance. Here the concern is practical: which coefficient to use, and what each one can and cannot detect. Choosing the wrong one produces confident conclusions about relationships that are not there, or misses relationships that are.

Pearson's r measures linear association. It assumes both variables are roughly interval-scaled and is sensitive to outliers, since it is built from means and standard deviations. A single extreme point can create or destroy it.

Spearman's rho is Pearson's correlation computed on ranks instead of values. Because ranks discard magnitude, it detects any monotonic relationship, linear or not, and is robust to outliers. It is the right default for ordinal data.

Kendall's tau counts concordant and discordant pairs. It is more interpretable as a probability statement, more robust with small samples, and slower to compute.

For categorical associations, Cramer's V generalizes the chi-square statistic to a 0-to-1 measure, and the point-biserial correlation handles one binary against one continuous variable.

The universal caution remains: correlation is not causation, and a coefficient collapses a scatter plot into one number that cannot distinguish a clean linear trend from a plot with entirely different structure.

  ANSCOMBE'S QUARTET: four datasets, identical statistics

   I            II           III          IV
    .  '         . . .          .           |
   . '         .       .       .            |
  . '         .         .     .             |
 . '          .         .    .              .
 -------      ---------      ------      ---------
 clean        curved        outlier      one point
 linear       (parabola)    drives it    drives everything

 ALL FOUR: mean x = 9, mean y = 7.5, Pearson r = 0.816, same
 regression line. Only one of them should be modelled with a line.

  WHAT EACH COEFFICIENT SEES

   relationship        Pearson    Spearman    Kendall
   perfectly linear      1.00       1.00       1.00
   monotonic curve       ~0.9       1.00       1.00     <- rank-based win
   U-shaped              ~0.0       ~0.0       ~0.0     <- all blind; plot it
   linear + 1 outlier    ~0.0       ~0.9       ~0.9     <- robustness wins

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
x = np.arange(1, 21, dtype=float)

cases = {
    "linear":          2 * x + rng.normal(0, 1, 20),
    "monotonic curve": np.exp(x / 4),            # increasing but very curved
    "u-shaped":        (x - 10.5) ** 2,
}

for name, y in cases.items():
    r = stats.pearsonr(x, y)
    rho = stats.spearmanr(x, y)
    tau = stats.kendalltau(x, y)
    print(f"{name:<16} pearson={r.statistic:>6.3f}  "
          f"spearman={rho.statistic:>6.3f}  kendall={tau.statistic:>6.3f}")
# linear           pearson~ 0.998  spearman~ 0.998  kendall~ 0.98
# monotonic curve  pearson~ 0.83   spearman= 1.000  kendall= 1.000
# u-shaped         pearson= 0.000  spearman= 0.000  kendall= 0.000
#                  ^ the curve is PERFECTLY monotonic; only rank methods say so
#                  ^ the U is a perfect deterministic relation; none detect it

# Outlier sensitivity, demonstrated.
a = np.arange(20, dtype=float)
b = a + rng.normal(0, 0.5, 20)
print(f"clean:   pearson={stats.pearsonr(a, b).statistic:.3f}  "
      f"spearman={stats.spearmanr(a, b).statistic:.3f}")
a2, b2 = np.append(a, 100.0), np.append(b, 0.0)      # one wild point
print(f"+outlier: pearson={stats.pearsonr(a2, b2).statistic:.3f}  "
      f"spearman={stats.spearmanr(a2, b2).statistic:.3f}")
# clean:    pearson=0.996  spearman=0.995
# +outlier: pearson~0.0    spearman~0.86
# One point destroyed Pearson and barely dented Spearman.

# Categorical association: Cramer's V from a contingency table.
def cramers_v(table):
    chi2 = stats.chi2_contingency(table, correction=False).statistic
    n = table.sum()
    return np.sqrt(chi2 / (n * (min(table.shape) - 1)))

table = np.array([[30, 10],
                  [ 5, 55]])
print(f"Cramer's V = {cramers_v(table):.3f}")        # strong association

Coefficient	Detects	Robust to outliers	Data type
Pearson r	Linear	No	Two continuous
Spearman rho	Monotonic	Yes	Continuous or ordinal
Kendall tau	Monotonic	Yes	Continuous or ordinal
Point-biserial	Linear	No	One binary, one continuous
Cramer's V	Any association	N/A	Two categorical

Key Takeaways

• Pearson detects linear association only and is destroyed by a single outlier; Spearman and Kendall work on ranks and are robust.
• A perfectly monotonic curve gives Spearman 1.0 and Pearson noticeably less, which is a useful diagnostic pair.
• No correlation coefficient detects a symmetric non-monotonic relationship such as a parabola.
• Anscombe's quartet shows that identical coefficients can describe completely different data, so always plot.

🧪 Practice

1. Compute Pearson and Spearman on an exponential relationship and explain why they differ.
2. Reproduce Anscombe's quartet and verify the four datasets share summary statistics.
3. Interview: Ice cream sales correlate with drownings at r = 0.9. What do you conclude? (Hint: look for a variable that drives both.)

Contingency Tables

When both variables are categorical, scatter plots and correlation coefficients do not apply. A contingency table, also called a cross-tab, is the right summary: a grid of counts for every combination of categories. It is the categorical analogue of the joint distribution from Chapter 2, filled with observed counts rather than probabilities.

The table itself is just bookkeeping. What makes it powerful is the comparison between what you observed and what you would expect if the two variables were independent. Under independence, the expected count in a cell is the row total times the column total divided by the grand total. Large gaps between observed and expected are evidence of association, and the chi-square test formalizes that comparison.

Reading a cross-tab correctly requires choosing the right percentages. Row percentages, column percentages, and percentages of the total answer different questions, and confusing them reverses conclusions. This is the same trap as confusing P(A | B) with P(B | A).

The final caution is Simpson's paradox from Chapter 2: a two-way table can mislead whenever a third variable differs across the groups. Always ask what is being collapsed.

  A CONTINGENCY TABLE AND ITS THREE PERCENTAGE VIEWS

   OBSERVED COUNTS
                churned  retained | total
      free          120      380  |   500
      pro            40      460  |   500
      ----------------------------+------
      total         160      840  |  1000

   EXPECTED IF INDEPENDENT  (row total * col total / grand total)
      free: 500*160/1000 = 80     500*840/1000 = 420
      pro:  500*160/1000 = 80     500*840/1000 = 420

   Observed 120 vs expected 80 for free churn -> far more than chance.

   ROW %      "of free users, what share churned?"     120/500 = 24%
   COLUMN %   "of churners, what share were free?"     120/160 = 75%
   TOTAL %    "what share of all users are free and churned?"  12%

   These are three different numbers answering three different
   questions. Quoting the wrong one is the base rate fallacy again.

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

df = pd.DataFrame({
    "plan":   ["free"] * 500 + ["pro"] * 500,
    "churned": [1] * 120 + [0] * 380 + [1] * 40 + [0] * 460,
})

table = pd.crosstab(df["plan"], df["churned"])
print(table)
#  churned    0    1
#  plan
#  free     380  120
#  pro      460   40

# The three percentage views, each answering a different question.
print(pd.crosstab(df["plan"], df["churned"], normalize="index").round(3))
#   of each PLAN, what fraction churned?     free 0.24, pro 0.08
print(pd.crosstab(df["plan"], df["churned"], normalize="columns").round(3))
#   of each OUTCOME, what fraction was each plan?   churners: 75% free
print(pd.crosstab(df["plan"], df["churned"], normalize="all").round(3))
#   share of the whole dataset in each cell

# Is the association real, or could it be chance?
result = stats.chi2_contingency(table)
print(f"chi2 = {result.statistic:.1f}, p = {result.pvalue:.2e}")
print("expected counts under independence:\n", result.expected_freq.round(1))
# expected: [[420, 80], [420, 80]] against observed [[380,120],[460,40]]
# chi2 is large and p is tiny: the plans genuinely differ in churn.

# Effect size, because with enough rows everything is "significant".
n = table.values.sum()
cramers_v = np.sqrt(result.statistic / (n * (min(table.shape) - 1)))
print(f"Cramer's V = {cramers_v:.3f}")        # ~0.22, a moderate association

# Always add margins when presenting: readers need the denominators.
print(pd.crosstab(df["plan"], df["churned"], margins=True))

Key Takeaways

• A contingency table cross-tabulates counts and is the right summary for two categorical variables.
• Association is judged by comparing observed counts to those expected under independence.
• Row, column, and total percentages answer different questions, and quoting the wrong one inverts the conclusion.
• Report an effect size such as Cramer's V alongside a chi-square p-value, since large samples make trivial associations significant.

🧪 Practice

1. Build a cross-tab of two categorical columns and produce all three percentage views, stating the question each answers.
2. Compute the expected counts under independence by hand for a 2-by-2 table and compare them to the observed values.
3. Interview: A cross-tab shows a treatment works overall, but you learn patients chose their own treatment. What now? (Hint: ask what else differs between the two groups.)

<a id="32-sampling-and-estimation"></a>

3.2 Sampling and Estimation

Every dataset is a sample, so every conclusion is an inference about a population you cannot see; this section covers how samples are drawn, how they mislead, and how to attach honest uncertainty to what you estimate.

Population vs Sample

The distinction between a population and a sample sounds pedantic until it silently invalidates a conclusion, which it frequently does.

The population is the complete set of units you want to describe: all customers, all transactions, all possible future requests. It is usually unobservable, either because it is too large, because it includes the future, or because measuring all of it is impossible. Its true characteristics are parameters, conventionally written with Greek letters.

The sample is the subset you actually observe. Quantities computed from it are statistics, and they are estimates of the parameters. The entire field of inferential statistics exists to quantify the gap between the two.

The crucial and frequently missed point is that the population must be defined by the question, not by the data you happen to have. If you want to know how all users will react and you only have data from beta opt-ins, your population is beta opt-ins, no matter how large the sample. More data does not fix the mismatch. It just makes the wrong answer more precise.

An analogy: tasting soup. A spoonful from a well-stirred pot tells you about the pot. A spoonful skimmed off the top tells you about the top, and drinking a litre of it still tells you only about the top.

  POPULATION                        SAMPLE

  +-----------------------------+
  |  all customers (unknown)    |      draw
  |   parameters:               |  ----------->   n = 500 observed
  |     mu     (true mean)      |                 statistics:
  |     sigma  (true sd)        |                   x-bar  (sample mean)
  |     p      (true rate)      |  <-----------     s      (sample sd)
  |                             |     infer         p-hat  (sample rate)
  +-----------------------------+

  NOTATION CONVENTION
   parameter (population)   statistic (sample)   estimates
   mu                       x-bar                the mean
   sigma                    s                    the sd
   p                        p-hat                a proportion

  THE FATAL MISMATCH

   question's population:  "all our users"
   sample's population:    "users who opted into the beta"
                            \_____________________________/
                             no sample size fixes this gap

import numpy as np

rng = np.random.default_rng(0)

# A population we can see only because we invented it.
population = rng.normal(loc=100, scale=15, size=1_000_000)
print(f"PARAMETERS: mu={population.mean():.2f} sigma={population.std():.2f}")

# Each sample gives a different estimate. That variation IS the subject.
for n in (10, 100, 1000):
    estimates = [rng.choice(population, n, replace=False).mean()
                 for _ in range(1000)]
    print(f"n={n:>5}: estimates range {min(estimates):.1f} to "
          f"{max(estimates):.1f}, sd of estimates = {np.std(estimates):.2f}")
# n=   10: wide spread, sd of estimates ~4.7
# n=  100: narrower,     sd ~1.5
# n= 1000: narrow,       sd ~0.47
# The sd of the estimates is the STANDARD ERROR, and it shrinks as
# sigma/sqrt(n) exactly as the central limit theorem predicts.

# The mismatch that sample size cannot fix.
biased_frame = population[population > 110]        # only high-value users
big_biased_sample = rng.choice(biased_frame, 100_000, replace=False)
print(f"huge but biased sample mean = {big_biased_sample.mean():.2f}")
# ~118, versus a true mu of 100. The sample is enormous and the answer
# is still wrong, because the sampling FRAME excluded most of the
# population. Precision and accuracy are different things.

Key Takeaways

• Parameters describe the unobservable population; statistics are computed from the sample and estimate them.
• The population is defined by the question being asked, not by the data that happens to be available.
• Sample-to-sample variation in an estimate is the standard error, and it shrinks with the square root of n.
• A biased sampling frame is not fixed by more data; it produces a precise wrong answer.

🧪 Practice

1. Simulate a population, draw 1000 samples at three sizes, and plot the spread of the sample means.
2. For a project you know, write one sentence naming the population the question implies and one naming the population the data represents.
3. Interview: Your survey has 50,000 responses and a 2% response rate. Is it reliable? (Hint: ask who the 98% are and whether they differ.)

Sampling Methods

Given a population, how you choose the sample determines what the sample can tell you. The methods differ in cost, in the precision they deliver per unit of effort, and in the assumptions they require.

Simple random sampling gives every unit an equal chance of selection. It is the reference standard because it is unbiased and the mathematics is straightforward. Its weakness is practical: it requires a complete list of the population and can miss small subgroups entirely by chance.

Stratified sampling divides the population into groups, called strata, and samples within each. When the strata differ meaningfully on the outcome, this produces more precise estimates than simple random sampling at the same cost, because it removes between-strata variation from the error.

Cluster sampling divides the population into naturally occurring groups and samples entire clusters. It is cheaper when travel or access is costly, and less precise, because units within a cluster resemble each other and so carry less independent information.

Systematic sampling takes every k-th unit. It is simple and usually fine, but it fails badly when the list has a periodicity matching k.

Convenience sampling takes whatever is easy to reach. It is the most common method in practice and supports no valid inference to a wider population.

  SIMPLE RANDOM            STRATIFIED              CLUSTER

  o o . o . . o o          [AAA][BBB][CCC]        [ooo][...][ooo]
  . o . . o o . .           o.o   .o.   oo.        ooo   ...   ooo
  o . o o . . o .           sample within          sample WHOLE
  . . o . o o . o           every stratum          clusters
                                                   
  each unit equally        guarantees coverage     cheap access,
  likely; may miss         of each group;          units within a
  small subgroups          MORE precise            cluster are alike,
                           when strata differ      so LESS precise

  WHY STRATIFICATION HELPS

   population: 90% free users (churn 20%), 10% pro users (churn 5%)

   simple random n=100:  expect 10 pro users, but could get 4 or 17,
                         and that swing moves the overall estimate

   stratified n=100:     force exactly 90 free and 10 pro, removing
                         that source of variation entirely

import numpy as np
import pandas as pd

rng = np.random.default_rng(0)
N = 100_000
# A population where a small stratum behaves very differently.
plan = np.where(rng.random(N) < 0.10, "pro", "free")
spend = np.where(plan == "pro",
                 rng.normal(500, 60, N),     # small group, high spend
                 rng.normal(50, 15, N))      # large group, low spend
pop = pd.DataFrame({"plan": plan, "spend": spend})
true_mean = pop["spend"].mean()
print(f"true mean spend = {true_mean:.2f}")     # ~95

def simple_random(df, n):
    return df.sample(n, random_state=None)["spend"].mean()

def stratified(df, n):
    # Allocate the sample across strata in proportion to their size.
    parts = [g.sample(max(1, int(round(n * len(g) / len(df)))))
             for _, g in df.groupby("plan")]
    return pd.concat(parts)["spend"].mean()

for method, fn in [("simple random", simple_random), ("stratified", stratified)]:
    estimates = [fn(pop, 200) for _ in range(500)]
    print(f"{method:<14} mean={np.mean(estimates):7.2f}  "
          f"sd of estimates={np.std(estimates):6.2f}")
# simple random  mean~95   sd of estimates ~9.7
# stratified     mean~95   sd of estimates ~1.7
# Both are unbiased. Stratification cuts the error roughly five-fold at
# the same sample size, purely by removing the randomness in how many
# expensive "pro" users happen to be drawn. With 10% pro users and
# n=200, a simple random sample might contain 12 of them or 28, and
# each one carries ten times the spend of a free user.

Method	Requires a frame	Cost	Precision	Main risk
Simple random	Yes, complete	Medium	Baseline	Misses small subgroups
Stratified	Yes, plus strata	Medium	Best	Needs a good stratifier
Cluster	Cluster list only	Low	Worst	Within-cluster similarity
Systematic	Ordered list	Low	Baseline	Periodicity in the list
Convenience	None	Lowest	Invalid	Unquantifiable bias

Key Takeaways

• Stratified sampling beats simple random sampling when strata genuinely differ on the outcome, at the same cost.
• Cluster sampling trades precision for access cost, because units inside a cluster carry redundant information.
• Systematic sampling is fine unless the list has a period matching the sampling interval.
• Convenience samples support description of the sample only, never inference to a wider population.

🧪 Practice

1. Simulate a population with two very different strata and compare the error of simple random against stratified sampling.
2. Construct a list with a weekly period and show how systematic sampling every 7th record produces a biased estimate.
3. Interview: You must estimate average household income in a city and visiting each home is expensive. Which design do you choose? (Hint: weigh travel cost against the similarity of neighbours.)

Sampling Bias

Sampling bias is systematic error in who ends up in your data, and it is the single most dangerous problem in applied statistics, because it is invisible in the data itself. Every diagnostic you can compute uses the biased sample, so nothing in the numbers warns you.

Selection bias occurs when the sampling mechanism relates to the outcome. Survivorship bias is the special case where units that failed have disappeared from the data before you look. Non-response bias arises when those who decline to answer differ from those who respond. Undercoverage means part of the population had no chance of selection at all.

The classic illustration is wartime aircraft armour. Analysts examined returning planes, mapped the bullet holes, and proposed armouring the places with the most holes. Abraham Wald pointed out the inversion: those were the survivable hits. The planes hit elsewhere never returned, so the data contained no examples of the fatal locations. The correct action was to armour where the returning planes showed no damage.

This structure recurs constantly in data science. A credit model trained only on approved applicants never observes how rejected applicants would have performed. A churn model trained on current customers is missing everyone who already left.

  SURVIVORSHIP BIAS: THE MISSING DATA IS THE ANSWER

   what you see:  planes that RETURNED
    +---------------+
    |  x      x  x  |      x = bullet holes, concentrated on wings
    | x  [ ]     x  |      [ ] = engines, suspiciously clean
    |  x    x   x   |
    +---------------+

   what you do NOT see: planes hit in the engines. They crashed.

   naive conclusion: armour the wings (most holes)
   correct conclusion: armour the engines (no holes among SURVIVORS)

  THE SAME SHAPE IN DATA SCIENCE

   observed          invisible                   consequence
   ---------------   -------------------------   --------------------
   approved loans    how rejects would perform   model cannot learn
                                                 the reject region
   current users     users who already churned   survivor-only model
   completed orders  abandoned carts             conversion overstated
   responded surveys the silent majority         satisfaction inflated

import numpy as np

rng = np.random.default_rng(0)
n = 20_000

# A lending population: credit score drives default risk.
score = rng.normal(650, 80, n)
default_prob = 1 / (1 + np.exp((score - 600) / 40))     # lower score, more risk
defaulted = rng.random(n) < default_prob

# The bank only ever observes applicants it approved.
approved = score > 620
print(f"approval rate = {approved.mean():.1%}")             # ~65%

print(f"true default rate, everyone   = {defaulted.mean():.3f}")
print(f"observed rate, approved only  = {defaulted[approved].mean():.3f}")
# true ~0.29 vs observed ~0.13 -- the observed data understates risk
# by more than half, and nothing inside the approved data reveals this.

# The subtler damage: the relationship itself looks weaker.
lo, hi = score < 620, score >= 620
print(f"score range in full population: {score.min():.0f} to {score.max():.0f}")
print(f"score range in approved data:   {score[approved].min():.0f} to "
      f"{score[approved].max():.0f}")
# Truncating the predictor's range flattens the estimated relationship,
# so the model learns that score matters LESS than it really does.

# Detection requires data from outside the sample. There is no internal
# test. Options in practice:
#   - approve a small random sample regardless of score (a "reject
#     inference" holdout), accepting the losses as the cost of learning
#   - compare your sample's marginals to a trusted external census
#   - reweight by known population proportions if coverage is merely uneven
external_population_mean_score = 650
print(f"sample mean score {score[approved].mean():.0f} vs external "
      f"{external_population_mean_score} -> coverage gap detected")

Key Takeaways

• Sampling bias is systematic and invisible from inside the data, so no internal diagnostic can detect it.
• Survivorship bias means the most informative cases are the ones missing from the dataset.
• Truncating a predictor's range flattens the apparent relationship, making important variables look weak.
• Detection requires external information: a census, a random holdout, or deliberate exploration outside the current policy.

🧪 Practice

1. Simulate a selection rule that keeps only high values of a predictor and show how the estimated slope changes.
2. List three sources of survivorship bias in a dataset you work with.
3. Interview: A model predicts loan default and is trained only on approved loans. What is the flaw and how would you mitigate it? (Hint: think about what outcome data exists for the applicants you turned down.)

Point Estimators and Their Properties

An estimator is a recipe for computing a guess at a parameter from a sample. The sample mean is an estimator of the population mean. Since many different recipes could estimate the same parameter, you need criteria for comparing them.

Bias is the difference between the estimator's average value across repeated samples and the true parameter. An unbiased estimator is centred on the truth. Variance is how much the estimate jumps around from sample to sample. Consistency means the estimator converges to the truth as the sample grows. Efficiency compares variance among unbiased estimators: the efficient one uses the data best.

The critical insight is that bias and variance trade off, and unbiasedness is not automatically preferable. An estimator that is slightly biased but much less variable will be closer to the truth on average, as measured by mean squared error, which decomposes exactly into bias squared plus variance. This is the same trade-off that reappears in Chapter 9 as the bias-variance decomposition for models, and it is why regularization works.

An analogy: a scale that reads 200 grams heavy every time is biased but precise, and you can correct it. A scale that averages correctly but swings by 5 kilograms is unbiased and useless.

  BIAS AND VARIANCE AS DARTBOARD ACCURACY

   low bias, low var      low bias, high var
      +-------+              +-------+
      |  .o.  |              | o   o |
      |  (x)  |              |  (x)o |
      |  .o.  |              |o    o |
      +-------+              +-------+
      ideal                  correct on average, unreliable

   high bias, low var     high bias, high var
      +-------+              +-------+
      |   (x) |              | o (x) |
      |    ooo|              |o     o|
      |    ooo|              |   o  o|
      +-------+              +-------+
      consistently off       hopeless
      but correctable

  MSE = BIAS^2 + VARIANCE

   A biased estimator can beat an unbiased one whenever the variance
   it saves exceeds the bias squared it adds. That is the entire
   justification for shrinkage and regularization.

import numpy as np

rng = np.random.default_rng(0)
TRUE_MU = 100.0
TRUE_SIGMA = 15.0

# Three estimators of the population mean.
estimators = {
    "sample mean":   lambda s: s.mean(),
    "median":        lambda s: np.median(s),
    "shrunk to 90":  lambda s: 0.8 * s.mean() + 0.2 * 90.0,   # biased on purpose
}

results = {name: [] for name in estimators}
for _ in range(20_000):
    sample = rng.normal(TRUE_MU, TRUE_SIGMA, size=10)
    for name, fn in estimators.items():
        results[name].append(fn(sample))

print(f"{'estimator':<14} {'bias':>8} {'variance':>10} {'MSE':>10}")
for name, vals in results.items():
    vals = np.array(vals)
    bias = vals.mean() - TRUE_MU
    var = vals.var()
    mse = np.mean((vals - TRUE_MU) ** 2)
    print(f"{name:<14} {bias:>8.3f} {var:>10.3f} {mse:>10.3f}")
# sample mean     ~0.00    ~22.5      ~22.5    unbiased, efficient
# median          ~0.00    ~33.7      ~33.7    unbiased, LESS efficient
# shrunk to 90    ~-2.00   ~14.4      ~18.4    BIASED but lowest MSE
#
# The shrunk estimator is wrong on average by 2 and still beats the
# unbiased sample mean on total error, because it is far less variable.
# Verify the decomposition: MSE = bias^2 + variance = 4.0 + 14.4 = 18.4.

# Consistency: bias and variance both vanish as n grows.
for n in (10, 100, 1000):
    est = [rng.normal(TRUE_MU, TRUE_SIGMA, n).mean() for _ in range(5000)]
    print(f"n={n:>5}: mean={np.mean(est):.3f}  sd={np.std(est):.3f}")
# The sd falls by roughly sqrt(10) each time n grows tenfold.

Property	Definition	Why it matters
Unbiased	Average over samples equals the parameter	No systematic error
Consistent	Converges to the parameter as n grows	More data actually helps
Efficient	Lowest variance among unbiased estimators	Best use of a fixed sample
Robust	Insensitive to outliers and assumption breaks	Survives messy real data

Key Takeaways

• Estimators are compared on bias, variance, consistency, and efficiency, not on any single one of these.
• Mean squared error decomposes exactly into bias squared plus variance.
• A biased estimator with low variance can beat an unbiased one, which is the justification for shrinkage and regularization.
• The sample median is unbiased for a symmetric distribution but less efficient than the mean.

🧪 Practice

1. Compare the sample mean and median as estimators of a normal mean by simulating bias, variance, and mean squared error.
2. Show that dividing by n rather than n - 1 gives a biased variance estimator and quantify the bias.
3. Interview: When would you prefer a biased estimator? (Hint: write out the mean squared error decomposition and ask which term dominates at small sample sizes.)

Maximum Likelihood Estimation

Most estimators so far were intuitive: to estimate a mean, average the data. But for a complicated model there is no obvious recipe. Maximum likelihood is the general-purpose principle that supplies one, and it underlies logistic regression, most distribution fitting, and the training objective of nearly every probabilistic model.

The idea inverts the usual question. Rather than asking "given the parameters, how probable is this data?", it asks "which parameter values would make the observed data most probable?" The likelihood is the probability of the data viewed as a function of the parameters, and the maximum likelihood estimate is the parameter value that maximizes it.

In practice you maximize the log likelihood instead. The log turns products into sums, which is numerically essential because multiplying thousands of probabilities underflows to zero in floating point. Since the logarithm is monotonic, the maximizing parameter is unchanged. Optimizers minimize by convention, so the quantity actually minimized is the negative log likelihood, which is where the name of that loss function comes from.

An analogy: a detective with a fingerprint asks which suspect makes the evidence least surprising. Maximum likelihood is that reasoning made arithmetic.

  THE LIKELIHOOD FUNCTION

   data: 7 heads in 10 flips. Which p makes that most probable?

   L(p)
    |            _--_
    |          .'    '.              peak at p = 0.7,
    |        .'        '.            which is exactly 7/10
    |      .'            '.
    |   .-'                '-.
    +--+----+----+----+----+---> p
      0   0.25  0.5  0.7  1.0

  WHY LOG: PRODUCTS UNDERFLOW

   L(p) = prod of 1000 probabilities, each about 0.5
        = 0.5^1000 = 9.3e-302 ... and at 1100 terms it is exactly 0.0

   log L = sum of 1000 log-probabilities = -693.1
        perfectly representable, and the argmax is identical

  RELATION TO FAMILIAR LOSSES
   normal errors + MLE       ->  minimize squared error
   Bernoulli labels + MLE    ->  minimize cross-entropy
   Both "loss functions" ARE negative log likelihoods in disguise.

import numpy as np
from scipy import stats, optimize

# Analytic case: MLE for a coin's bias is just the observed proportion.
flips = np.array([1, 1, 0, 1, 1, 0, 1, 1, 1, 0])      # 7 heads of 10

def neg_log_likelihood_bernoulli(p, data):
    p = np.clip(p, 1e-12, 1 - 1e-12)                  # guard log(0)
    return -np.sum(data * np.log(p) + (1 - data) * np.log(1 - p))

grid = np.linspace(0.01, 0.99, 99)
nll = [neg_log_likelihood_bernoulli(p, flips) for p in grid]
print(f"MLE by grid search = {grid[int(np.argmin(nll))]:.2f}")   # 0.70
print(f"closed form        = {flips.mean():.2f}")                # 0.70

# Numerical case: fit a gamma distribution where no simple formula exists.
rng = np.random.default_rng(0)
data = rng.gamma(shape=2.5, scale=3.0, size=2000)

def neg_log_likelihood_gamma(params, x):
    shape, scale = params
    if shape <= 0 or scale <= 0:
        return np.inf                                  # reject invalid values
    return -np.sum(stats.gamma.logpdf(x, a=shape, scale=scale))

fit = optimize.minimize(neg_log_likelihood_gamma, x0=[1.0, 1.0],
                        args=(data,), method="Nelder-Mead")
print(f"MLE shape={fit.x[0]:.3f} scale={fit.x[1]:.3f}  (true 2.5, 3.0)")

# SciPy does the same thing through a fitting interface.
print(stats.gamma.fit(data, floc=0))      # floc=0 fixes the location

# The connection that makes MLE central: least squares IS maximum
# likelihood when the errors are assumed normal with constant variance.
x = rng.normal(size=200)
y = 2 * x + 1 + rng.normal(scale=0.5, size=200)

def nll_linear(params, x, y):
    a, b, sigma = params
    if sigma <= 0:
        return np.inf
    return -np.sum(stats.norm.logpdf(y, loc=a * x + b, scale=sigma))

mle = optimize.minimize(nll_linear, [0, 0, 1], args=(x, y)).x
ols = np.polyfit(x, y, 1)
print(f"MLE slope={mle[0]:.4f}  OLS slope={ols[0]:.4f}")   # identical

Key Takeaways

• Maximum likelihood chooses the parameters that make the observed data most probable.
• Always work with the log likelihood, because products of many probabilities underflow to zero.
• Minimizing negative log likelihood is the same operation, phrased for optimizers that minimize.
• Squared error loss and cross-entropy loss are negative log likelihoods under normal and Bernoulli assumptions respectively.

🧪 Practice

1. Derive the maximum likelihood estimate for a Poisson rate and confirm it equals the sample mean.
2. Fit a gamma distribution numerically and compare your result to scipy.stats.gamma.fit.
3. Interview: Why is log likelihood used instead of likelihood? (Hint: compute the product of a thousand values near 0.5 in floating point.)

Method of Moments

Maximum likelihood requires a likelihood function and usually a numerical optimizer. The method of moments is an older, simpler alternative that often gives a usable answer in closed form, and it remains valuable as a starting point for optimization and as a sanity check.

The principle is direct. A distribution's parameters determine its theoretical moments: the mean, the variance, and higher quantities. The method sets the theoretical moments equal to the moments computed from the sample, then solves the resulting equations for the parameters. If a distribution has two parameters, you match two moments.

Compared with maximum likelihood, it is usually easier to compute and often available analytically, but it is generally less efficient, meaning it has higher variance. It can also produce impossible estimates, such as a negative variance parameter, because nothing in the procedure enforces the parameter space.

The standard practical use is to supply starting values for a maximum likelihood optimizer, which converges faster and more reliably from a sensible starting point than from an arbitrary one.

  MATCHING MOMENTS

   theoretical (from the parameters)     sample (from the data)
   ------------------------------------  ----------------------
   gamma:  mean     = shape * scale  ==  x-bar
           variance = shape * scale^2 == s^2

   Two equations, two unknowns. Solve:

       scale = s^2 / x-bar
       shape = x-bar^2 / s^2

   No optimizer, no iteration, no likelihood function required.

  METHOD OF MOMENTS vs MAXIMUM LIKELIHOOD

              MoM                        MLE
   form       often closed form          usually numerical
   effort     low                        higher
   efficiency lower (more variance)      highest achievable
   validity   can give impossible values parameters stay in range
   typical use starting values, checks   the final estimate

import numpy as np
from scipy import stats, optimize

rng = np.random.default_rng(0)
TRUE_SHAPE, TRUE_SCALE = 2.5, 3.0
data = rng.gamma(TRUE_SHAPE, TRUE_SCALE, size=2000)

# Method of moments, in closed form.
m, v = data.mean(), data.var(ddof=1)
mom_scale = v / m
mom_shape = m ** 2 / v
print(f"MoM: shape={mom_shape:.3f} scale={mom_scale:.3f}")   # near 2.5, 3.0

# Maximum likelihood, numerically, started from the MoM answer.
def nll(params, x):
    shape, scale = params
    if shape <= 0 or scale <= 0:
        return np.inf
    return -np.sum(stats.gamma.logpdf(x, a=shape, scale=scale))

fit = optimize.minimize(nll, x0=[mom_shape, mom_scale], args=(data,),
                        method="Nelder-Mead")
print(f"MLE: shape={fit.x[0]:.3f} scale={fit.x[1]:.3f}")

# Compare efficiency: repeat both on many small samples.
mom_est, mle_est = [], []
for _ in range(300):
    s = rng.gamma(TRUE_SHAPE, TRUE_SCALE, size=50)          # small sample
    mm, vv = s.mean(), s.var(ddof=1)
    mom_est.append(mm ** 2 / vv)
    mle_est.append(optimize.minimize(nll, x0=[mm**2/vv, vv/mm], args=(s,),
                                     method="Nelder-Mead").x[0])
print(f"MoM shape: mean={np.mean(mom_est):.3f} sd={np.std(mom_est):.3f}")
print(f"MLE shape: mean={np.mean(mle_est):.3f} sd={np.std(mle_est):.3f}")
# MLE has the smaller standard deviation, which is what "more efficient"
# means. At n=2000 the two agree closely; at n=50 the gap is visible.

Key Takeaways

• The method of moments equates theoretical moments to sample moments and solves for the parameters.
• It is frequently available in closed form, requiring no optimizer or likelihood function.
• It is less efficient than maximum likelihood and can return values outside the valid parameter space.
• Its best practical use is supplying starting values for a maximum likelihood fit.

🧪 Practice

1. Derive method-of-moments estimators for the two parameters of a uniform distribution on [a, b].
2. Compare method-of-moments and maximum likelihood estimates for a gamma at sample sizes 30 and 3000.
3. Interview: When would method of moments give an invalid estimate? (Hint: consider a parameter that must be positive and a sample whose moments imply otherwise.)

Confidence Intervals

A point estimate alone is dishonest, because it conveys no sense of how uncertain it is. A sample mean of 52.3 might be reliable to within 0.1 or within 20, and the point estimate looks identical in both cases. A confidence interval attaches that uncertainty.

A 95% confidence interval is constructed by a procedure with this property: if you repeated the whole study many times, 95% of the intervals produced would contain the true parameter. The confidence attaches to the procedure, not to any particular interval.

This is the part everyone gets wrong. Once computed, a specific interval either contains the parameter or does not; there is no probability left in it. Saying "there is a 95% probability the true mean lies in this interval" is a Bayesian statement about a credible interval, covered later in this chapter. The frequentist statement is about long-run coverage of the method.

The standard construction is the estimate plus or minus a critical value times the standard error. Use the t-distribution rather than the normal when the population standard deviation is unknown and estimated from the sample, which is essentially always. The difference matters below roughly 30 observations and becomes negligible above it.

  WHAT "95% CONFIDENCE" ACTUALLY MEANS

   true mu
      |
   ---+--------------------------------------------
      |  [------]                 interval 1  covers
      |      [------]             interval 2  covers
      |[------]                   interval 3  covers
      |             [------]      interval 4  MISSES
      |   [------]                interval 5  covers
      ...
   Across many repetitions, 95% of the intervals cover mu.
   Any SINGLE interval either does or does not. No probability remains.

  THE CONSTRUCTION

   estimate  +/-  critical value  *  standard error
     x-bar          t(0.975, df)        s / sqrt(n)
      ^                 ^                    ^
   your best      how many SEs wide    how noisy the
   guess          for the confidence   estimate is
                  level you chose

  WIDTH SHRINKS ONLY AS sqrt(n): to halve it, quadruple the data.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
TRUE_MU, TRUE_SIGMA = 100.0, 15.0

def t_interval(sample, confidence=0.95):
    n = len(sample)
    se = sample.std(ddof=1) / np.sqrt(n)             # standard error
    # t critical value with n-1 degrees of freedom, not 1.96, when n is small
    crit = stats.t.ppf((1 + confidence) / 2, df=n - 1)
    half = crit * se
    return sample.mean() - half, sample.mean() + half

sample = rng.normal(TRUE_MU, TRUE_SIGMA, size=30)
lo, hi = t_interval(sample)
print(f"mean={sample.mean():.2f}, 95% CI = ({lo:.2f}, {hi:.2f})")

# Verify the coverage claim by simulation -- this IS the definition.
covered = 0
for _ in range(10_000):
    s = rng.normal(TRUE_MU, TRUE_SIGMA, size=30)
    lo, hi = t_interval(s)
    covered += (lo <= TRUE_MU <= hi)
print(f"empirical coverage = {covered / 10_000:.3%}")     # ~95.0%

# Using 1.96 instead of the t critical value undercovers at small n.
def wrong_interval(sample, z=1.96):
    se = sample.std(ddof=1) / np.sqrt(len(sample))
    return sample.mean() - z * se, sample.mean() + z * se

for n in (5, 10, 30, 100):
    cov_t = cov_z = 0
    for _ in range(10_000):
        s = rng.normal(TRUE_MU, TRUE_SIGMA, size=n)
        lo, hi = t_interval(s);      cov_t += (lo <= TRUE_MU <= hi)
        lo, hi = wrong_interval(s);  cov_z += (lo <= TRUE_MU <= hi)
    print(f"n={n:>4}: t-based {cov_t/100:.1f}%   z-based {cov_z/100:.1f}%")
# n=   5: t-based 95.0%   z-based 87.6%   <- the z interval is too narrow
# n=  10: t-based 95.0%   z-based 91.9%
# n=  30: t-based 95.0%   z-based 93.9%
# n= 100: t-based 95.0%   z-based 94.8%   <- the gap has closed

# Width scales as 1/sqrt(n): quadrupling the sample halves the interval.
for n in (100, 400, 1600):
    s = rng.normal(TRUE_MU, TRUE_SIGMA, n)
    lo, hi = t_interval(s)
    print(f"n={n:>5}: width = {hi - lo:.3f}")

Key Takeaways

• Confidence is a property of the procedure's long-run coverage, not of any single computed interval.
• Use the t-distribution when the standard deviation is estimated from the sample, which is nearly always.
• Using 1.96 at small sample sizes produces intervals that are too narrow and undercover.
• Interval width shrinks with the square root of n, so halving it requires four times the data.

🧪 Practice

1. Build a 95% confidence interval for a sample mean and verify the coverage rate by simulation.
2. Compare t-based and z-based intervals at n = 5 and n = 200 and explain the difference in coverage.
3. Interview: A colleague says "there is a 95% chance the true mean is in this interval." Correct them. (Hint: ask what is random once the data has been collected.)

Bootstrap Resampling

Analytic confidence intervals require knowing the sampling distribution of your statistic. That is available for a mean, thanks to the central limit theorem, but not for a median, a 95th percentile, a correlation, a ratio of two medians, or the output of a complicated pipeline. The bootstrap supplies a confidence interval for essentially any statistic, using computation in place of mathematics.

The idea is to treat the sample as a stand-in for the population. Resample from it with replacement to the same size, producing a new pseudo-sample, and compute the statistic. Repeat thousands of times. The spread of those values approximates the sampling distribution of the statistic, from which an interval is read off directly as percentiles.

Sampling with replacement is essential. Without it, every resample would be the original sample reordered, and the statistic would never vary.

An analogy: you cannot rerun the study a thousand times, so you rerun it on the best available stand-in for the population, which is your sample. The approximation is good when the sample is representative, and inherits any bias when it is not.

  THE BOOTSTRAP LOOP

   original sample (n = 8)
   [ 3  7  2  9  4  6  8  5 ]
            |
            | resample WITH replacement, size 8, thousands of times
            v
   [ 7  2  7  9  9  3  5  2 ]  -> statistic = 5.50
   [ 3  3  4  6  8  8  9  7 ]  -> statistic = 6.00
   [ 9  5  2  4  4  6  7  3 ]  -> statistic = 5.00
   ...  (repeat 10,000 times)

   distribution of the bootstrap statistics
        |     ___
        |   .'   '.            the 2.5th and 97.5th percentiles
        |  /       \           of THIS distribution are the
        | /         \          95% confidence interval
        +--+-------+---
          lo       hi

  WHEN IT WORKS AND WHEN IT DOES NOT
   works:  means, medians, quantiles, correlations, ratios, model metrics
   fails:  the maximum or minimum (the extreme is bounded by the sample)
           very small samples (nothing to resample from)
           dependent data, unless you resample BLOCKS not points

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
# Skewed data where the median is the sensible summary and has no
# convenient analytic standard error.
data = rng.lognormal(mean=3.0, sigma=1.0, size=200)

def bootstrap_ci(sample, statistic, n_boot=10_000, confidence=0.95, seed=0):
    r = np.random.default_rng(seed)
    n = len(sample)
    # Resample indices WITH replacement, n_boot times, vectorized.
    idx = r.integers(0, n, size=(n_boot, n))
    stats_boot = statistic(sample[idx], axis=1)
    alpha = (1 - confidence) / 2
    return np.percentile(stats_boot, [100 * alpha, 100 * (1 - alpha)])

lo, hi = bootstrap_ci(data, np.median)
print(f"median = {np.median(data):.2f}, 95% CI = ({lo:.2f}, {hi:.2f})")

# Works for statistics with no standard formula at all.
lo, hi = bootstrap_ci(data, lambda x, axis: np.percentile(x, 90, axis=axis))
print(f"p90 = {np.percentile(data, 90):.2f}, 95% CI = ({lo:.2f}, {hi:.2f})")

# Validate against the analytic interval where one exists.
boot_lo, boot_hi = bootstrap_ci(data, np.mean)
se = data.std(ddof=1) / np.sqrt(len(data))
crit = stats.t.ppf(0.975, df=len(data) - 1)
print(f"bootstrap mean CI = ({boot_lo:.2f}, {boot_hi:.2f})")
print(f"analytic  mean CI = ({data.mean() - crit*se:.2f}, "
      f"{data.mean() + crit*se:.2f})")
# The two agree closely, which is the reassurance that the method works.

# Where it FAILS: the maximum can never exceed the observed maximum,
# so the bootstrap interval is one-sided and badly wrong.
lo, hi = bootstrap_ci(data, np.max)
print(f"observed max = {data.max():.1f}, bootstrap CI = ({lo:.1f}, {hi:.1f})")
# The upper bound EQUALS the observed maximum. The true population
# maximum is certainly larger, and the bootstrap cannot see that.

# SciPy provides a well-tested implementation with better interval types.
res = stats.bootstrap((data,), np.median, confidence_level=0.95,
                      n_resamples=10_000, random_state=0)
print(res.confidence_interval)

Key Takeaways

• The bootstrap approximates a sampling distribution by resampling the sample with replacement.
• It provides intervals for statistics that have no convenient analytic standard error.
• Resampling must be with replacement, or every resample would be identical to the original.
• It fails for extreme-value statistics and for dependent data unless you resample blocks rather than individual points.

🧪 Practice

1. Bootstrap a 95% confidence interval for the median of a skewed sample and compare its width to the interval for the mean.
2. Bootstrap the correlation between two variables and report the interval.
3. Interview: Why does the bootstrap fail for the maximum of a distribution? (Hint: ask what the largest value in any resample can possibly be.)

<a id="33-hypothesis-testing"></a>

3.3 Hypothesis Testing

Hypothesis testing is the formal procedure for deciding whether an observed difference is large enough to be worth believing, and most of its difficulty lies in understanding precisely what it does and does not claim.

Null and Alternative Hypotheses

Suppose a new checkout flow converts at 12.4% against the old flow's 11.8%. Is that a real improvement or the kind of wobble you would see from two identical flows? Hypothesis testing answers by setting up a deliberately skeptical framework.

The null hypothesis, written H0, is the position of no effect: the two flows convert identically, the coin is fair, the drug does nothing. The alternative hypothesis, H1, is what you suspect instead. The test asks a single question: if the null were true, how surprising would this data be?

The asymmetry is deliberate and important. You never prove the null. You either reject it, concluding the data is too improbable under it, or fail to reject it, concluding the data is compatible with it. "Fail to reject" is not "accept". Absence of evidence is not evidence of absence, and a test on 20 users failing to detect a 5% effect tells you almost nothing.

Tests are one-tailed if the alternative specifies a direction and two-tailed if it does not. Two-tailed is the honest default. Choosing one-tailed after seeing which way the data went doubles your effective false positive rate, and it is a recognized form of cheating.

An analogy: a criminal trial. The null is innocence, and it is presumed. The verdict is "guilty" or "not guilty", never "innocent". Acquittal means the evidence was insufficient, not that innocence was demonstrated.

  THE LOGIC OF THE TEST

   1. Assume H0 is true (no effect)
   2. Compute how extreme the observed data is under that assumption
   3. If it is extreme enough, reject H0
   4. Otherwise, fail to reject H0  -- NOT "H0 is true"

  ONE-TAILED vs TWO-TAILED at alpha = 0.05

   TWO-TAILED                        ONE-TAILED (right)
    reject        reject                        reject
     2.5%          2.5%                          5%
    \\\\|          |////                         |////////
   -----+----------+-----              ----------+-----------
      -1.96      +1.96                         +1.645
   alternative: "different"           alternative: "greater"

   Note the one-tailed threshold is LOWER, so it is easier to reject.
   That is legitimate only if the direction was committed to in advance.

  WHAT YOU MAY AND MAY NOT CONCLUDE
   p < alpha      -> "the data is unlikely under H0" -> reject H0
   p >= alpha     -> "the data is compatible with H0" -> say nothing more
   NEVER          -> "we proved there is no effect"

import numpy as np
from scipy import stats

# Two checkout flows. Is the difference real?
control_n, control_conv = 5000, 590        # 11.80%
variant_n, variant_conv = 5000, 620        # 12.40%

# H0: the two conversion rates are equal.
# H1: they differ (two-tailed -- we did not pre-commit to a direction).
count = np.array([control_conv, variant_conv])
nobs = np.array([control_n, variant_n])

# Pooled two-proportion z-test, written out so the logic is visible.
p_pool = count.sum() / nobs.sum()
se = np.sqrt(p_pool * (1 - p_pool) * (1 / control_n + 1 / variant_n))
z = (variant_conv / variant_n - control_conv / control_n) / se
p_two = 2 * (1 - stats.norm.cdf(abs(z)))          # two tails
p_one = 1 - stats.norm.cdf(z)                     # one tail, right

print(f"observed lift = {(620/5000 - 590/5000) * 100:.2f} pp")   # 0.60 pp
print(f"z = {z:.3f}")                                            # ~0.92
print(f"two-tailed p = {p_two:.4f}")                             # ~0.358
print(f"one-tailed p = {p_one:.4f}")                             # ~0.179
# Neither reaches 0.05. We FAIL TO REJECT. That is not evidence the
# variant is useless; the test simply cannot distinguish 0.6pp from
# noise at this sample size.

# The danger of switching tails after looking: under a true null, the
# post-hoc one-tailed test rejects about twice as often as advertised.
rng = np.random.default_rng(0)
false_positives = 0
for _ in range(20_000):
    a = rng.binomial(5000, 0.12)
    b = rng.binomial(5000, 0.12)                  # H0 is TRUE by construction
    pp = (a + b) / 10_000
    s = np.sqrt(pp * (1 - pp) * (2 / 5000))
    zz = (b / 5000 - a / 5000) / s
    # Cheat: pick the tail that favours whichever direction appeared.
    p_cheat = 1 - stats.norm.cdf(abs(zz))
    false_positives += (p_cheat < 0.05)
print(f"false positive rate when tail chosen post-hoc: "
      f"{false_positives / 20_000:.1%}")           # ~10%, not 5%

Key Takeaways

• The null hypothesis is the skeptical default of no effect, and the test asks how surprising the data would be if it held.
• You reject or fail to reject; you never prove the null.
• Failing to reject at a small sample size carries almost no information about whether an effect exists.
• Choose one-tailed or two-tailed before seeing the data; switching afterwards doubles the false positive rate.

🧪 Practice

1. State the null and alternative hypotheses for three business questions, including whether each should be one- or two-tailed.
2. Simulate a true null and confirm that a correctly specified test rejects about 5% of the time.
3. Interview: A test returns p = 0.35. Can you conclude the two versions perform identically? (Hint: consider what the result would look like if your sample were ten times larger.)

p-Values and Significance Levels

The p-value is the most used and most misunderstood number in statistics. Getting its definition exactly right prevents a long list of errors.

The p-value is the probability of observing data at least as extreme as what you got, assuming the null hypothesis is true. It is a statement about the data given the null, written P(data | H0).

It is emphatically not the probability that the null is true. That would be P(H0 | data), and Chapter 2 showed with Bayes' theorem that these two quantities can differ by orders of magnitude. Converting between them requires a prior, which the p-value does not have.

The significance level, alpha, is the threshold chosen in advance, conventionally 0.05. It is the false positive rate you are willing to accept: if the null is true, you will wrongly reject it 5% of the time, by construction.

The conventional 0.05 is arbitrary, a historical artifact rather than a law. Particle physics requires the equivalent of about 0.0000003. A low-stakes product experiment might reasonably use 0.10. The right threshold depends on the relative cost of the two error types.

The most consequential practical point: statistical significance is not practical importance. With a large enough sample, any difference however trivial becomes significant, because the standard error shrinks toward zero. Always report the effect size alongside the p-value.

  WHAT A p-VALUE IS

   distribution of the test statistic IF H0 IS TRUE

        |        _--_
        |      .'    '.
        |    .'        '.
        |  .'            '.
        +-'----------------'\\\\\\-----> test statistic
                             ^      ^
                          observed  |
                                    |
                     p-value = the shaded area:
                     P(a result this extreme or more | H0)

  WHAT IT IS NOT

   p = 0.03 does NOT mean:
     - there is a 3% chance the null is true      (that needs a prior)
     - there is a 97% chance the effect is real   (same problem)
     - the effect is large                        (says nothing of size)
     - the result will replicate 97% of the time  (unrelated question)

  SIGNIFICANCE vs IMPORTANCE

   n = 100      difference 0.5%   p = 0.72   not significant
   n = 10,000   difference 0.5%   p = 0.04   significant
   n = 1,000,000 difference 0.05% p = 0.001  significant AND pointless

   The effect size never changed meaningfully. Only the sample did.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# The same trivial effect at three sample sizes.
for n in (100, 10_000, 1_000_000):
    a = rng.normal(100.0, 15, n)
    b = rng.normal(100.3, 15, n)              # a real but tiny 0.3 difference
    t = stats.ttest_ind(a, b)
    # Cohen's d: the difference in units of the pooled standard deviation.
    pooled_sd = np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
    d = (b.mean() - a.mean()) / pooled_sd
    print(f"n={n:>9}: diff={b.mean()-a.mean():>6.3f}  p={t.pvalue:<10.2e} "
          f"Cohen's d={d:.4f}")
# n=      100: p large        d~0.02   not significant
# n=   10,000: p smaller      d~0.02   borderline
# n=1,000,000: p tiny         d~0.02   highly significant, still trivial
# The effect size d is stable near 0.02, which is negligible by any
# convention. Only the p-value moved, and it moved because of n alone.

# Under a true null, p-values are UNIFORMLY distributed. This is the
# cleanest way to internalize what a p-value is.
p_values = []
for _ in range(20_000):
    x = rng.normal(0, 1, 50)
    y = rng.normal(0, 1, 50)                  # H0 is true: same distribution
    p_values.append(stats.ttest_ind(x, y).pvalue)
p_values = np.array(p_values)

for lo, hi in [(0.0, 0.1), (0.4, 0.5), (0.9, 1.0)]:
    print(f"share of p in [{lo}, {hi}] = "
          f"{((p_values >= lo) & (p_values < hi)).mean():.3f}")   # each ~0.10
print(f"share below 0.05 = {(p_values < 0.05).mean():.4f}")       # ~0.05
# Exactly alpha, by construction. A "significant" result under a true
# null is not a malfunction; it is the advertised 5% error rate.

# Cohen's conventions for interpreting d, for reference:
#   0.2 small | 0.5 medium | 0.8 large

Key Takeaways

• A p-value is P(data this extreme | H0), never P(H0 | data).
• Alpha is the false positive rate you accept in advance, and 0.05 is a convention, not a law.
• Under a true null, p-values are uniformly distributed, so 5% of honest tests reject at alpha = 0.05.
• Statistical significance is driven by sample size; always report an effect size such as Cohen's d alongside it.

🧪 Practice

1. Simulate 10,000 tests under a true null and plot the histogram of p-values to confirm uniformity.
2. Take a fixed tiny effect and find the sample size at which it becomes significant, then compute Cohen's d at that size.
3. Interview: A result has p = 0.001 on 2 million rows. Should the company act on it? (Hint: ask for the effect size and what it is worth in revenue.)

Type I and Type II Errors

Any decision rule based on incomplete information makes two kinds of mistake, and they trade off against each other. Naming them makes the trade-off explicit and forces the question of which one costs more.

A Type I error is a false positive: rejecting a true null, declaring an effect that does not exist. Its probability is exactly alpha, because that is how alpha is defined.

A Type II error is a false negative: failing to reject a false null, missing an effect that is real. Its probability is called beta, and unlike alpha it is not chosen directly. It depends on the true effect size, the sample size, the variability in the data, and alpha.

The trade-off is unavoidable. Lowering alpha to reduce false positives makes the rejection threshold stricter, which increases beta. The only way to reduce both simultaneously is to collect more data.

Which error is worse is a business question, never a statistical one. Falsely declaring a drug effective harms patients; missing an effective drug also harms patients. Falsely flagging fraud annoys a customer; missing fraud costs money. The costs must be named before alpha is chosen.

                        REALITY
                 H0 true          H0 false
              +----------------+-----------------+
      reject  |  TYPE I ERROR  |   correct       |
       H0     |  false positive|   (power)       |
              |  probability   |   probability   |
  DECISION    |  = alpha       |   = 1 - beta    |
              +----------------+-----------------+
   fail to    |  correct       |  TYPE II ERROR  |
   reject H0  |  probability   |  false negative |
              |  = 1 - alpha   |  probability    |
              |                |  = beta         |
              +----------------+-----------------+

  THE TRADE-OFF, SHOWN AS OVERLAPPING DISTRIBUTIONS

    H0 true              H1 true
       _--_               _--_
     .'    '.           .'    '.
    /        \         /        \
   /       \\\\\\\\\\\/          \
  ------------+----------------------
              ^ threshold
        \\\\  = alpha (Type I)
        ////  = beta  (Type II)

  Slide the threshold RIGHT:  alpha falls, beta rises
  Slide the threshold LEFT:   alpha rises, beta falls
  Move the curves APART (bigger effect) or narrow them (more data):
    BOTH fall. That is the only free lunch, and it costs sample size.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

def error_rates(true_diff, n, alpha, trials=20_000):
    """Empirical Type I (if diff=0) or power (if diff>0)."""
    rejects = 0
    for _ in range(trials):
        a = rng.normal(0, 1, n)
        b = rng.normal(true_diff, 1, n)
        if stats.ttest_ind(a, b).pvalue < alpha:
            rejects += 1
    return rejects / trials

# Type I error equals alpha, by construction.
for alpha in (0.01, 0.05, 0.10):
    rate = error_rates(0.0, 50, alpha, trials=20_000)
    print(f"alpha={alpha:<5} -> empirical Type I rate = {rate:.4f}")
# 0.01 -> ~0.010 | 0.05 -> ~0.050 | 0.10 -> ~0.100

# The trade-off: stricter alpha buys fewer false positives and more misses.
print()
for alpha in (0.01, 0.05, 0.10):
    power = error_rates(0.5, 50, alpha, trials=20_000)
    print(f"alpha={alpha:<5} -> power = {power:.3f}, "
          f"Type II (beta) = {1 - power:.3f}")
# alpha=0.01 -> power ~0.46, beta ~0.54    strict: misses over half
# alpha=0.05 -> power ~0.70, beta ~0.30
# alpha=0.10 -> power ~0.80, beta ~0.20    lenient: catches more, risks more

# The only way to improve both: more data.
print()
for n in (50, 100, 200):
    power = error_rates(0.5, n, 0.01, trials=10_000)
    print(f"n={n:>4}, alpha=0.01 -> power = {power:.3f}")
# n=  50 -> ~0.46
# n= 100 -> ~0.76
# n= 200 -> ~0.97     strict alpha AND high power, bought with sample size

# Choosing alpha from costs rather than convention.
COST_FALSE_POSITIVE = 50_000      # launching a feature that does nothing
COST_FALSE_NEGATIVE = 5_000       # missing a small improvement
# When a false positive costs ten times a false negative, a stricter
# alpha than 0.05 is the rational choice, and the reverse when it is
# cheap to be wrong and expensive to miss out.

Key Takeaways

• Type I is a false positive with probability alpha; Type II is a false negative with probability beta.
• Alpha is chosen directly; beta is a consequence of effect size, sample size, variability, and alpha.
• Tightening alpha always raises beta, so the two can only be improved together by collecting more data.
• Which error is worse is determined by business costs, not by statistical convention.

🧪 Practice

1. Simulate tests under a true null at three alpha levels and confirm the Type I rate matches each alpha.
2. Hold the effect size fixed and show how beta falls as n increases.
3. Interview: For a cancer screening test, which error would you tolerate more and how does that change your threshold? (Hint: compare the cost of an unnecessary follow-up against a missed diagnosis.)

Statistical Power and Sample Size

Power is the probability of detecting an effect that genuinely exists, 1 - beta. It is the question every experiment should answer before it runs: if there is a real effect of the size we care about, how likely are we to see it?

Running an underpowered experiment is worse than running none. It consumes time and traffic, and then produces an ambiguous result that cannot distinguish "no effect" from "we could not tell". Worse, among the underpowered studies that do reach significance, the estimated effects are systematically inflated, because only the luckiest, largest-looking samples cross the threshold. This is the winner's curse, and it is a major driver of the replication crisis.

Power depends on four quantities, and fixing any three determines the fourth. That is what makes sample size calculation possible. The convention is 80% power, meaning you accept a one-in-five chance of missing a real effect.

The most important input is the minimum detectable effect: the smallest difference that would change your decision. It is a business judgement, not a statistical one, and it must be set before the experiment.

  THE FOUR INTERLOCKING QUANTITIES

        alpha  <-->  power  <-->  effect size  <-->  sample size

   Fix any three, and the fourth is determined. The usual question is:
   "given alpha = 0.05, power = 0.80, and the smallest effect worth
    detecting, how many observations do I need?"

  POWER AS A FUNCTION OF SAMPLE SIZE (effect = 0.3 sd, alpha = 0.05)

   power
    1.0 |                        ...--------
        |                ..---''
    0.8 |- - - - - - -.-'- - - - - - - - - -   <- the 80% convention
        |          .-'
    0.5 |        .'
        |     .-'
    0.0 +---'--------------------------------> n per group
        0   100  200  300  400  500  600

   Note the flattening: going from 80% to 95% power costs far more
   than going from 50% to 80%.

  THE WINNER'S CURSE

   true effect = 0.2, underpowered study (power = 0.20)
   Among the 20% of runs that reach significance, the AVERAGE
   estimated effect is roughly 0.5 -- inflated two to three times,
   because only unusually large sample fluctuations clear the bar.

import numpy as np
from scipy import stats
from statsmodels.stats.power import TTestIndPower, NormalIndPower
from statsmodels.stats.proportion import proportion_effectsize

# How many users per group to detect a 0.3 standard deviation shift?
analysis = TTestIndPower()
n = analysis.solve_power(effect_size=0.3, alpha=0.05, power=0.80,
                         alternative="two-sided")
print(f"n per group = {np.ceil(n):.0f}")          # 176

# The cost of detecting smaller effects grows as 1/effect^2.
for d in (0.8, 0.5, 0.3, 0.1):
    n = analysis.solve_power(effect_size=d, alpha=0.05, power=0.80)
    print(f"effect {d}: n per group = {np.ceil(n):>6.0f}")
# effect 0.8: n =     26
# effect 0.5: n =     64
# effect 0.3: n =    176
# effect 0.1: n =   1571      halving the effect QUADRUPLES the sample

# For conversion rates, convert to an effect size first.
baseline, target = 0.12, 0.13                    # a 1pp absolute lift
effect = proportion_effectsize(target, baseline)
n = NormalIndPower().solve_power(effect_size=effect, alpha=0.05, power=0.80)
print(f"\n12% -> 13% needs n = {np.ceil(n):,.0f} per group")   # ~8,700

# Demonstrate the winner's curse directly.
rng = np.random.default_rng(0)
TRUE_EFFECT = 0.2
significant_estimates = []
all_estimates = []
for _ in range(20_000):
    a = rng.normal(0, 1, 20)                     # deliberately underpowered
    b = rng.normal(TRUE_EFFECT, 1, 20)
    est = b.mean() - a.mean()
    all_estimates.append(est)
    # Count only significant results pointing the RIGHT way, which is
    # what a team would actually act on and report.
    if stats.ttest_ind(a, b).pvalue < 0.05 and est > 0:
        significant_estimates.append(est)

print(f"\ntrue effect                = {TRUE_EFFECT}")
print(f"mean of ALL estimates      = {np.mean(all_estimates):.3f}")
print(f"mean of SIGNIFICANT ones   = {np.mean(significant_estimates):.3f}")
print(f"share significant, correct sign = "
      f"{len(significant_estimates)/20_000:.3f}")
# All estimates average ~0.20, correctly recovering the truth. But the
# significant subset averages ~0.77, roughly four times the real effect,
# and only about 9% of runs land there. If you publish or act only on
# significant results, you systematically overstate what you found.

Key Takeaways

• Power is the probability of detecting a real effect, and 80% is the usual convention.
• Alpha, power, effect size, and sample size are interlocked; fixing three determines the fourth.
• Required sample size scales with the inverse square of the effect, so halving the detectable effect quadruples the cost.
• Underpowered studies that reach significance systematically overstate the effect, which is the winner's curse.

🧪 Practice

1. Compute the sample size needed to detect a conversion lift from 5% to 5.5% at 80% power.
2. Plot power against sample size for a fixed effect and mark where the curve flattens.
3. Interview: Your A/B test was underpowered but came back significant with a huge effect. How do you interpret it? (Hint: think about which samples were able to clear the threshold.)

t-Tests

The t-test compares means and is the most frequently used test in applied work. It exists because of a specific problem: when the population standard deviation is unknown and must be estimated from a small sample, the estimate itself is uncertain, and pretending otherwise produces overconfident results.

William Gosset, publishing as "Student" while working at Guinness, derived the distribution that accounts for this. The t-distribution resembles a normal but has heavier tails, and those heavier tails are exactly the extra caution the extra uncertainty demands. As the sample grows, the estimate firms up and the t-distribution converges to the normal.

Three variants cover the common cases. The one-sample t-test compares a sample mean to a fixed reference value. The independent two-sample test compares two separate groups. The paired test compares two measurements on the same units, such as before and after, and is far more powerful when applicable because it eliminates between-subject variation.

The assumptions are approximate normality of the sampling distribution of the mean, which the central limit theorem usually supplies, and independence, which it does not. Independence is the assumption that actually gets violated. For unequal variances, always use Welch's version, which is the default in SciPy and costs essentially nothing when variances are in fact equal.

  t-DISTRIBUTION vs NORMAL

        normal       .-'''-.
        t (df=5)    /:     :\        the t has FATTER TAILS,
                   / :     : \       so a more extreme observation
                  /  :     :  \      is needed to reach significance
        .........'   :     :   '........
      ----------------------------------
                  ^         ^
              t critical  z critical
              (further out)

   df=1   very heavy tails      df=30  nearly normal
   df=5   noticeably heavier    df=inf exactly normal

  WHICH t-TEST

   one sample   : is this mean different from a known value?
                  "is average latency above our 200ms SLA?"

   independent  : do these two separate groups differ?
                  "do control and variant convert differently?"

   paired       : did the same units change?
                  "did latency drop after the optimization, per endpoint?"
                  much more powerful: removes per-unit differences

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# ONE SAMPLE: is mean latency above the 200ms service level agreement?
latency = rng.normal(207, 25, size=40)
res = stats.ttest_1samp(latency, popmean=200, alternative="greater")
print(f"mean={latency.mean():.1f}  t={res.statistic:.3f}  p={res.pvalue:.4f}")

# INDEPENDENT TWO SAMPLE: always prefer Welch (equal_var=False).
control = rng.normal(100, 15, size=60)
variant = rng.normal(106, 22, size=55)          # different mean AND variance
student = stats.ttest_ind(control, variant, equal_var=True)
welch = stats.ttest_ind(control, variant, equal_var=False)
print(f"Student p={student.pvalue:.4f}   Welch p={welch.pvalue:.4f}")
# Welch adjusts the degrees of freedom for unequal variances. When the
# variances ARE equal it loses almost nothing, so it is the safe default.

# PAIRED: the same endpoints measured before and after an optimization.
endpoints = rng.normal(300, 120, size=25)        # endpoints differ a LOT
before = endpoints + rng.normal(0, 8, 25)
after = endpoints - 12 + rng.normal(0, 8, 25)    # a consistent 12ms gain

unpaired = stats.ttest_ind(before, after)
paired = stats.ttest_rel(before, after)
print(f"\nunpaired p = {unpaired.pvalue:.4f}   paired p = {paired.pvalue:.6f}")
# The unpaired test is swamped by the huge variation BETWEEN endpoints
# and may miss the effect entirely. The paired test looks only at the
# within-endpoint differences, where the 12ms gain is unmistakable.
print(f"per-endpoint differences: mean={np.mean(before-after):.2f}, "
      f"sd={np.std(before-after, ddof=1):.2f}")

# Why Welch is the default: Student's test breaks with unequal variances
# and unequal group sizes.
false_pos_student = false_pos_welch = 0
for _ in range(10_000):
    a = rng.normal(0, 5, 10)                     # SMALL group, LARGE variance
    b = rng.normal(0, 1, 100)                    # large group, small variance
    false_pos_student += stats.ttest_ind(a, b, equal_var=True).pvalue < 0.05
    false_pos_welch += stats.ttest_ind(a, b, equal_var=False).pvalue < 0.05
print(f"\nType I rate -- Student: {false_pos_student/10_000:.3f}  "
      f"Welch: {false_pos_welch/10_000:.3f}")
# Student rejects far more often than the advertised 5%; Welch holds at
# about 0.05. The pooled variance is dominated by the large group, so it
# badly understates the uncertainty in the small group's mean.
#
# The direction matters: pair a SMALL group with a LARGE variance and
# Student's test is too liberal, as here. Reverse it, giving the small
# group the small variance, and it becomes too conservative instead.
# Welch is correct either way, which is why it is the default.

Key Takeaways

• The t-distribution has heavier tails than the normal to account for estimating the standard deviation from the sample.
• Use Welch's test by default; it handles unequal variances and costs almost nothing when they are equal.
• Paired tests remove between-unit variation and are much more powerful when the design allows them.
• The assumption that actually matters is independence, not normality of the raw data.

🧪 Practice

1. Run one-sample, independent, and paired t-tests on suitable data and interpret each result.
2. Construct data where the paired test is significant and the unpaired test is not, and explain why.
3. Interview: When is Student's t-test inappropriate and what do you use instead? (Hint: consider two groups with very different spreads and very different sizes.)

Chi-Square Tests

The t-test compares means of numeric data. When both variables are categorical, means are undefined and a different tool is needed. Chi-square tests compare observed counts against the counts expected under a hypothesis, building directly on the contingency tables from earlier in this chapter.

Two variants cover most uses. The test of independence asks whether two categorical variables are associated, comparing the observed cross-tab to what independence would produce. The goodness-of-fit test asks whether one categorical variable follows a specified distribution, such as whether a die is fair or whether traffic is evenly split across variants.

Both compute the same statistic: the sum over cells of the squared difference between observed and expected, divided by expected. Dividing by expected scales each discrepancy by how surprising it is, since a gap of 10 matters enormously when 5 were expected and hardly at all when 5000 were.

The assumptions are counts rather than percentages, independent observations, and expected counts of at least 5 in most cells. Small expected counts break the approximation, and Fisher's exact test should be used instead.

As always with large samples, significance arrives easily, so report Cramer's V or another effect size alongside.

  THE CHI-SQUARE STATISTIC

   chi2 = sum over cells of  (observed - expected)^2 / expected
                                                       \_________/
                             scaling by expected is what makes a
                             discrepancy of 10 mean something
                             different at E=5 than at E=5000

  TEST OF INDEPENDENCE            GOODNESS OF FIT

   observed cross-tab              observed counts vs a claimed
   vs expected under                distribution
   independence
                                    die face: 1   2   3   4   5   6
        A    B                      observed: 8  14  11  22  18  27
   X   [  ][  ]                     expected:16.7 ... (all equal)
   Y   [  ][  ]                     H0: the die is fair

   H0: the row and column
       variables are independent

  ASSUMPTION CHECK
   expected count >= 5 in at least 80% of cells, and never below 1.
   Otherwise use Fisher's exact test (2x2) or pool sparse categories.

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

# TEST OF INDEPENDENCE: does device type relate to churn?
observed = np.array([
    [200, 50],      # desktop: retained, churned
    [180, 70],      # mobile
    [120, 80],      # tablet
])
result = stats.chi2_contingency(observed)
print(f"chi2 = {result.statistic:.2f}, df = {result.dof}, "
      f"p = {result.pvalue:.2e}")
print("expected under independence:\n", result.expected_freq.round(1))

# Always check the assumption BEFORE trusting the p-value.
small_expected = (result.expected_freq < 5).sum()
print(f"cells with expected < 5: {small_expected}")      # 0, so we are fine

# Effect size, because chi-square grows with n regardless of strength.
n = observed.sum()
v = np.sqrt(result.statistic / (n * (min(observed.shape) - 1)))
print(f"Cramer's V = {v:.3f}")        # ~0.16, a weak-to-moderate association

# Standardized residuals show WHERE the association lives.
residuals = (observed - result.expected_freq) / np.sqrt(result.expected_freq)
print("standardized residuals:\n", residuals.round(2))
# Values beyond about +/- 2 mark the cells driving the result. Here the
# tablet-churn cell is well above expectation and desktop-churn below.

# GOODNESS OF FIT: is this die fair?
rolls = np.array([8, 14, 11, 22, 18, 27])          # 100 rolls
gof = stats.chisquare(rolls)                       # expects uniform by default
print(f"\ndie test: chi2 = {gof.statistic:.2f}, p = {gof.pvalue:.4f}")
# p < 0.05 -> reject fairness. Face 6 appears far more than 16.7 times.

# GOODNESS OF FIT against a NON-uniform expectation: traffic allocation.
observed_split = np.array([4900, 2600, 2500])
intended_split = np.array([0.50, 0.25, 0.25]) * observed_split.sum()
srm = stats.chisquare(observed_split, f_exp=intended_split)
print(f"sample ratio mismatch: p = {srm.pvalue:.4f}")
# This exact check is standard practice in A/B testing. A significant
# result means the randomization is broken, and the experiment's other
# results cannot be trusted until it is explained.

# When expected counts are small, switch to Fisher's exact test.
sparse = np.array([[8, 2], [1, 5]])
print(f"\nchi2 p = {stats.chi2_contingency(sparse).pvalue:.4f} (unreliable)")
print(f"Fisher p = {stats.fisher_exact(sparse).pvalue:.4f} (correct here)")

Key Takeaways

• Chi-square tests compare observed counts to expected counts and require raw counts, never percentages.
• The test of independence checks association between two categorical variables; goodness of fit checks one against a claimed distribution.
• Expected counts below 5 break the approximation, so use Fisher's exact test or pool categories.
• Standardized residuals identify which cells drive a significant result, which the overall p-value cannot.

🧪 Practice

1. Test whether two categorical columns in a dataset are independent and report Cramer's V alongside the p-value.
2. Run a sample ratio mismatch check on an A/B test's traffic split.
3. Interview: A chi-square test on 500,000 rows returns p below 1e-10. What do you report to the business? (Hint: significance is guaranteed at that size, so what actually informs the decision?)

ANOVA

A t-test compares two groups. With four groups, you might be tempted to run all six pairwise t-tests, and that is exactly the mistake ANOVA exists to prevent. Each test carries a 5% false positive risk, and six tests carry roughly a 26% chance of at least one spurious result.

Analysis of variance tests all groups at once with a single question: are these group means all equal? It works by decomposing total variation into two parts. Between-group variation reflects how far the group means sit from the overall mean. Within-group variation reflects the spread inside each group. The F-statistic is their ratio.

The logic is intuitive. If groups genuinely differ, the spread between group means will be large relative to the noise inside each group, and F will be large. If the groups are identical, both quantities estimate the same noise and F sits near 1.

ANOVA is an omnibus test: a significant result says "at least one group differs" without saying which. Identifying the specific pairs requires a follow-up post-hoc test such as Tukey's HSD, which builds in the correction for multiple comparisons.

Assumptions are independence, approximately normal residuals, and roughly equal variances across groups. When variances differ badly, use Welch's ANOVA; when normality fails badly, use the Kruskal-Wallis test.

  THE F-RATIO: SIGNAL OVER NOISE

              variation BETWEEN group means
    F  =  ---------------------------------------
              variation WITHIN groups (noise)

  LARGE F (groups differ)          SMALL F (groups do not)

   A  ..                            A  ...........
   B      ..                        B  ...........
   C          ..                    C  ...........
      +----------+                     +----------+
   tight within, spread between     wide within, means nearly identical
   F is large                       F is near 1

  WHY NOT SIX t-TESTS

   comparisons  P(at least one false positive) at alpha = 0.05
        1                    5.0%
        3                   14.3%
        6                   26.5%      <- four groups
       10                   40.1%      <- five groups

   ANOVA asks ONE question and holds the overall rate at 5%.

  THE WORKFLOW
   ANOVA significant?  --no-->  stop, no evidence of any difference
          |
         yes
          v
   post-hoc test (Tukey HSD) -> which specific pairs differ

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# Four landing pages, one of which is genuinely better.
groups = {
    "A": rng.normal(100, 15, 50),
    "B": rng.normal(102, 15, 50),
    "C": rng.normal(101, 15, 50),
    "D": rng.normal(112, 15, 50),      # the real winner
}

f_stat, p = stats.f_oneway(*groups.values())
print(f"ANOVA: F = {f_stat:.3f}, p = {p:.2e}")
# Significant -> at least one group differs, but WHICH is unanswered.

# Why not pairwise t-tests: demonstrate the inflated error rate.
false_positive_runs = 0
for _ in range(5_000):
    # Four IDENTICAL groups: every rejection is a false positive.
    gs = [rng.normal(100, 15, 50) for _ in range(4)]
    any_sig = False
    for i in range(4):
        for j in range(i + 1, 4):
            if stats.ttest_ind(gs[i], gs[j]).pvalue < 0.05:
                any_sig = True
    false_positive_runs += any_sig
print(f"pairwise t-tests, 4 identical groups: "
      f"{false_positive_runs / 5_000:.1%} of runs had a false positive")
# ~20%, four times the nominal 5%. This is below the 26.5% that
# 1 - 0.95^6 predicts, because the six comparisons share groups and are
# therefore correlated rather than independent.

# Post-hoc: which pairs actually differ, with the correction built in.
from statsmodels.stats.multicomp import pairwise_tukeyhsd
values = np.concatenate(list(groups.values()))
labels = np.repeat(list(groups.keys()), 50)
print(pairwise_tukeyhsd(values, labels, alpha=0.05))
# Tukey reports each pair with an adjusted p-value and confidence
# interval. Expect D vs A, D vs B, D vs C to be flagged and the rest not.

# Assumption checks before trusting any of it.
print(f"\nLevene (equal variances): p = "
      f"{stats.levene(*groups.values()).pvalue:.3f}")     # want p > 0.05
for name, g in groups.items():
    print(f"  {name}: Shapiro p = {stats.shapiro(g).pvalue:.3f}")

# Robust alternatives when assumptions fail.
print(f"\nKruskal-Wallis (no normality needed): "
      f"p = {stats.kruskal(*groups.values()).pvalue:.2e}")

Key Takeaways

• ANOVA tests whether any of several group means differ, holding the overall false positive rate at alpha.
• Running all pairwise t-tests inflates the error rate sharply: four groups give roughly a 20% chance of a spurious result.
• The F-statistic is between-group variation over within-group variation, and sits near 1 when groups are identical.
• A significant ANOVA requires a post-hoc test such as Tukey's HSD to identify which pairs differ.

🧪 Practice

1. Run a one-way ANOVA on four groups and follow it with Tukey's HSD, interpreting both outputs.
2. Simulate identical groups and measure how often at least one pairwise t-test is significant.
3. Interview: ANOVA returns p = 0.03. What can and cannot you tell the business? (Hint: the test's alternative hypothesis is quite weak.)

Non-Parametric Tests

Every test so far assumes something about the shape of the underlying distribution. Non-parametric tests drop those assumptions, buying robustness at the cost of some statistical power. They are the right choice for ordinal data, small samples from unknown distributions, and data with outliers severe enough to distort means.

The trick most of them use is to replace values with their ranks. Ranks discard magnitude and keep only order, which makes the test immune to outliers and to the distribution's shape, since the ranks of any continuous sample are just a permutation of 1 to n.

The main correspondences are worth memorizing. The Mann-Whitney U test replaces the independent two-sample t-test. The Wilcoxon signed-rank test replaces the paired t-test. Kruskal-Wallis replaces one-way ANOVA. Spearman's correlation replaces Pearson's.

A common misconception is that these tests compare medians. In general they compare distributions, and they answer "is one sample systematically larger than the other?" They compare medians only under the extra assumption that the two distributions have the same shape.

The cost is modest. When normality genuinely holds, Mann-Whitney retains about 95% of the t-test's power, which is a small price for immunity to assumption failure.

Parametric test	Non-parametric equivalent	Compares
One-sample t-test	Wilcoxon signed-rank	Location vs a value
Independent t-test	Mann-Whitney U	Two distributions
Paired t-test	Wilcoxon signed-rank	Paired differences
One-way ANOVA	Kruskal-Wallis	Several groups
Pearson correlation	Spearman correlation	Monotonic relation

  HOW RANK TESTS NEUTRALIZE OUTLIERS

   raw values   A: 12  15  18  21         B: 14  17  20  9000
   the 9000 dominates any mean-based comparison

   pooled ranks
     9000 -> 8     21 -> 6     20 -> 5     18 -> 4
       17 -> 3     15 -> 2     14 -> 1     12 -> 0  (0-indexed idea)

   A ranks: 1, 3, 5, 7        B ranks: 2, 4, 6, 8
   The monstrous outlier contributes exactly one rank, no more than
   any other observation. That is the entire mechanism.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# Heavily skewed data with outliers: session durations in seconds.
a = np.concatenate([rng.exponential(30, 48), [4000, 5200]])
b = np.concatenate([rng.exponential(45, 48), [3800, 6100]])

t_res = stats.ttest_ind(a, b, equal_var=False)
u_res = stats.mannwhitneyu(a, b, alternative="two-sided")
print(f"means:   A={a.mean():.1f}  B={b.mean():.1f}")
print(f"medians: A={np.median(a):.1f}  B={np.median(b):.1f}")
print(f"Welch t p = {t_res.pvalue:.4f}")
print(f"Mann-Whitney p = {u_res.pvalue:.4f}")
# The means are dominated by four extreme sessions and the t-test is
# swamped. Mann-Whitney sees the consistent rank shift and detects it.

# Paired version: Wilcoxon signed-rank on before/after measurements.
before = rng.exponential(100, 30)
after = before * rng.uniform(0.6, 0.95, 30)        # a consistent improvement
print(f"\npaired t p    = {stats.ttest_rel(before, after).pvalue:.5f}")
print(f"Wilcoxon p    = {stats.wilcoxon(before, after).pvalue:.5f}")

# Several groups without normality: Kruskal-Wallis.
gs = [rng.exponential(scale, 40) for scale in (20, 22, 35)]
print(f"\nANOVA p        = {stats.f_oneway(*gs).pvalue:.4f}")
print(f"Kruskal p      = {stats.kruskal(*gs).pvalue:.4f}")

# The cost of robustness, measured. When normality DOES hold, how much
# power does Mann-Whitney give up against the t-test?
t_wins = u_wins = 0
for _ in range(5_000):
    x = rng.normal(0, 1, 30)
    y = rng.normal(0.6, 1, 30)                     # a real effect
    t_wins += stats.ttest_ind(x, y).pvalue < 0.05
    u_wins += stats.mannwhitneyu(x, y).pvalue < 0.05
print(f"\npower on clean normal data -- t: {t_wins/5000:.3f}  "
      f"Mann-Whitney: {u_wins/5000:.3f}")
# The rank test retains roughly 95% of the t-test's power. Cheap
# insurance against assumptions you cannot verify.

Key Takeaways

• Non-parametric tests replace values with ranks, which removes both distributional assumptions and outlier sensitivity.
• Each parametric test has a standard rank-based counterpart worth memorizing.
• They compare distributions generally, and compare medians only when the distributions share a shape.
• On genuinely normal data they retain about 95% of the parametric test's power, making them cheap insurance.

🧪 Practice

1. Compare a t-test and Mann-Whitney on skewed data with outliers and explain the divergence.
2. Run Wilcoxon signed-rank on paired measurements and compare to the paired t-test.
3. Interview: When would you choose Mann-Whitney over a t-test? (Hint: consider what you can verify about the distribution at n = 15.)

Multiple Comparisons Correction

If you run one test at alpha = 0.05, you accept a 5% false positive rate. If you run 100 tests on data with no real effects, you should expect about five significant results. Nothing has malfunctioned. This is the advertised rate, and it is why testing many hypotheses without adjustment produces confident nonsense.

The problem appears everywhere in data science, usually without anyone labelling it: testing a treatment across 20 user segments, screening 10,000 genes, comparing every pair of groups, checking an A/B test daily for a month, or trying many features until one looks predictive.

Two families of correction exist, and they control different things. The family-wise error rate is the probability of any false positive. Bonferroni controls it by dividing alpha by the number of tests, which is simple and extremely conservative. Holm's method achieves the same guarantee with more power and should generally be preferred to plain Bonferroni.

The false discovery rate controls the expected proportion of false positives among the results you call significant. Benjamini-Hochberg is the standard method. This is usually the better choice in exploratory data science, where tolerating a known fraction of false leads is sensible and demanding zero is paralyzing.

  THE PROBLEM, QUANTIFIED

   P(at least one false positive) = 1 - (1 - alpha)^m

   tests m:    1      5      10     20     100
   P(error): 5.0%  22.6%  40.1%  64.2%   99.4%

   At 100 tests you are essentially GUARANTEED a false positive.

  TWO DIFFERENT GUARANTEES

   FAMILY-WISE ERROR RATE (FWER)
     "P(any false positive among all tests) <= 0.05"
     strict; use when a single false claim is costly
     methods: Bonferroni (simple), Holm (uniformly better)

   FALSE DISCOVERY RATE (FDR)
     "expected share of false positives AMONG REJECTIONS <= 0.05"
     permissive; use for screening and exploration
     method: Benjamini-Hochberg

   Of 100 rejections, FDR at 0.05 expects about 5 to be wrong.
   FWER at 0.05 expects a 5% chance that even ONE is wrong.

  BENJAMINI-HOCHBERG PROCEDURE
   1. sort the m p-values ascending
   2. find the largest k with  p(k) <= (k / m) * alpha
   3. reject hypotheses 1 through k

import numpy as np
from scipy import stats
from statsmodels.stats.multitest import multipletests

rng = np.random.default_rng(0)

# Test a treatment across 20 segments where it does NOTHING anywhere.
p_values = []
for _ in range(20):
    a = rng.normal(0, 1, 100)
    b = rng.normal(0, 1, 100)            # no effect, by construction
    p_values.append(stats.ttest_ind(a, b).pvalue)
p_values = np.array(p_values)

print(f"raw significant at 0.05: {(p_values < 0.05).sum()} of 20")
# Typically 1 or 2 "discoveries", all of them false.

for method, label in [("bonferroni", "Bonferroni"),
                      ("holm", "Holm"),
                      ("fdr_bh", "Benjamini-Hochberg")]:
    reject, adjusted, _, _ = multipletests(p_values, alpha=0.05, method=method)
    print(f"{label:<20} significant: {reject.sum()}")
# All three correctly report 0 discoveries.

# Now a realistic mix: 100 tests, 10 with a genuine effect.
p_mixed = []
truth = []
for i in range(100):
    effect = 0.8 if i < 10 else 0.0           # first ten are real
    a = rng.normal(0, 1, 50)
    b = rng.normal(effect, 1, 50)
    p_mixed.append(stats.ttest_ind(a, b).pvalue)
    truth.append(effect > 0)
p_mixed, truth = np.array(p_mixed), np.array(truth)

print(f"\n{'method':<22}{'found':>7}{'true':>7}{'false':>7}")
for method, label in [("none", "no correction"),
                      ("bonferroni", "Bonferroni"),
                      ("holm", "Holm"),
                      ("fdr_bh", "Benjamini-Hochberg")]:
    if method == "none":
        reject = p_mixed < 0.05
    else:
        reject, _, _, _ = multipletests(p_mixed, alpha=0.05, method=method)
    tp = (reject & truth).sum()
    fp = (reject & ~truth).sum()
    print(f"{label:<22}{reject.sum():>7}{tp:>7}{fp:>7}")
# no correction      ~14     ~10      ~5    many false leads
# Bonferroni          ~7      ~7       0    strict, misses real effects
# Holm                ~7      ~7       0    same guarantee, more power
# Benjamini-Hochberg ~10      ~9      ~1    best balance for screening

# Bonferroni by hand, to show there is no magic in it.
print(f"\nBonferroni threshold = {0.05 / 100:.5f}")
print(f"tests below it: {(p_mixed < 0.05 / 100).sum()}")

Key Takeaways

• Running m tests at alpha = 0.05 makes at least one false positive nearly certain once m reaches about 100.
• Family-wise methods control the chance of any false positive; Holm dominates plain Bonferroni and should replace it.
• False discovery rate control bounds the proportion of false positives among rejections and suits exploratory screening.
• Checking an experiment repeatedly over time is a multiple comparisons problem, even though it looks like one test.

🧪 Practice

1. Run 50 tests on pure noise and count significant results before and after Benjamini-Hochberg correction.
2. Implement the Benjamini-Hochberg procedure by hand and verify it against multipletests.
3. Interview: An analyst tested 30 segments and reports the two that came back significant. What do you ask? (Hint: work out how many significant results 30 null tests would produce on their own.)

<a id="34-bayesian-statistics"></a>

3.4 Bayesian Statistics

Bayesian statistics treats parameters as uncertain quantities with distributions rather than fixed unknowns, which lets you state directly what everyone wants to know: how probable each possible value is, given the data.

Priors, Likelihoods, Posteriors

The frequentist machinery of the previous section answers a slightly awkward question. A confidence interval tells you about the long-run behaviour of a procedure, when what you wanted was a probability statement about the parameter. Bayesian inference provides that directly, at the cost of requiring you to state your beliefs before seeing the data.

The framework has three components, connected by Bayes' theorem from Chapter 2. The prior is your belief about the parameter before observing data, expressed as a probability distribution. The likelihood is how probable the observed data would be under each possible parameter value. The posterior is the updated belief after combining them.

The mechanism is simple to describe: the posterior is proportional to the likelihood times the prior. Parameter values that both fit the data well and were plausible beforehand come out with high posterior probability.

The prior is the controversial part, and the objection is that it injects subjectivity. Two responses matter in practice. First, with enough data the likelihood dominates and the prior's influence vanishes, so reasonable people with different priors converge. Second, a prior is an explicit, auditable assumption, whereas frequentist analyses carry implicit assumptions that are harder to inspect.

An analogy: a doctor seeing a patient with a cough does not start from scratch. Prior knowledge of what is circulating this season plus the specific symptoms yields a diagnosis. Ignoring the prior would mean treating a common cold and a rare tropical disease as equally likely.

  POSTERIOR = LIKELIHOOD x PRIOR  (then normalized)

   PRIOR                LIKELIHOOD            POSTERIOR
   (belief before)      (what data says)      (updated belief)

      _--_                     /\                  /\
    .'    '.                  /  \                /  \
   /        \        x       /    \      =       /    \
  ----------------         --------------      -------------
   broad, uncertain         data is             narrower, and
                            informative         pulled toward
                                                the data

  THE PRIOR'S INFLUENCE FADES WITH DATA

   n = 5            n = 50           n = 500
    _--_              /\               |
  .'    '.           /  \              |
 /        \         /    \            /|\
 ----------        --------          --+--
 prior still       data taking       prior irrelevant;
 dominates         over              likelihood wins

  Two analysts with different priors converge as evidence accumulates.
  If they do not, the disagreement is about the data, not the prior.

import numpy as np
from scipy import stats

# Estimating a conversion rate. The Beta distribution is the natural
# prior for a probability, since it lives on [0, 1].

# Three analysts with different starting beliefs.
priors = {
    "uninformative": (1, 1),        # Beta(1,1) is uniform: anything goes
    "optimistic":    (8, 2),        # believes the rate is around 80%
    "skeptical":     (2, 18),       # believes the rate is around 10%
}

observed_successes, observed_trials = 30, 100      # 30% observed

print(f"{'prior':<15}{'prior mean':>12}{'posterior mean':>16}{'95% CI':>22}")
for name, (a, b) in priors.items():
    # Beta is CONJUGATE to the binomial: the posterior is just
    # Beta(a + successes, b + failures). No integration required.
    post_a = a + observed_successes
    post_b = b + (observed_trials - observed_successes)
    post = stats.beta(post_a, post_b)
    lo, hi = post.interval(0.95)
    print(f"{name:<15}{a/(a+b):>12.3f}{post.mean():>16.3f}"
          f"{f'({lo:.3f}, {hi:.3f})':>22}")
# uninformative   0.500          0.304        (0.221, 0.394)
# optimistic      0.800          0.345        (0.257, 0.438)
# skeptical       0.100          0.267        (0.187, 0.354)
# The priors disagreed wildly (0.10 vs 0.80) and the posteriors are
# already close after only 100 observations.

# With more data they converge completely.
print()
for n in (10, 100, 1000, 10_000):
    s = int(0.30 * n)                              # same 30% rate throughout
    means = []
    for name, (a, b) in priors.items():
        means.append(stats.beta(a + s, b + n - s).mean())
    print(f"n={n:>6}: posterior means = "
          f"{[round(m, 3) for m in means]}  spread={max(means)-min(means):.3f}")
# n=    10: spread ~0.38   the prior is doing most of the work
# n=   100: spread ~0.08
# n=  1000: spread ~0.008
# n= 10000: spread ~0.001  the prior is now irrelevant

# The posterior answers the question people actually ask.
post = stats.beta(1 + 30, 1 + 70)
print(f"\nP(rate > 0.25 | data) = {1 - post.cdf(0.25):.3f}")     # ~0.88
print(f"P(rate > 0.35 | data) = {1 - post.cdf(0.35):.3f}")       # ~0.15
# These are direct probability statements about the parameter, which a
# frequentist p-value cannot provide.

Key Takeaways

• The posterior is proportional to the likelihood times the prior, combining prior belief with observed evidence.
• Priors are explicit and auditable assumptions rather than hidden ones.
• With sufficient data the likelihood dominates and different reasonable priors converge to the same posterior.
• The posterior supports direct probability statements about parameters, which frequentist methods cannot provide.

🧪 Practice

1. Update a Beta(1,1) prior with 8 successes in 20 trials and report the posterior mean and interval.
2. Show how three different priors converge as the sample grows from 10 to 10,000 observations.
3. Interview: A colleague objects that priors are subjective. How do you respond? (Hint: consider what assumptions a frequentist analysis makes and whether they are stated anywhere.)

Conjugate Priors

Bayes' theorem requires dividing by the probability of the data, which is an integral over all parameter values. For most models that integral has no closed form, which is why Bayesian methods were impractical for most of their history. Conjugate priors are the special cases where the mathematics works out exactly.

A prior is conjugate to a likelihood when the resulting posterior belongs to the same distribution family as the prior. Updating then reduces to arithmetic on the parameters instead of integration. Beta is conjugate to the binomial, so a Beta prior updated with binomial data yields a Beta posterior whose parameters are simply the prior's plus the observed counts.

This produces an unusually clear interpretation: the prior parameters act as pseudo-counts, observations you are pretending to have seen already. A Beta(8, 2) prior behaves exactly like having previously observed 8 successes and 2 failures. That makes the prior's strength concrete and easy to justify or challenge.

Conjugate updates are also what make streaming Bayesian inference trivial. Since the posterior has the same form as the prior, it becomes the prior for the next batch, and the update is a few additions.

Likelihood	Conjugate prior	Posterior parameters
Binomial	Beta(a, b)	Beta(a + successes, b + failures)
Poisson	Gamma(a, b)	Gamma(a + sum of counts, b + n)
Normal, known var	Normal	Precision-weighted blend of both means
Multinomial	Dirichlet	Dirichlet(alpha + category counts)
Exponential	Gamma	Gamma(a + n, b + sum of observations)

  PRIOR PARAMETERS AS PSEUDO-COUNTS

   Beta(a, b)  behaves exactly like having already seen
               a - 1 successes and b - 1 failures

   Beta(1, 1)    = 0 pseudo-observations   -> uniform, no information
   Beta(3, 7)    = weak belief near 0.3    -> 8 pseudo-observations
   Beta(30, 70)  = strong belief near 0.3  -> 98 pseudo-observations

   The SAME mean, very different strength:

     Beta(3,7)                    Beta(30,70)
        _---_                          /\
      .'     '.                       /  \
     /         \                     /    \
    ------------                   ---------
    wide: easily moved             narrow: needs a lot of data to shift

  STREAMING UPDATES ARE FREE

   day 1: Beta(1,1)   + (12 conv, 88 no)  -> Beta(13, 89)
   day 2: Beta(13,89) + (15 conv, 85 no)  -> Beta(28, 174)
   day 3: Beta(28,174)+ (11 conv, 89 no)  -> Beta(39, 263)

   No refitting, no stored history. Just addition.

import numpy as np
from scipy import stats

# Pseudo-count interpretation, demonstrated.
for a, b in [(1, 1), (3, 7), (30, 70)]:
    prior = stats.beta(a, b)
    # Observe a run of data that disagrees with the prior: 70% success.
    post = stats.beta(a + 70, b + 30)
    print(f"Beta({a:>2},{b:>2}): prior mean={prior.mean():.3f} "
          f"sd={prior.std():.3f} -> posterior mean={post.mean():.3f}")
# Beta( 1, 1): prior mean=0.500 sd=0.289 -> posterior mean=0.696
# Beta( 3, 7): prior mean=0.300 sd=0.138 -> posterior mean=0.664
# Beta(30,70): prior mean=0.300 sd=0.046 -> posterior mean=0.500
# The strong prior resists the data: 100 observations at 70% only move
# it halfway, because it already carries 98 pseudo-observations.

# Streaming: the posterior becomes the next prior, at no cost.
a, b = 1.0, 1.0
daily = [(12, 88), (15, 85), (11, 89), (18, 82)]
for day, (conv, no_conv) in enumerate(daily, start=1):
    a, b = a + conv, b + no_conv
    post = stats.beta(a, b)
    lo, hi = post.interval(0.95)
    print(f"day {day}: Beta({a:.0f}, {b:.0f}) mean={post.mean():.4f} "
          f"95% CI=({lo:.4f}, {hi:.4f})")
# The interval tightens each day and no past data needs to be retained.

# Poisson-Gamma: estimating a support ticket arrival rate.
observed_counts = np.array([3, 5, 2, 7, 4, 6])     # tickets per hour
prior_a, prior_b = 2.0, 1.0                        # weak belief near 2/hour
post_a = prior_a + observed_counts.sum()
post_b = prior_b + len(observed_counts)
rate_post = stats.gamma(a=post_a, scale=1 / post_b)
print(f"\nrate posterior mean = {rate_post.mean():.3f} per hour")
print(f"95% credible interval = "
      f"({rate_post.ppf(0.025):.2f}, {rate_post.ppf(0.975):.2f})")

# Verify conjugacy numerically against a brute-force grid posterior.
grid = np.linspace(0.01, 15, 20_000)
log_post = (stats.gamma.logpdf(grid, a=prior_a, scale=1 / prior_b)
            + np.array([stats.poisson.logpmf(observed_counts, g).sum()
                        for g in grid]))
dens = np.exp(log_post - log_post.max())
dens /= np.trapezoid(dens, grid)
grid_mean = np.trapezoid(grid * dens, grid)
print(f"grid-computed mean  = {grid_mean:.3f}")    # matches the closed form

Key Takeaways

• A conjugate prior makes the posterior fall in the same family, reducing Bayesian updating to arithmetic.
• Prior parameters act as pseudo-counts, which makes the prior's strength explicit and easy to justify.
• A strong prior resists new data in proportion to the pseudo-observations it encodes.
• Conjugacy makes streaming updates free, since the posterior becomes the next prior with no stored history.

🧪 Practice

1. Show that a Beta(2,2) prior with 10 successes in 15 trials gives a Beta(12,7) posterior, and report its mean.
2. Implement a streaming Beta-binomial update over a week of daily data.
3. Interview: How would you choose between Beta(1,1) and Beta(30,70) as a prior for a new feature's conversion rate? (Hint: ask how many comparable features you have already launched.)

Bayesian Inference vs Frequentist Inference

The two schools answer different questions, and most confusion comes from expecting one to answer the other's. Neither is universally correct; they make different trade-offs about what probability means and what you are allowed to say.

The frequentist view treats probability as long-run frequency. A parameter is a fixed unknown constant, not a random variable, so it cannot have a distribution. Randomness lives in the data, and therefore in the procedure's behaviour across hypothetical repetitions. That is why confidence attaches to the interval-producing method rather than to any particular interval.

The Bayesian view treats probability as a degree of belief. A parameter can have a distribution because your knowledge about it is uncertain, even though its true value is fixed. That permits direct statements like "there is a 94% probability the rate exceeds 10%".

The practical differences matter most in three places. Bayesian methods incorporate prior information, which is valuable with small samples and awkward when priors are contested. Bayesian methods permit continuous monitoring, since the posterior is valid whenever you look, whereas repeated frequentist testing inflates the false positive rate. Frequentist methods have guarantees that do not depend on a prior being reasonable, and they are computationally cheaper.

Aspect	Frequentist	Bayesian
Probability means	Long-run frequency	Degree of belief
Parameter is	Fixed unknown constant	Random variable with a distribution
Prior information	Not used	Required and explicit
Interval	Confidence: covers 95% of the time	Credible: 95% probability inside
Answers	P(data given H)	P(H given data)
Peeking at data	Inflates the error rate	Valid at any time
Computation	Usually cheap	Often expensive
Small samples	Weak	Strong, if the prior is good

  THE SAME EXPERIMENT, TWO ANSWERS

   data: 30 conversions in 100 trials

   FREQUENTIST (Wilson interval)
     point estimate   0.300
     95% CI           (0.219, 0.396)
     statement        "if I repeated this study many times, 95% of
                       intervals built this way would cover the truth"
     cannot say       "there is a 95% chance the rate is in here"

   BAYESIAN (uniform prior)
     posterior        Beta(31, 71)
     posterior mean   0.304
     95% credible     (0.221, 0.394)
     statement        "given the data and my prior, there is a 95%
                       probability the rate lies in here"
     can also say     "P(rate > 0.25) = 0.87"

   The numbers are nearly identical. The INTERPRETATIONS differ
   completely, and only one of them matches how people speak.

import numpy as np
from scipy import stats

successes, trials = 30, 100

# Frequentist: Wilson interval, which behaves better than the normal
# approximation for proportions near 0 or 1.
from statsmodels.stats.proportion import proportion_confint
lo_f, hi_f = proportion_confint(successes, trials, alpha=0.05, method="wilson")
print(f"frequentist 95% CI       = ({lo_f:.3f}, {hi_f:.3f})")

# Bayesian with a uniform prior.
post = stats.beta(1 + successes, 1 + trials - successes)
lo_b, hi_b = post.interval(0.95)
print(f"Bayesian 95% credible    = ({lo_b:.3f}, {hi_b:.3f})")
# Very close numerically, but they mean different things.

# The peeking problem: where the two genuinely diverge in practice.
rng = np.random.default_rng(0)

def frequentist_peeking(trials=2000, checks=20):
    """Check significance repeatedly under a TRUE null."""
    false_positives = 0
    for _ in range(trials):
        a, b = [], []
        for _ in range(checks):
            a.extend(rng.binomial(1, 0.12, 100))   # identical rates
            b.extend(rng.binomial(1, 0.12, 100))
            table = [[sum(a), len(a) - sum(a)], [sum(b), len(b) - sum(b)]]
            if stats.chi2_contingency(table).pvalue < 0.05:
                false_positives += 1
                break                               # stop at first "win"
    return false_positives / trials

print(f"\nfrequentist false positive rate with 20 peeks: "
      f"{frequentist_peeking():.1%}")
# ~25-30%, far above the nominal 5%. Each look is another chance to
# cross the threshold by luck, and stopping when you win locks it in.
#
# A Bayesian posterior does not have this problem: it is a summary of
# the evidence you currently hold, and looking at it does not change
# what it says. But note the honest caveat -- if you STOP collecting
# data when the posterior looks good, your decision rule still has
# frequentist properties worth checking.

Key Takeaways

• Frequentists treat the parameter as fixed and the data as random; Bayesians treat the parameter as uncertain and the data as given.
• A confidence interval describes the procedure; a credible interval describes the parameter.
• Repeated frequentist testing inflates the false positive rate sharply, while a posterior is valid whenever you inspect it.
• Neither framework is universally correct; they trade prior information against assumption-free guarantees.

🧪 Practice

1. Compute a frequentist confidence interval and a Bayesian credible interval for the same proportion and state each interpretation precisely.
2. Simulate the peeking problem and measure how the false positive rate grows with the number of checks.
3. Interview: When would you strongly prefer a Bayesian approach? (Hint: consider a launch with 40 users and ten comparable past launches.)

Credible Intervals

A credible interval is the Bayesian answer to "where is the parameter likely to be?" and it says exactly what most people incorrectly believe a confidence interval says.

A 95% credible interval is a range containing 95% of the posterior probability mass. Given the data and the prior, there is a 95% probability the parameter lies inside it. No reference to hypothetical repeated experiments is needed.

There are two common constructions. The equal-tailed interval cuts 2.5% from each end of the posterior, which is simple and the usual default. The highest density interval is the shortest interval containing 95% of the mass, which guarantees every point inside is more probable than every point outside. For a symmetric posterior they coincide. For a skewed one they differ, and the highest density interval is preferable because the equal-tailed version can include less probable values while excluding more probable ones.

The practical advantage over confidence intervals is that credible intervals compose naturally. You can ask directly for the probability that a parameter exceeds a threshold, or that one variant beats another, by integrating the posterior. These are the quantities decisions actually need.

  EQUAL-TAILED vs HIGHEST DENSITY on a SKEWED posterior

        posterior density
         |     _
         |    / \
         |   /   \__
         |  /       \___
         | /            \______
         +----------------------------
           |   |              |     |
          HDI_lo             HDI_hi |
             |                      |
            ET_lo                 ET_hi

   EQUAL-TAILED: 2.5% cut from each side. Simple, but on a skewed
                 posterior it excludes high-density values on the left
                 while including low-density values on the right.

   HIGHEST DENSITY: the SHORTEST interval holding 95%. Every point
                 inside is more probable than every point outside.

  WHAT A CREDIBLE INTERVAL LETS YOU SAY

   "P(conversion rate > 0.25) = 0.87"
   "P(variant B beats variant A) = 0.96"
   "P(lift is between 1% and 3%) = 0.62"

   All are integrals of the posterior. None has a frequentist analogue
   that means the same thing.

import numpy as np
from scipy import stats

# A skewed posterior, where the two interval types visibly differ.
post = stats.beta(2, 20)          # low rate, right-skewed

# Equal-tailed: 2.5% from each end.
et_lo, et_hi = post.ppf(0.025), post.ppf(0.975)
print(f"equal-tailed  95%: ({et_lo:.4f}, {et_hi:.4f})  "
      f"width={et_hi - et_lo:.4f}")

def highest_density_interval(dist, mass=0.95, grid_size=200_000):
    """Shortest interval containing `mass` of the probability."""
    lows = np.linspace(0, 1 - mass, grid_size)          # candidate left tails
    los = dist.ppf(lows)
    his = dist.ppf(lows + mass)
    widths = his - los
    k = int(np.argmin(widths))                          # the narrowest one
    return los[k], his[k]

hdi_lo, hdi_hi = highest_density_interval(post)
print(f"highest density 95%: ({hdi_lo:.4f}, {hdi_hi:.4f})  "
      f"width={hdi_hi - hdi_lo:.4f}")
# The HDI is narrower and shifted left, because the posterior's mass
# piles up near zero. The equal-tailed version wastes width on a long
# thin right tail while cutting off denser values on the left.

# Direct decision probabilities, which is the real payoff.
posterior_a = stats.beta(1 + 120, 1 + 880)      # variant A: 120/1000
posterior_b = stats.beta(1 + 150, 1 + 850)      # variant B: 150/1000

rng = np.random.default_rng(0)
samples_a = posterior_a.rvs(200_000, random_state=rng)
samples_b = posterior_b.rvs(200_000, random_state=rng)

print(f"\nP(B > A)              = {(samples_b > samples_a).mean():.4f}")
print(f"P(B beats A by >2pp)  = {((samples_b - samples_a) > 0.02).mean():.4f}")
lift = samples_b - samples_a
print(f"expected lift         = {lift.mean():.4f}")
print(f"95% credible on lift  = "
      f"({np.percentile(lift, 2.5):.4f}, {np.percentile(lift, 97.5):.4f})")
# P(B > A) is near 0.98, and the expected lift is about 3pp. The
# business can act on these numbers directly. A p-value would only
# have said "the difference is unlikely under the null".

# Expected loss from choosing wrongly -- the quantity that should
# actually drive a ship / no-ship decision.
loss_if_ship_b = np.maximum(samples_a - samples_b, 0).mean()
loss_if_keep_a = np.maximum(samples_b - samples_a, 0).mean()
print(f"\nexpected loss shipping B = {loss_if_ship_b:.5f}")
print(f"expected loss keeping  A = {loss_if_keep_a:.5f}")
# Ship whichever has the smaller expected loss. This is a decision rule,
# not a significance threshold.

Key Takeaways

• A credible interval contains a stated probability of the parameter, which is what people wrongly assume a confidence interval means.
• Equal-tailed intervals are the simple default; highest density intervals are shorter and preferable on skewed posteriors.
• Posteriors compose, so you can compute the probability one variant beats another directly.
• Expected loss derived from the posterior is a better ship decision rule than a significance threshold.

🧪 Practice

1. Compute both interval types for a Beta(2,20) posterior and explain why they differ.
2. Given two Beta posteriors, estimate the probability that one exceeds the other by sampling.
3. Interview: Your posterior says P(B beats A) = 0.91. Do you ship? (Hint: ask what shipping costs if you are in the other 9%.)

Markov Chain Monte Carlo

Conjugate priors cover a small set of tidy models. For anything realistic, with several parameters, hierarchical structure, or a non-standard likelihood, the posterior has no closed form and the normalizing integral is intractable. Markov chain Monte Carlo is the method that made Bayesian inference practical for real models.

The key realization is that you rarely need the posterior's formula. You need to compute means, intervals, and probabilities from it, and all of those can be estimated from samples. MCMC generates samples from the posterior without ever computing the normalizing constant.

It works by constructing a random walk through parameter space whose long-run visiting frequency matches the posterior. Steps toward higher-probability regions are always accepted; steps toward lower probability are accepted with a probability equal to the density ratio. That ratio cancels the unknown normalizing constant, which is the entire trick.

Because samples arrive as a correlated chain, diagnostics are mandatory. You must discard an initial burn-in period before the chain reaches the posterior, run multiple chains from different starting points, and check convergence. The R-hat statistic compares within-chain to between-chain variance and should sit below about 1.01. Effective sample size reports how many independent samples your correlated draws are worth.

  THE CHAIN EXPLORES THE POSTERIOR

   parameter
      ^
      |   burn-in        stationary sampling
      |  /\  /                                     the chain forgets
      | /  \/  ,-,_,-'-,_,'-,_,-'-,_,-'-,_,'       its starting point
      |/        (samples from here are usable)     and then wanders
      +-----|----------------------------------->  the posterior in
      0   discard              iteration           proportion to density

  METROPOLIS-HASTINGS ACCEPTANCE

   propose theta'  from  theta
   ratio = p(theta' | data) / p(theta | data)
         = [ L(theta') prior(theta') ] / [ L(theta) prior(theta) ]
                                           the normalizing constant
                                           CANCELS -- this is the trick
   if ratio >= 1  -> always accept  (uphill)
   else           -> accept with probability = ratio  (sometimes downhill,
                                                       which is what makes
                                                       it explore, not just
                                                       climb)

  DIAGNOSTICS THAT ARE NOT OPTIONAL
   R-hat < 1.01          chains agree with each other
   ESS > ~400 per chain  enough effective independent draws
   trace plot            should look like a fuzzy caterpillar, not a trend

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# Data: normal with unknown mean and unknown standard deviation.
true_mu, true_sigma = 5.0, 2.0
data = rng.normal(true_mu, true_sigma, size=200)

def log_posterior(mu, log_sigma):
    """Unnormalized log posterior. The constant is never needed."""
    sigma = np.exp(log_sigma)                  # sample log sigma to keep it >0
    log_prior = (stats.norm.logpdf(mu, 0, 10)          # weak prior on mu
                 + stats.norm.logpdf(log_sigma, 0, 2)) # weak prior on log sd
    log_like = stats.norm.logpdf(data, mu, sigma).sum()
    return log_prior + log_like

def metropolis(n_steps=40_000, step=0.15, seed=0):
    r = np.random.default_rng(seed)
    mu, log_sigma = 0.0, 0.0                   # deliberately poor start
    current = log_posterior(mu, log_sigma)
    chain, accepted = np.empty((n_steps, 2)), 0
    for i in range(n_steps):
        prop_mu = mu + r.normal(0, step)
        prop_ls = log_sigma + r.normal(0, step)
        proposed = log_posterior(prop_mu, prop_ls)
        # Compare in LOG space: log(ratio) = proposed - current.
        if np.log(r.random()) < proposed - current:
            mu, log_sigma, current = prop_mu, prop_ls, proposed
            accepted += 1
        chain[i] = (mu, np.exp(log_sigma))
    return chain, accepted / n_steps

chain, acceptance = metropolis()
print(f"acceptance rate = {acceptance:.3f}")    # aim for roughly 0.2 to 0.5

burn = 10_000                                   # discard the approach phase
samples = chain[burn:]
print(f"mu:    mean={samples[:,0].mean():.3f}  "
      f"95% CI=({np.percentile(samples[:,0], 2.5):.3f}, "
      f"{np.percentile(samples[:,0], 97.5):.3f})   true={true_mu}")
print(f"sigma: mean={samples[:,1].mean():.3f}  "
      f"95% CI=({np.percentile(samples[:,1], 2.5):.3f}, "
      f"{np.percentile(samples[:,1], 97.5):.3f})   true={true_sigma}")

# R-hat: run several chains from different seeds and compare variances.
chains = [metropolis(seed=s)[0][burn:, 0] for s in range(4)]
chains = np.array(chains)
m, n = chains.shape
W = chains.var(axis=1, ddof=1).mean()                    # within-chain
B = n * chains.mean(axis=1).var(ddof=1)                  # between-chain
var_est = ((n - 1) / n) * W + B / n
r_hat = np.sqrt(var_est / W)
print(f"\nR-hat for mu = {r_hat:.4f}")          # want < 1.01
# If R-hat is well above 1.01, the chains have not mixed and every
# posterior summary computed from them is untrustworthy.

# In practice use a probabilistic programming library (PyMC, Stan,
# NumPyro). They implement Hamiltonian Monte Carlo, which explores far
# more efficiently, and compute these diagnostics automatically.

Key Takeaways

• MCMC estimates posterior quantities from samples, avoiding the intractable normalizing integral.
• The Metropolis acceptance ratio cancels the unknown constant, which is what makes the method possible.
• Samples are correlated, so burn-in must be discarded and convergence must be checked.
• Check R-hat below about 1.01 and adequate effective sample size before trusting any posterior summary.

🧪 Practice

1. Implement Metropolis-Hastings for a single-parameter model and compare the result to the conjugate answer.
2. Run four chains from different starting points and compute R-hat.
3. Interview: Your MCMC trace plot shows a slow upward drift. What is wrong? (Hint: ask whether the chain has reached its stationary distribution yet.)

Probabilistic Programming Basics

Writing your own sampler, as above, is instructive once and impractical thereafter. Probabilistic programming languages let you declare a model as a generative story and leave inference to a well-tested engine.

The shift is from how to what. You specify which parameters exist, what priors they have, and how the data is generated from them. The library handles gradients, sampling, adaptation, and diagnostics. Modern engines use Hamiltonian Monte Carlo and the No-U-Turn Sampler, which use the posterior's gradient to propose distant, high-quality moves rather than random small steps, exploring far more efficiently than plain Metropolis.

The capability that makes this genuinely valuable is hierarchical models. When you have many related groups, such as conversion rates across 50 landing pages, you can model each group's parameter as drawn from a shared population distribution. Groups with little data are then pulled toward the overall mean, an effect called partial pooling or shrinkage. This is automatic regularization with a principled justification, and it outperforms both analysing each group separately and lumping them all together.

  THREE WAYS TO HANDLE MANY GROUPS

   NO POOLING            COMPLETE POOLING       PARTIAL POOLING
   each group alone      one rate for all       hierarchical

   p1  p2  p3 ... p50         p                    mu, tau
    |   |   |      |          |                  /  |  |  \
   y1  y2  y3 ... y50      y1..y50             p1  p2 p3 .. p50
                                                |   |  |      |
                                               y1  y2 y3 .. y50

   overfits small         ignores real          small groups borrow
   groups; a page         differences           strength from the
   with 3 visits and                            population; large
   1 conversion gets                            groups keep their own
   an estimate of 33%                           estimate

  SHRINKAGE IN ACTION

   raw rate    n      partially pooled estimate
     0.33      3   ->   0.14    pulled hard toward the mean
     0.30     50   ->   0.22    pulled moderately
     0.28    900   ->   0.277   barely moved; it has its own evidence

   The amount of shrinkage is learned from the data, not chosen.

# Sketch in PyMC. The model reads as the generative story itself.
import pymc as pm
import numpy as np

# 50 landing pages with very unequal traffic.
rng = np.random.default_rng(0)
n_pages = 50
visits = rng.integers(3, 1000, n_pages)
true_rates = rng.beta(2, 18, n_pages)               # true rate near 0.10
conversions = rng.binomial(visits, true_rates)

with pm.Model() as hierarchical:
    # POPULATION level: what rates look like across pages in general.
    mu = pm.Beta("mu", alpha=2, beta=20)            # typical rate
    kappa = pm.HalfNormal("kappa", sigma=50)        # how tightly grouped

    # GROUP level: each page's rate is drawn from that population.
    rates = pm.Beta("rates",
                    alpha=mu * kappa,
                    beta=(1 - mu) * kappa,
                    shape=n_pages)

    # OBSERVED data, linking parameters to the counts we actually saw.
    pm.Binomial("obs", n=visits, p=rates, observed=conversions)

    # NUTS runs automatically; diagnostics are computed for you.
    trace = pm.sample(2000, tune=1000, chains=4, target_accept=0.9,
                      random_seed=0, progressbar=False)

import arviz as az
summary = az.summary(trace, var_names=["mu", "kappa"])
print(summary[["mean", "hdi_3%", "hdi_97%", "r_hat", "ess_bulk"]])
# Always read r_hat (want < 1.01) and ess_bulk (want > 400) before
# interpreting anything else in the table.

# Shrinkage, made visible.
posterior_rates = trace.posterior["rates"].mean(dim=["chain", "draw"]).values
raw_rates = conversions / visits
for i in np.argsort(visits)[[0, 1, -2, -1]]:       # two smallest, two largest
    print(f"visits={visits[i]:>4}  raw={raw_rates[i]:.3f}  "
          f"pooled={posterior_rates[i]:.3f}  true={true_rates[i]:.3f}")
# Low-traffic pages are pulled sharply toward the population mean and
# end up CLOSER to their true rate than the raw estimate was.
# High-traffic pages barely move, because their own data dominates.

Library	Language	Sampler	Notable for
PyMC	Python	NUTS	Readable models, strong ecosystem
Stan	Own DSL	NUTS	Reference implementation, fast
NumPyro	Python	NUTS on JAX	GPU acceleration, very fast
Turing	Julia	Several	Composable with Julia code

Key Takeaways

• Probabilistic programming lets you declare the generative model and delegates inference to a tested engine.
• Hamiltonian Monte Carlo and NUTS use posterior gradients and explore far more efficiently than plain Metropolis.
• Hierarchical models let groups borrow strength, shrinking small-sample estimates toward the population mean.
• The degree of shrinkage is learned from the data rather than chosen, which makes it principled regularization.

🧪 Practice

1. Fit a hierarchical Beta-binomial model over several groups and compare pooled estimates to raw rates.
2. Show that groups with the fewest observations are shrunk the most.
3. Interview: Why does partial pooling usually beat fitting each group separately? (Hint: consider a group with three observations and what its unpooled estimate looks like.)

<a id="35-causal-inference"></a>

3.5 Causal Inference

Stakeholders never ask what is correlated with revenue; they ask what will happen to revenue if they act, and this section covers the assumptions and designs that let observational data answer that question.

Correlation vs Causation

Every statistics course repeats that correlation does not imply causation, and every analyst nods and then proceeds to build a model whose coefficients they describe causally. The slogan is easy; the discipline is not.

The reason the distinction matters is entirely practical. Correlation supports prediction: knowing X helps you guess Y. Causation supports intervention: changing X will change Y. A model can be excellent at the first and worthless for the second, and it will not tell you which.

Four structures produce correlation without a direct causal link. Confounding: a third variable causes both. Reverse causation: Y causes X rather than the reverse. Selection effects: the sampling process induces an association that does not exist in the population. Coincidence: with enough variables tested, some correlate by chance.

The canonical example is that hospitals correlate with death. Modelling mortality from hospital visits yields a strong positive coefficient. As prediction it is genuinely useful. As a guide to action, "avoid hospitals" is lethal, because sickness causes both the visit and the risk.

An analogy: a thermometer predicts fever accurately. Putting it in the freezer does not cure anyone. Prediction reads the world; intervention changes it, and the relationship that held while you were only watching need not survive your interference.

  FOUR WAYS TO GET CORRELATION WITHOUT DIRECT CAUSATION

  CONFOUNDING              REVERSE CAUSATION
       Z                        X <--- Y
      / \                   observed as X ~ Y, but the arrow
     v   v                  points the other way
     X   Y                  e.g. "ads cause sales" when in fact
   X ~ Y with no arrow      budgets are set FROM last month's sales
   between them
   e.g. sickness -> both
        hospital and death

  SELECTION / COLLIDER              COINCIDENCE
     X ---> C <--- Y                  test 1000 pairs at alpha=0.05
           |                          and about 50 will "correlate"
      condition on C                  with no relationship at all
   Conditioning on a common
   EFFECT creates a spurious
   link between its causes.
   e.g. among admitted students,
   test scores and grades appear
   negatively related

  THE TEST THAT MATTERS
   "If I intervene on X and change nothing else, does Y move?"
   Only a randomized experiment answers this directly. Everything
   else in this section is a way to approximate it.

import numpy as np
import pandas as pd
import statsmodels.api as sm

rng = np.random.default_rng(0)
n = 5000

# CONFOUNDING, built deliberately so we know the truth.
# Severity causes BOTH hospital admission and death. Hospital care
# actually REDUCES death, but the naive analysis will say otherwise.
severity = rng.normal(0, 1, n)
hospital = (rng.random(n) < 1 / (1 + np.exp(-2 * severity))).astype(int)
death_risk = -3 + 2.0 * severity - 0.8 * hospital     # hospital HELPS
death = (rng.random(n) < 1 / (1 + np.exp(-death_risk))).astype(int)

df = pd.DataFrame({"severity": severity, "hospital": hospital, "death": death})

naive = sm.Logit(df["death"], sm.add_constant(df[["hospital"]])).fit(disp=0)
print(f"naive hospital coefficient     = {naive.params['hospital']:+.3f}")
# POSITIVE: "hospitals kill". The true effect is -0.8.

adjusted = sm.Logit(df["death"],
                    sm.add_constant(df[["hospital", "severity"]])).fit(disp=0)
print(f"adjusted hospital coefficient  = {adjusted.params['hospital']:+.3f}")
print(f"true coefficient               = -0.800")
# Controlling for the confounder recovers roughly the truth. This works
# ONLY because we know severity is the confounder and we measured it.

# COLLIDER BIAS: adjusting for the WRONG variable creates bias.
# Talent and luck are independent causes of getting hired.
talent = rng.normal(0, 1, n)
luck = rng.normal(0, 1, n)
hired = (talent + luck) > 1.5                    # a common EFFECT

print(f"\ncorr(talent, luck) overall        = "
      f"{np.corrcoef(talent, luck)[0,1]:+.3f}")            # ~0.00
print(f"corr(talent, luck) among the hired = "
      f"{np.corrcoef(talent[hired], luck[hired])[0,1]:+.3f}")  # ~-0.60
# Conditioning on a collider MANUFACTURED a strong negative correlation
# between two independent variables. "Controlling for more variables"
# is not automatically safer; it can create bias where none existed.

Key Takeaways

• Correlation supports prediction; only causation supports intervention, and a model cannot tell you which it has captured.
• Confounding, reverse causation, selection effects, and coincidence all produce correlation without a direct causal link.
• Adjusting for a confounder removes bias, but only if you know about it and measured it.
• Adjusting for a collider, a common effect of two variables, creates bias rather than removing it.

🧪 Practice

1. Simulate a confounded relationship and show how the naive and adjusted coefficients differ from the truth.
2. Construct a collider and demonstrate that conditioning on it induces a correlation between independent variables.
3. Interview: A model shows customers who contact support churn more. Should you reduce support availability? (Hint: ask what makes a person contact support in the first place.)

Confounding Variables

A confounder is a variable that causes both the treatment and the outcome. Its presence means the comparison you make is not the comparison you think you are making, because the treated and untreated groups differ in ways beyond the treatment.

The formal criterion is worth stating precisely: a variable confounds if it causes the treatment and it causes the outcome, and it is not itself caused by the treatment. That last clause matters enormously. A variable on the causal path from treatment to outcome is a mediator, and controlling for it removes part of the very effect you are measuring.

This gives three categories with three different correct actions. Control for confounders. Do not control for mediators, unless you specifically want the direct effect only. Never control for colliders.

The uncomfortable truth is that you cannot verify from data alone which category a variable falls into. All three can produce identical correlation patterns. The classification comes from domain knowledge about what causes what, which is why causal inference is not a purely statistical exercise.

Omitted variable bias quantifies the damage from missing a confounder. Its size is the product of two things: how strongly the confounder relates to the treatment, and how strongly it relates to the outcome. Both must be substantial for the bias to be large, which is the basis of sensitivity analysis.

  THREE ROLES, THREE DIFFERENT ACTIONS

  CONFOUNDER            MEDIATOR              COLLIDER
      Z                                            
     / \              T -> M -> Y            T -> C <- Y
    v   v                                          
    T   Y                                          
                                                    
  CONTROL FOR IT       do NOT control        NEVER control
  removes bias         (it is part of        (creates bias)
                       the effect)

  EXAMPLE, ALL THREE AT ONCE
   T = exercise programme,  Y = heart disease

   age            causes both exercise and disease  -> CONFOUNDER, control
   blood pressure exercise lowers it, it lowers     -> MEDIATOR, do not
                  disease
   gym membership caused by exercise AND by health  -> COLLIDER, never
                  consciousness which affects Y

  OMITTED VARIABLE BIAS

   bias  =  (effect of Z on T)  x  (effect of Z on Y)

   Both factors must be large for the bias to matter. A confounder
   strongly linked to treatment but weakly to outcome does little
   damage, and vice versa.

import numpy as np
import pandas as pd
import statsmodels.api as sm

rng = np.random.default_rng(0)
n = 10_000

# Build a world with a confounder, a mediator, and a collider.
age = rng.normal(50, 12, n)                                   # CONFOUNDER
exercise = 0.5 * (age < 45) + rng.normal(0, 1, n)             # age -> exercise
blood_pressure = 130 - 3 * exercise + 0.4 * age + rng.normal(0, 5, n)  # MEDIATOR
TRUE_DIRECT = -1.5
disease = (TRUE_DIRECT * exercise                              # direct effect
           + 0.10 * blood_pressure                             # via mediator
           + 0.08 * age                                        # confounder
           + rng.normal(0, 1, n))
gym = 0.8 * exercise - 0.3 * disease + rng.normal(0, 1, n)     # COLLIDER

df = pd.DataFrame({"exercise": exercise, "age": age, "disease": disease,
                   "blood_pressure": blood_pressure, "gym": gym})

def coef(controls):
    X = sm.add_constant(df[["exercise"] + controls])
    return sm.OLS(df["disease"], X).fit().params["exercise"]

total_effect = TRUE_DIRECT + 0.10 * (-3)        # direct plus via mediator
print(f"true DIRECT effect          = {TRUE_DIRECT:+.3f}")
print(f"true TOTAL effect           = {total_effect:+.3f}")
print()
print(f"no controls                 = {coef([]):+.3f}   confounded by age")
print(f"+ age (confounder)          = {coef(['age']):+.3f}   recovers TOTAL")
print(f"+ age + blood_pressure      = "
      f"{coef(['age', 'blood_pressure']):+.3f}   gives DIRECT, not total")
print(f"+ age + gym (collider)      = "
      f"{coef(['age', 'gym']):+.3f}   BIASED by the collider")
# Controlling for age recovers the total effect of -1.8.
# Adding the mediator changes the question to the direct effect of -1.5.
# Adding the collider corrupts the estimate in a direction nothing in
# the data warns you about.

# Omitted variable bias is a product of two associations.
print("\nbias from omitting age depends on BOTH links:")
for age_to_exercise in (0.0, 0.5, 2.0):
    ex = age_to_exercise * (age < 45) + rng.normal(0, 1, n)
    dis = TRUE_DIRECT * ex + 0.08 * age + rng.normal(0, 1, n)
    naive = sm.OLS(dis, sm.add_constant(ex)).fit().params[1]
    print(f"  age->exercise strength {age_to_exercise}: "
          f"naive estimate = {naive:+.3f}")
# With no age-to-exercise link, age is not a confounder at all and the
# naive estimate is unbiased despite age still affecting the outcome.

Key Takeaways

• A confounder causes both treatment and outcome and is not caused by the treatment; control for it.
• A mediator lies on the causal path, so controlling for it estimates the direct rather than the total effect.
• You cannot tell confounders, mediators, and colliders apart from data alone; the classification requires domain knowledge.
• Omitted variable bias is the product of the confounder's link to the treatment and its link to the outcome.

🧪 Practice

1. Simulate a mediator and show how controlling for it changes the estimated effect from total to direct.
2. Vary the strength of a confounder's two links and plot how the naive bias responds.
3. Interview: Should you control for every available variable to be safe? (Hint: consider what happens when one of them is a collider.)

Potential Outcomes Framework

To define a causal effect precisely, you need language for something that did not happen. The potential outcomes framework, also called the Rubin causal model, supplies it, and it is the foundation on which the remaining methods in this section rest.

For each unit, imagine two potential outcomes. Y(1) is what would happen if treated; Y(0) is what would happen if not. The individual treatment effect is their difference. The problem is immediate and unavoidable: you only ever observe one of them. The other is the counterfactual, and it is missing, permanently. This is the fundamental problem of causal inference, and it means causal inference is a missing data problem.

Individual effects are therefore unidentifiable. What can be estimated are averages, most commonly the average treatment effect, the mean of Y(1) minus Y(0) over the population.

Randomization is what makes that average estimable. When treatment is assigned at random, the treated and untreated groups are statistically identical in every respect, observed and unobserved. The untreated group therefore serves as a valid stand-in for what would have happened to the treated group, and the difference in group means estimates the average effect without bias.

Two assumptions are required whenever randomization is absent. Ignorability says that, conditional on the measured covariates, treatment assignment is as good as random. Positivity says every unit has a non-zero chance of either treatment. Neither is testable from the data.

  THE FUNDAMENTAL PROBLEM

   unit   Y(1)    Y(0)    effect    observed
   ----   -----   -----   ------    --------------------
    1      120      ?       ?       treated,   saw Y(1)=120
    2       ?       95      ?       untreated, saw Y(0)=95
    3      140      ?       ?       treated,   saw Y(1)=140
    4       ?      110      ?       untreated, saw Y(0)=110

   Exactly half the table is unobservable, forever. Every causal
   method is a strategy for filling in the question marks credibly.

  WHY RANDOMIZATION WORKS

   random assignment  =>  treated and untreated groups are
                          identical in expectation on EVERYTHING,
                          measured or not

      treated group      untreated group
      [ mixed ages ]  ~  [ mixed ages ]
      [ mixed risk ]  ~  [ mixed risk ]
      [ unknown X  ]  ~  [ unknown X  ]   <- the crucial one

   So the untreated group's average outcome is a valid estimate of
   what the treated group WOULD have had without treatment.

  ESTIMANDS
   ATE   average treatment effect over everyone
   ATT   average effect on those actually treated
   CATE  conditional effect for a subgroup, the basis of personalization

import numpy as np
import pandas as pd

rng = np.random.default_rng(0)
n = 20_000

# A simulated world where BOTH potential outcomes are known, so the
# true answer is available for checking. Real data never allows this.
age = rng.normal(45, 12, n)
y0 = 50 + 0.5 * age + rng.normal(0, 5, n)              # outcome if untreated
effect = 8 - 0.1 * age                                  # effect varies by age
y1 = y0 + effect                                        # outcome if treated

TRUE_ATE = effect.mean()
print(f"TRUE ATE (unobservable in practice) = {TRUE_ATE:.3f}")

# SCENARIO 1: randomized assignment.
treated = rng.random(n) < 0.5
observed = np.where(treated, y1, y0)                    # only one is seen
ate_rct = observed[treated].mean() - observed[~treated].mean()
print(f"randomized estimate                 = {ate_rct:.3f}   unbiased")

# SCENARIO 2: self-selection. Younger people opt in more often, and
# younger people also benefit more. The naive difference is biased.
p_treat = 1 / (1 + np.exp((age - 45) / 8))
treated_bias = rng.random(n) < p_treat
observed_bias = np.where(treated_bias, y1, y0)
ate_naive = observed_bias[treated_bias].mean() - observed_bias[~treated_bias].mean()
print(f"self-selected naive estimate        = {ate_naive:.3f}   BIASED")

# Under ignorability given age, adjusting within age strata recovers it.
df = pd.DataFrame({"age": age, "treated": treated_bias, "y": observed_bias})
df["stratum"] = pd.qcut(df["age"], 20, labels=False)
strata = df.groupby("stratum").apply(
    lambda g: pd.Series({
        "diff": g.loc[g["treated"], "y"].mean() - g.loc[~g["treated"], "y"].mean(),
        "weight": len(g),
    }), include_groups=False)
ate_adjusted = np.average(strata["diff"], weights=strata["weight"])
print(f"stratified adjusted estimate        = {ate_adjusted:.3f}   recovered")

# ATT differs from ATE when the effect varies and selection is not random.
att = effect[treated_bias].mean()
print(f"\nTRUE ATE = {TRUE_ATE:.3f}   TRUE ATT = {att:.3f}")
# The treated are younger and benefit more, so ATT exceeds ATE. Which
# one answers your question depends on who you plan to treat next.

# POSITIVITY, violated: nobody over 70 ever gets treated.
never_treated = age > 70
print(f"\nunits with age > 70: {never_treated.sum()}, "
      f"of which treated: {(treated_bias & never_treated).sum()}")
# With no treated units in a region of covariate space, no amount of
# statistical adjustment can estimate the effect there. The honest
# move is to restrict the claim to the range where both groups exist.

Key Takeaways

• A causal effect is the difference between two potential outcomes, only one of which is ever observed.
• Causal inference is a missing data problem, and individual effects are never identifiable.
• Randomization balances the groups on everything, measured or not, which is why it estimates the average effect without bias.
• Ignorability and positivity are required without randomization, and neither can be tested from the data.

🧪 Practice

1. Simulate both potential outcomes, then compare a randomized estimate to a self-selected one against the known truth.
2. Construct a case where the average treatment effect and the effect on the treated differ substantially, and explain why.
3. Interview: Why can you never measure an individual's treatment effect? (Hint: write down what you would need to observe for one person.)

Causal Graphs and DAGs

Deciding what to control for requires making your causal assumptions explicit. A directed acyclic graph is the standard notation: variables are nodes, direct causal effects are arrows, and no path loops back on itself.

The value of a DAG is that it converts an argument about intuitions into a question with a mechanical answer. Once the graph is drawn, rules determine exactly which variables must be controlled to identify a given effect, and which must not be. The graph is an assumption, not a discovery, but it is a stated assumption others can challenge.

Three path structures are the building blocks. A chain, A -> B -> C, transmits association; controlling for B blocks it. A fork, A <- B -> C, also transmits association; controlling for B blocks it. A collider, A -> B <- C, blocks association by default; controlling for B opens it. That inversion is what makes colliders so dangerous.

The backdoor criterion uses these rules. To identify the effect of T on Y, find a set of variables that blocks every backdoor path, meaning every path from T to Y that starts with an arrow into T, while opening no collider paths. Control for exactly that set: no more, no less.

  THE THREE BUILDING BLOCKS

  CHAIN            FORK               COLLIDER
  A -> B -> C      A <- B -> C        A -> B <- C

  A and C are      A and C are        A and C are
  associated       associated         INDEPENDENT

  control B ->     control B ->       control B ->
  blocks it        blocks it          OPENS a spurious
                                      association

  A WORKED BACKDOOR EXAMPLE

              Age
             /   \
            v     v
      Exercise -> Disease        the effect we want: Exercise -> Disease
            \     ^
             v   /
          BloodPressure          (a mediator, on the causal path)

   paths from Exercise to Disease:
     1. Exercise -> Disease                    the direct causal path, KEEP
     2. Exercise -> BP -> Disease              also causal (indirect), KEEP
                                               for the TOTAL effect
     3. Exercise <- Age -> Disease             BACKDOOR (fork), must BLOCK

   backdoor set = {Age}.  Control for Age and nothing else.
   Adding BloodPressure would block path 2 and give the direct effect.
   Adding a collider would open a new spurious path.

# DAGs are declared, then queried. This uses the standard networkx
# representation plus explicit path logic, so the rules stay visible.
import networkx as nx
from itertools import combinations

g = nx.DiGraph()
g.add_edges_from([
    ("Age", "Exercise"), ("Age", "Disease"),
    ("Exercise", "BloodPressure"), ("BloodPressure", "Disease"),
    ("Exercise", "Disease"),
    ("Exercise", "Gym"), ("Disease", "Gym"),      # Gym is a COLLIDER
])
print("acyclic:", nx.is_directed_acyclic_graph(g))    # must be True

def backdoor_paths(graph, treatment, outcome):
    """Undirected paths from treatment to outcome starting with an
    arrow INTO the treatment. Those are the ones needing to be blocked."""
    undirected = graph.to_undirected()
    out = []
    for path in nx.all_simple_paths(undirected, treatment, outcome):
        # A backdoor path begins with an edge pointing INTO the treatment.
        if graph.has_edge(path[1], treatment):
            out.append(path)
    return out

for p in backdoor_paths(g, "Exercise", "Disease"):
    print("backdoor path:", " - ".join(p))
# backdoor path: Exercise - Age - Disease
# backdoor path: Exercise - Gym - Disease      <- via the COLLIDER

# The Gym path is already BLOCKED, because Gym is a collider on it and
# we are not conditioning on Gym. The Age path is open and must be
# blocked by controlling for Age.

# The textbook backdoor test: delete every edge LEAVING the treatment,
# then ask whether treatment and outcome are d-separated by the
# candidate set. If they are, that set blocks all backdoor paths.
def satisfies_backdoor(graph, treatment, outcome, adjust):
    stripped = graph.copy()
    stripped.remove_edges_from(list(graph.out_edges(treatment)))
    return nx.is_d_separator(stripped, {treatment}, {outcome}, set(adjust))

for candidate in [set(), {"Age"}, {"Age", "Gym"}]:
    ok = satisfies_backdoor(g, "Exercise", "Disease", candidate)
    print(f"adjust for {str(candidate or '{}'):<18} backdoor satisfied: {ok}")
# adjust for {}                 backdoor satisfied: False
# adjust for {'Age'}            backdoor satisfied: True
# adjust for {'Age', 'Gym'}     backdoor satisfied: False
#
# The last line is the important one. Adding the collider Gym does not
# merely waste a degree of freedom; it opens a path and breaks the
# criterion that {Age} alone satisfied.

# The practical payoff, stated plainly:
#   adjust for {Age}                 -> TOTAL effect of Exercise
#   adjust for {Age, BloodPressure}  -> DIRECT effect only
#   adjust for {Age, Gym}            -> BIASED, opens the collider path
#
# In production work, the `dowhy` library takes a DAG and a dataset and
# returns the valid adjustment sets automatically, then estimates the
# effect and runs refutation tests against it.

Key Takeaways

• A DAG makes causal assumptions explicit and turns "what should I control for?" into a question with a mechanical answer.
• Chains and forks transmit association and are blocked by conditioning; colliders block it and are opened by conditioning.
• The backdoor criterion says to block every path into the treatment while opening no collider paths.
• The graph is an assumption others can challenge, which is its main advantage over unstated intuition.

🧪 Practice

1. Draw a DAG for a marketing campaign's effect on purchases and list every backdoor path.
2. Identify the minimal adjustment set for your DAG and name one variable that must not be included.
3. Interview: Your DAG says not to control for a variable your stakeholder insists on including. How do you explain it? (Hint: show what happens to an estimate you both agree on when the variable is added.)

Propensity Score Matching

With observational data and a credible set of confounders, you could control for them in a regression. Matching takes a different approach that is often more transparent: construct a comparison group that resembles the treated group on the measured covariates, then compare outcomes directly.

Matching on many covariates at once quickly becomes impossible, because exact matches stop existing as dimensions grow. The propensity score solves this. It is the predicted probability of receiving treatment given the covariates, usually estimated by logistic regression, and it compresses all covariates into a single number. Rosenbaum and Rubin proved that matching on this one score balances all the covariates that went into it.

The workflow has four mandatory steps. Estimate the propensity score. Check overlap, meaning both groups span a common range of scores. Match, weight, or stratify on the score. Then check balance, confirming the covariates are actually similar after matching. Skipping the last step is the most common error; a matching procedure that failed to balance is worse than no matching, because it looks rigorous.

The fundamental limitation must be stated clearly. Propensity scores adjust only for measured covariates. Unlike randomization, they do nothing about unmeasured confounding. A propensity-matched study is an observational study with better bookkeeping, not an experiment.

  THE PROPENSITY SCORE COLLAPSES MANY COVARIATES INTO ONE

   age, income, tenure, region, plan, usage, ...
                     |
                     v   logistic regression on treatment
              e(X) = P(treated | X)
                     |
                     v
       match treated and control units with similar e(X)

  OVERLAP (common support) -- CHECK THIS FIRST

   good overlap                     poor overlap
    control  treated                 control      treated
     ___      ___                     ___             ___
    /   \    /   \                   /   \           /   \
   |     \  /     |                 |     |         |     |
   ------------------              -----------     -----------
   0        e(X)      1            0                          1
   most units can be               almost no comparable
   matched                         pairs exist; matching
                                   here is extrapolation

  BALANCE CHECK: standardized mean difference per covariate

              before matching   after matching
   age            0.62              0.04     good (< 0.10)
   income         0.48              0.07     good
   tenure         0.55              0.31     STILL IMBALANCED -> fix

  WHAT IT CANNOT DO
   unmeasured confounder U -> treatment and outcome
   No propensity model built from measured X touches this at all.

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import NearestNeighbors

rng = np.random.default_rng(0)
n = 4000

# Observational data: who enrols in the loyalty programme is not random.
age = rng.normal(45, 12, n)
tenure = rng.exponential(24, n)
prior_spend = rng.lognormal(4, 0.6, n)

# Enrolment depends on the covariates -> confounding.
logit_p = -3 + 0.03 * age + 0.02 * tenure + 0.004 * prior_spend
enrolled = rng.random(n) < 1 / (1 + np.exp(-logit_p))

TRUE_EFFECT = 25.0
spend = (100 + 1.2 * age + 0.8 * tenure + 0.5 * prior_spend
         + TRUE_EFFECT * enrolled + rng.normal(0, 20, n))

df = pd.DataFrame({"age": age, "tenure": tenure, "prior_spend": prior_spend,
                   "enrolled": enrolled.astype(int), "spend": spend})
covariates = ["age", "tenure", "prior_spend"]

naive = df.loc[df.enrolled == 1, "spend"].mean() - \
        df.loc[df.enrolled == 0, "spend"].mean()
print(f"true effect        = {TRUE_EFFECT:.2f}")
print(f"naive difference   = {naive:.2f}   badly overstated\n")

# STEP 1: estimate the propensity score.
ps_model = LogisticRegression(max_iter=1000).fit(df[covariates], df["enrolled"])
df["ps"] = ps_model.predict_proba(df[covariates])[:, 1]

# STEP 2: check overlap BEFORE matching anything.
t_ps, c_ps = df.loc[df.enrolled == 1, "ps"], df.loc[df.enrolled == 0, "ps"]
print(f"treated ps range {t_ps.min():.3f}-{t_ps.max():.3f}")
print(f"control ps range {c_ps.min():.3f}-{c_ps.max():.3f}")
lo, hi = max(t_ps.min(), c_ps.min()), min(t_ps.max(), c_ps.max())
df = df[(df.ps >= lo) & (df.ps <= hi)]          # restrict to common support

# STEP 3: nearest-neighbour matching, 1-to-1 without replacement.
treated = df[df.enrolled == 1].reset_index(drop=True)
control = df[df.enrolled == 0].reset_index(drop=True)
nn = NearestNeighbors(n_neighbors=1).fit(control[["ps"]])
distance, idx = nn.kneighbors(treated[["ps"]])

# A caliper discards poor matches rather than accepting them silently.
caliper = 0.2 * df["ps"].std()
keep = distance.ravel() <= caliper
matched_t = treated[keep]
matched_c = control.iloc[idx.ravel()[keep]]
print(f"\nmatched {keep.sum()} of {len(treated)} treated units")

# STEP 4: check BALANCE. Never skip this.
def smd(a, b):
    """Standardized mean difference: want |smd| < 0.10 after matching."""
    return (a.mean() - b.mean()) / np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)

print(f"\n{'covariate':<14}{'before':>9}{'after':>9}")
for c in covariates:
    before = smd(treated[c], control[c])
    after = smd(matched_t[c], matched_c[c])
    flag = "" if abs(after) < 0.10 else "  <- still imbalanced"
    print(f"{c:<14}{before:>9.3f}{after:>9.3f}{flag}")

matched_effect = matched_t["spend"].mean() - matched_c["spend"].mean()
print(f"\nmatched estimate   = {matched_effect:.2f}   close to {TRUE_EFFECT}")

# Inverse probability weighting is the common alternative: reweight so
# both groups resemble the full population rather than discarding units.
w = np.where(df.enrolled == 1, 1 / df.ps, 1 / (1 - df.ps))
w = np.clip(w, None, np.percentile(w, 99))          # trim extreme weights
ipw = sm.WLS(df["spend"], sm.add_constant(df[["enrolled"]]),
             weights=w).fit()
print(f"IPW estimate       = {ipw.params['enrolled']:.2f}")

Key Takeaways

• The propensity score compresses many covariates into one number, and matching on it balances all the covariates it was built from.
• Check overlap before matching and balance after; an unbalanced match is worse than none because it looks rigorous.
• Use a caliper so that poor matches are discarded rather than silently accepted.
• Propensity methods adjust only for measured covariates and provide no protection against unmeasured confounding.

🧪 Practice

1. Estimate propensity scores on an observational dataset and produce a before-and-after balance table.
2. Compare nearest-neighbour matching with inverse probability weighting on the same data.
3. Interview: Your matched analysis shows perfect balance on all measured covariates. How confident are you in the causal claim? (Hint: list what balance on measured variables does not tell you.)

Instrumental Variables

When an unmeasured confounder exists, adjustment methods fail, because you cannot control for something you did not measure. Instrumental variables offer an escape by finding a source of variation in the treatment that is unrelated to the confounder.

An instrument Z must satisfy three conditions. Relevance: it actually affects the treatment, and this one is testable. Exclusion: it affects the outcome only through the treatment, never directly. Independence: it is unrelated to the unmeasured confounders. The last two are assumptions and cannot be verified from the data, which is why a good instrument requires a persuasive substantive argument, not a statistical test.

The logic is that the instrument creates quasi-random variation in the treatment. The portion of treatment variation driven by the instrument is clean, and the effect is recovered by comparing how much the instrument moves the treatment against how much it moves the outcome.

Two-stage least squares implements this. The first stage predicts treatment from the instrument. The second stage regresses the outcome on those predictions, which contain only the clean variation.

A weak instrument, one only loosely related to the treatment, is worse than none: it amplifies any small violation of the other assumptions and produces badly biased estimates with misleadingly tight intervals. The conventional check is a first-stage F-statistic above 10.

  THE INSTRUMENT'S ROLE

        U (unmeasured confounder)
         / \
        v   v
   Z -> T -> Y
   ^
   |
   instrument: affects T, and reaches Y ONLY through T,
   and is unrelated to U

   FORBIDDEN:  Z -> Y directly      (exclusion violated)
               U -> Z               (independence violated)

  TWO-STAGE LEAST SQUARES

   stage 1:  T_hat = a + b*Z            keep only the part of T that
                                        the instrument explains
   stage 2:  Y     = c + d*T_hat        d estimates the causal effect

   The confounded part of T was discarded in stage 1.

  CLASSIC INSTRUMENTS
   distance to a college    -> years of education -> earnings
   draft lottery number     -> military service   -> earnings
   random encouragement     -> feature adoption   -> retention
   supply-side cost shocks  -> price              -> demand

  WEAK INSTRUMENT WARNING
   first-stage F < 10  ->  the estimate is unreliable and its
                           confidence interval is dishonestly narrow

import numpy as np
import pandas as pd
import statsmodels.api as sm
from linearmodels.iv import IV2SLS

rng = np.random.default_rng(0)
n = 5000

# An UNMEASURED confounder drives both treatment and outcome.
ability = rng.normal(0, 1, n)                     # never observed
distance_to_college = rng.normal(0, 1, n)         # the INSTRUMENT

# Education depends on the instrument AND on the confounder.
education = (12 + 1.5 * (-distance_to_college) + 2.0 * ability
             + rng.normal(0, 1, n))

TRUE_RETURN = 0.08
log_wage = (2.0 + TRUE_RETURN * education
            + 0.5 * ability                        # confounding path
            + rng.normal(0, 0.3, n))

df = pd.DataFrame({"education": education, "log_wage": log_wage,
                   "distance": distance_to_college, "ability": ability})

# NAIVE: biased upward, because ability raises both education and wages.
ols = sm.OLS(df["log_wage"], sm.add_constant(df[["education"]])).fit()
print(f"true return          = {TRUE_RETURN:.4f}")
print(f"naive OLS            = {ols.params['education']:.4f}   biased up")

# If we COULD measure ability, adjustment would work. We cannot, in
# reality; this line only exists to confirm the setup.
ctrl = sm.OLS(df["log_wage"],
              sm.add_constant(df[["education", "ability"]])).fit()
print(f"OLS + ability        = {ctrl.params['education']:.4f}   (infeasible)")

# RELEVANCE CHECK: run the first stage and inspect its F-statistic.
first = sm.OLS(df["education"], sm.add_constant(df[["distance"]])).fit()
print(f"\nfirst-stage F        = {first.fvalue:.1f}   (want > 10)")

# TWO-STAGE LEAST SQUARES, done manually so the mechanics are visible.
df["education_hat"] = first.fittedvalues
second = sm.OLS(df["log_wage"],
                sm.add_constant(df[["education_hat"]])).fit()
print(f"manual 2SLS          = {second.params['education_hat']:.4f}")
# Note: the manual standard errors here are WRONG, because the second
# stage ignores that education_hat was itself estimated. Use a proper
# IV routine for inference.

iv = IV2SLS.from_formula(
    "log_wage ~ 1 + [education ~ distance]", df).fit()
print(f"IV2SLS               = {iv.params['education']:.4f}  "
      f"se={iv.std_errors['education']:.4f}")

# WEAK INSTRUMENT: the same code with a nearly irrelevant instrument.
df["noise_iv"] = rng.normal(0, 1, n)
weak_first = sm.OLS(df["education"],
                    sm.add_constant(df[["noise_iv"]])).fit()
weak_iv = IV2SLS.from_formula(
    "log_wage ~ 1 + [education ~ noise_iv]", df).fit()
print(f"\nweak first-stage F   = {weak_first.fvalue:.2f}   (far below 10)")
print(f"weak IV estimate     = {weak_iv.params['education']:.4f}  "
      f"se={weak_iv.std_errors['education']:.4f}")
# The estimate is wild and the standard error explodes. When it does
# NOT explode, a weak instrument is even more dangerous, because the
# result looks precise and is badly biased toward the OLS estimate.

Key Takeaways

• An instrument affects the treatment, reaches the outcome only through it, and is unrelated to the confounders.
• Only relevance is testable; exclusion and independence require a substantive argument.
• Two-stage least squares keeps only the treatment variation the instrument explains, discarding the confounded part.
• A weak instrument produces badly biased estimates, so check that the first-stage F-statistic exceeds 10.

🧪 Practice

1. Simulate a confounded relationship with a valid instrument and compare ordinary least squares to two-stage least squares.
2. Weaken the instrument gradually and plot how the estimate degrades.
3. Interview: Someone proposes using a user's device type as an instrument for app usage when estimating its effect on spending. Is it valid? (Hint: ask whether device type could affect spending any other way.)

Difference-in-Differences

Often a treatment arrives at a specific time for a specific group: a policy changes in one region, a feature launches for one segment. You cannot randomize, but you have observations before and after, for both the affected group and an unaffected one. Difference-in-differences exploits that structure.

Comparing the treated group before and after is confounded by anything else that changed over time. Comparing treated and untreated after the change is confounded by pre-existing differences between the groups. Taking the difference of those two differences cancels both problems: time-invariant differences between groups, and group-invariant changes over time.

Everything rests on the parallel trends assumption: absent the treatment, the two groups' outcomes would have moved in parallel. The levels may differ, only the trends must match. This is not testable directly, because it concerns a counterfactual, but it can be supported by showing the trends were parallel in several pre-treatment periods.

Two practical warnings. Check for anticipation, where units change behaviour before the official start. And when treatment rolls out at different times for different units, the standard two-way fixed effects estimator is biased, a result established relatively recently; use one of the modern staggered-adoption estimators instead.

  THE TWO DIFFERENCES

   outcome
      ^                                        * treated, actual
      |                                    ,-'
      |                                ,-'
      |  treated            ,,--''''''       <- the gap at the end
      |            ,,--''''                     is NOT the effect
      |     ,,--''                        o counterfactual: where the
      |  --'                     ,-'-'-'-'  treated group WOULD have
      |                    ,-'-'            gone on the parallel trend
      |  control     ,-'-'
      |        ,-'-'
      +----------------|------------------------> time
                   treatment

   EFFECT = (treated_after - treated_before)
          - (control_after - control_before)

   The control's change estimates what would have happened anyway.

  THE 2x2 TABLE

                  before    after    difference
    treated         100      140        +40
    control          80       95        +15
                                    ----------
    DiD estimate = 40 - 15 =  +25

  PARALLEL TRENDS, CHECKED IN THE PRE-PERIOD

   good                          bad
    treated  ,-'                  treated    ,-'
    control ,-'                   control  -'
            parallel before               already diverging
            -> assumption credible        -> DiD will be biased

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rng = np.random.default_rng(0)

# Panel: 200 stores over 12 months. Treatment starts in month 7 for
# half the stores. Treated stores are BIGGER on average (a level
# difference), but both groups share the same seasonal trend.
n_stores, n_months = 200, 12
TREATMENT_MONTH = 7
TRUE_EFFECT = 30.0

rows = []
for store in range(n_stores):
    treated_store = store < n_stores // 2
    base = 500 + (80 if treated_store else 0) + rng.normal(0, 40)
    for month in range(1, n_months + 1):
        trend = 6 * month                                 # shared over time
        post = month >= TREATMENT_MONTH
        effect = TRUE_EFFECT if (treated_store and post) else 0
        rows.append({
            "store": store, "month": month,
            "treated": int(treated_store), "post": int(post),
            "sales": base + trend + effect + rng.normal(0, 25),
        })
df = pd.DataFrame(rows)

# The 2x2 table, computed directly.
cell = df.groupby(["treated", "post"])["sales"].mean().unstack()
did_manual = ((cell.loc[1, 1] - cell.loc[1, 0])
              - (cell.loc[0, 1] - cell.loc[0, 0]))
print(cell.round(1))
print(f"\nmanual DiD        = {did_manual:.2f}   (true {TRUE_EFFECT})")

# The regression form, which also gives standard errors. The
# interaction term IS the difference-in-differences estimate.
model = smf.ols("sales ~ treated + post + treated:post", data=df).fit(
    cov_type="cluster", cov_kwds={"groups": df["store"]})
# Cluster standard errors by store: observations from the same store
# across months are NOT independent, and ignoring that understates
# the uncertainty badly.
print(f"regression DiD    = {model.params['treated:post']:.2f} "
      f"(se {model.bse['treated:post']:.2f})")

# TEST PARALLEL TRENDS with an event-study specification: one
# coefficient per month relative to treatment.
df["rel"] = df["month"] - TREATMENT_MONTH
event = smf.ols("sales ~ C(rel) * treated", data=df).fit(
    cov_type="cluster", cov_kwds={"groups": df["store"]})
print("\npre-treatment interaction terms (should be near zero):")
for term, value in event.params.items():
    if "treated" in term and ":" in term and "-" in term:
        print(f"  {term:<30} {value:>8.2f}")
# Pre-treatment coefficients near zero support parallel trends.
# A visible pre-trend means the assumption fails and the DiD estimate
# is absorbing that divergence as if it were the effect.

Key Takeaways

• Difference-in-differences cancels time-invariant group differences and group-invariant time changes at once.
• Parallel trends is the critical assumption: the groups' levels may differ, but their trends must match absent treatment.
• Support the assumption with an event study showing flat pre-treatment coefficients.
• Cluster standard errors by unit, and avoid plain two-way fixed effects when treatment timing is staggered.

🧪 Practice

1. Build a panel dataset and estimate a difference-in-differences effect both manually and by regression.
2. Introduce a divergent pre-trend and show how it biases the estimate.
3. Interview: A region raised its minimum wage and employment fell 2% while the neighbouring region's fell 3%. What is the estimated effect? (Hint: subtract the two changes, then ask what the comparison assumes.)

Regression Discontinuity

Some treatments are assigned by a strict threshold on a continuous variable. Scholarships go to students above a test score. Credit goes to applicants above a score cutoff. Free shipping applies above a basket value. This rule creates a natural experiment at the boundary.

The insight is that units just barely above and just barely below the cutoff are, in every respect other than treatment, essentially identical. Whether a student scored 69 or 70 is largely noise. Yet one receives the treatment and the other does not. Comparing outcomes in a narrow window around the cutoff therefore approximates a randomized experiment.

Sharp designs assign treatment deterministically at the cutoff. Fuzzy designs merely change the probability of treatment, and are handled by using the threshold as an instrument.

Two checks are mandatory. First, test for manipulation: if units can precisely control which side of the cutoff they land on, sorting destroys the comparability. A density test looks for a suspicious pile-up just above the threshold. Second, run placebo tests at fake cutoffs, where no discontinuity should appear.

The important limitation is that the estimate is local. It identifies the effect for units near the threshold and says nothing about units far from it, which may be the units you care about most.

  THE DISCONTINUITY IS THE ESTIMATE

   outcome
      ^                                o o
      |                            o o      <- treated side
      |                        o o
      |                    o              }  the JUMP at the cutoff
      |                    |              }  is the causal effect
      |                  . |
      |              . .   |
      |          . .       |                  <- untreated side
      |      . .           |
      +----------------------------------> running variable
                        cutoff

   Units just left and just right differ only in treatment.
   Far from the cutoff they differ in many ways, which is why the
   estimate is LOCAL to the threshold.

  MANIPULATION CHECK: density of the running variable

   clean                          manipulated
    ___________                    _______  ||
   /           \                  /       \ ||___
   -------|------                 --------|-----
        cutoff                        cutoff
   smooth through the cutoff      pile-up just above -> units are
                                  sorting themselves; comparability
                                  is destroyed

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rng = np.random.default_rng(0)
n = 5000
CUTOFF = 70.0
TRUE_EFFECT = 8.0

# Sharp design: a scholarship for everyone scoring at or above 70.
score = rng.normal(68, 12, n)
treated = (score >= CUTOFF).astype(int)
# The outcome depends smoothly on score, PLUS a jump for the treated.
outcome = (30 + 0.45 * score + TRUE_EFFECT * treated + rng.normal(0, 6, n))

df = pd.DataFrame({"score": score, "treated": treated, "outcome": outcome})
df["centered"] = df["score"] - CUTOFF

# Naive comparison of all treated vs all untreated is badly confounded,
# because the treated simply have higher scores.
naive = df.loc[df.treated == 1, "outcome"].mean() - \
        df.loc[df.treated == 0, "outcome"].mean()
print(f"true effect            = {TRUE_EFFECT:.2f}")
print(f"naive all-units diff   = {naive:.2f}   confounded by score\n")

# LOCAL LINEAR RD: fit within a bandwidth, allowing different slopes
# on each side. The `treated` coefficient is the jump at the cutoff.
for bandwidth in (20, 10, 5, 3):
    window = df[np.abs(df["centered"]) <= bandwidth]
    fit = smf.ols("outcome ~ centered * treated", data=window).fit()
    print(f"bandwidth {bandwidth:>2}: estimate = {fit.params['treated']:6.3f} "
          f"(se {fit.bse['treated']:.3f}, n={len(window)})")
# Narrow bandwidths are less biased and noisier; wide ones are the
# reverse. The estimate should be stable across sensible bandwidths.
# If it moves a lot, the functional form is doing the work, not the
# discontinuity.

# MANIPULATION CHECK: is the running variable's density smooth?
just_below = ((df.centered >= -2) & (df.centered < 0)).sum()
just_above = ((df.centered >= 0) & (df.centered < 2)).sum()
print(f"\ndensity just below cutoff: {just_below}")
print(f"density just above cutoff: {just_above}")
ratio = just_above / just_below
print(f"ratio = {ratio:.3f}   (a big imbalance suggests sorting)")
# In production use a McCrary density test or `rdrobust`'s rddensity.

# PLACEBO TEST: there should be NO jump at a fake cutoff.
for fake in (60.0, 80.0):
    tmp = df.copy()
    tmp["c2"] = tmp["score"] - fake
    tmp["t2"] = (tmp["score"] >= fake).astype(int)
    w = tmp[np.abs(tmp["c2"]) <= 8]
    f = smf.ols("outcome ~ c2 * t2", data=w).fit()
    print(f"placebo cutoff {fake}: estimate = {f.params['t2']:+.3f} "
          f"(p={f.pvalues['t2']:.3f})   expect ~0")
# A significant jump at a fake cutoff means the specification is
# picking up curvature rather than a real discontinuity.

Key Takeaways

• A threshold rule creates a natural experiment, because units just either side of it are otherwise comparable.
• Estimate the jump with a local linear fit inside a bandwidth, allowing different slopes on each side.
• Check for manipulation of the running variable and run placebo tests at fake cutoffs.
• The estimate is local to the threshold and does not generalize to units far from it.

🧪 Practice

1. Simulate a sharp design and estimate the effect at three bandwidths, commenting on the stability.
2. Run placebo tests at two fake cutoffs and confirm no effect appears.
3. Interview: A scholarship is awarded above a test score and you find a large effect on graduation. What would make you doubt it? (Hint: ask whether anyone can retake the test or appeal a borderline score.)

---

<a id="4-programming-and-tooling"></a>