2.1 Data Generating Processes and Distributions

Why we simulate: playing god with data

In real data the truth is invisible: you never observe the true effect of a programme, only data produced by it. In a simulation you write the rules of the world yourself, so you know the true answer and can check whether a method recovers it.

That is the teaching strategy of this course: code over calculus. Instead of deriving the variance of an estimator, you will draw thousands of samples and look at how the estimates spread. Instead of proving that a regression is biased, you will build a world with a confounder and watch the bias appear.

A data generating process (DGP) is that rulebook: which variables exist, how each one is produced, and which ones cause which.

Setup and the random number generator

!pip install -q numpy scipy polars matplotlib seaborn
import numpy as np
import scipy.stats as stats
import polars as pl
import matplotlib.pyplot as plt
import seaborn as sns

rng = np.random.default_rng(101)     # a seeded generator: same numbers every run
sns.set_theme(style="whitegrid")
print(rng.normal(size=3))

We use NumPy's modern default_rng. A seed makes the "random" numbers reproducible, so your results match your classmates' and your future self's.

Building a public-health DGP

We simulate a population of 50,000 patients. The true parameters are set by us, and in the real world they would be unknown.

POPULATION_SIZE = 50_000
TRUE_MEAN_AGE = 45.0
TRUE_SD_AGE = 12.0
TRUE_COMORBIDITY_RATE = 0.25

age = np.clip(rng.normal(TRUE_MEAN_AGE, TRUE_SD_AGE, POPULATION_SIZE), 18, 90)
comorbidity = rng.binomial(n=1, p=TRUE_COMORBIDITY_RATE, size=POPULATION_SIZE)
visits = rng.poisson(lam=2 + 2 * comorbidity)          # GP visits per year

population = pl.DataFrame({
    "patient_id": np.arange(POPULATION_SIZE),
    "age": age,
    "has_comorbidity": comorbidity,
    "gp_visits": visits,
})
print(population.head())
print(population.select(
    pl.col("age").mean().round(2).alias("mean_age"),
    pl.col("has_comorbidity").mean().round(3).alias("comorbidity_rate"),
    pl.col("gp_visits").mean().round(2).alias("mean_visits"),
))

Notice gp_visits depends on has_comorbidity: we have already written a small causal rule into the DGP.

fig, axes = plt.subplots(1, 3, figsize=(14, 4))
sns.histplot(population["age"].to_numpy(), bins=40, ax=axes[0]).set(title="Age (normal, clipped)")
sns.countplot(x=population["has_comorbidity"].to_numpy(), ax=axes[1]).set(title="Comorbidity (Bernoulli)")
sns.countplot(x=population["gp_visits"].to_numpy(), ax=axes[2]).set(title="GP visits (Poisson)")
plt.tight_layout(); plt.show()

Distributions you will meet all course long

scipy.stats gives the theoretical side of the same distributions, for example the exact probability that a normal variable exceeds a value, so you can compare simulated and theoretical answers.

theory = 1 - stats.norm.cdf(65, loc=TRUE_MEAN_AGE, scale=TRUE_SD_AGE)
simulated = (population["age"] > 65).mean()
print(f"P(age > 65): theory {theory:.4f} | simulation {simulated:.4f}")

Population vs sample. The 45.0 we typed is a parameter. The mean of the simulated age column is a statistic: close to 45, but not exactly 45, and slightly lower here because of clipping. Real studies only ever see statistics from samples. Lesson 2.2 is about how far they can stray.

Exercise: an income distribution

income = rng.lognormal(mean=10.8, sigma=0.6, size=POPULATION_SIZE)
print(f"Mean income:   ${income.mean():,.0f}")
print(f"Median income: ${np.median(income):,.0f}")
# The long right tail pulls the mean above the median. For a 'typical' household, report the median.

Checklist and recap

• I can define a DGP and write one in NumPy with a fixed seed.
• I know which distribution to use for treatments, counts, incomes and noise.
• I can tell a population parameter from a sample statistic.