7. Exploratory Data Analysis and Visualization

Exploratory data analysis is the process of looking at data closely enough to know what questions it can answer, before committing to a model or a claim. This chapter covers the systematic workflow for exploring one variable and then many, the design principles that make a chart honest and readable, the practical tooling for producing charts, and the reporting formats that carry a finding to the people who will act on it. The recurring theme is that exploration generates hypotheses rather than confirming them, and that a chart is a claim about the data which the reader will believe.

<a id="71-exploratory-analysis-workflow"></a>

7.1 Exploratory Analysis Workflow

Exploration has a sensible order, from one variable to many, and each step answers a question that constrains what you do next; this section covers that progression and the traps at each stage.

Univariate Exploration

Univariate exploration means looking at one variable at a time. It is the first step because almost every later technique assumes something about individual distributions, and those assumptions are cheap to check and expensive to get wrong.

The question a single variable answers is: what does a typical value look like, how much do values vary, and what shape do they form? Those three questions map onto centre, spread, and shape, and the reason all three matter is that no one of them determines the others. Two variables can share a mean and differ completely.

Start with centre. The mean is the balance point and responds to every value, including absurd ones. The median is the middle value and does not. When they disagree substantially, that disagreement is itself the finding: it tells you the distribution is skewed, and it tells you which summary to report. The mode matters for categorical data and for spotting sentinel values that repeat implausibly often.

Then spread. The standard deviation is in the variable's own units and is interpretable when the distribution is roughly symmetric. The interquartile range, the distance from the 25th to the 75th percentile, describes the middle half and ignores the tails. Reporting a mean without a spread is nearly meaningless, because it hides whether values cluster tightly or scatter wildly.

Then shape, which is where visualisation becomes essential rather than convenient. Numbers summarise; shapes reveal. Specifically you are looking for:

• Skew, a longer tail on one side, which tells you the mean is misleading.
• Modality, the number of peaks. A second peak almost always means two populations are mixed in one column, and that is a structural finding: you probably need to split the data rather than model it as one group.
• Bounds and spikes, a hard floor or ceiling, or a pile-up at one value. These indicate censoring, clipping, or a default value.
• Gaps, ranges with no observations at all, which often mean a filter was applied upstream that nobody told you about.

The essential tools are the histogram, which bins values and shows shape, and the box plot, which shows quartiles compactly and is ideal for comparing many groups at once. Both have a weakness worth knowing. A histogram's appearance depends heavily on bin width, so a single histogram can create or hide a second peak; always look at two or three bin widths. A box plot hides modality entirely, because a bimodal distribution and a uniform one can produce identical boxes. The violin plot and the empirical cumulative distribution function avoid that, the latter being the most honest single view of a distribution because it involves no binning choice at all.

For categorical variables the questions are different: how many distinct levels, how concentrated are the counts, and is there a long tail of rare values. A bar chart sorted by frequency answers all three at a glance.

import numpy as np
import pandas as pd

rng = np.random.default_rng(0)

# Four distributions that share a mean of ~100 and differ in every other way.
n = 2000
data = pd.DataFrame({
    "symmetric": rng.normal(100, 15, n),
    "right_skew": rng.lognormal(np.log(100) - 0.5 * 0.6**2, 0.6, n),
    "bimodal": np.where(rng.random(n) < 0.5,
                        rng.normal(70, 8, n), rng.normal(130, 8, n)),
    "with_spike": np.where(rng.random(n) < 0.25, 100.0, rng.normal(100, 25, n)),
})

summary = data.agg(["mean", "median", "std", "skew", "kurtosis"]).T
summary["iqr"] = data.quantile(0.75) - data.quantile(0.25)
summary["mean_minus_median"] = summary["mean"] - summary["median"]
print(summary.round(2).to_string())
# All four means are near 100. Only skew, kurtosis, and the IQR hint that
# these are different shapes, and NONE of them reveals the bimodality.

# --- THE DISAGREEMENT BETWEEN MEAN AND MEDIAN IS THE FINDING ---
for col in data.columns:
    s = data[col]
    pct = (s < s.mean()).mean()
    print(f"{col:12s} {pct:.0%} of values fall below the mean")
# For a symmetric variable this is ~50%. Anything far from 50% says the mean
# is not a typical value and you should report the median instead.

# --- A TEXT HISTOGRAM: enough to see shape without a plotting library. ---
def ascii_hist(s: pd.Series, bins: int = 20, width: int = 46) -> str:
    counts, edges = np.histogram(s.dropna(), bins=bins)
    scale = width / max(counts.max(), 1)
    out = []
    for c, lo in zip(counts, edges[:-1]):
        out.append(f"  {lo:7.1f} | {'#' * int(round(c * scale)):<{width}} {c}")
    return "\n".join(out)

print("\nbimodal, 20 bins:\n" + ascii_hist(data["bimodal"]))

# --- BIN WIDTH CHANGES THE STORY. Always look at more than one. ---
print("\nthe same column at three bin counts:")
for b in (5, 15, 60):
    counts, _ = np.histogram(data["bimodal"], bins=b)
    # Count local maxima as a crude proxy for "how many peaks do I see?"
    peaks = sum(1 for i in range(1, len(counts) - 1)
                if counts[i] > counts[i-1] and counts[i] > counts[i+1])
    print(f"  bins={b:>3}  apparent peaks: {peaks}")
# Too few bins smooth the two peaks into one; too many turn noise into peaks.

# Rules of thumb for a starting bin count, then adjust by eye.
x = data["bimodal"].dropna()
sturges = int(np.ceil(np.log2(len(x)) + 1))
fd_width = 2 * (x.quantile(.75) - x.quantile(.25)) / len(x) ** (1/3)
print(f"\nSturges suggests {sturges} bins;"
      f" Freedman-Diaconis suggests {int(np.ceil(x.ptp() / fd_width))}")

# --- WHAT A BOX PLOT HIDES ---
def five_number(s: pd.Series) -> dict:
    q1, q2, q3 = s.quantile([.25, .5, .75])
    return {"min": s.min(), "q1": q1, "median": q2, "q3": q3, "max": s.max()}

box = pd.DataFrame({c: five_number(data[c]) for c in ["bimodal", "symmetric"]})
print("\nfive-number summaries:\n", box.round(1).to_string())
print("these two boxes look similar; one has two peaks and the other one.")

# --- ECDF: no binning choice, so no way to mislead yourself. ---
def ascii_ecdf(s: pd.Series, rows: int = 9) -> str:
    qs = np.linspace(0, 1, rows)
    vals = s.quantile(qs)
    out = []
    for q, v in zip(qs, vals):
        out.append(f"  p{int(q*100):>3} = {v:8.1f}  {'*' * int(q * 40)}")
    return "\n".join(out)

print("\nECDF of right_skew:\n" + ascii_ecdf(data["right_skew"]))
# Read it as: "what fraction of values are at or below x?" A steep section is
# a dense region; a flat section is a gap.

# --- SPIKE DETECTION: a repeated value is a code, not a measurement. ---
for col in data.columns:
    vc = data[col].round(3).value_counts()
    if vc.iloc[0] > len(data) * 0.02:
        print(f"\n{col}: value {vc.index[0]} occurs {vc.iloc[0]} times"
              f" ({vc.iloc[0]/len(data):.0%}) -> suspect a default or sentinel")

# --- CATEGORICAL: levels, concentration, and the tail. ---
cat = pd.Series(rng.choice(list("abcdefghij"), 1500,
                           p=np.array([40, 25, 15, 8, 5, 3, 2, 1, 0.7, 0.3])/100))
counts = cat.value_counts()
print("\ncategory counts:\n", counts.to_string())
print(f"\ntop 2 levels cover {counts.head(2).sum()/len(cat):.0%} of rows")
print(f"levels under 1% of rows: {(counts/len(cat) < 0.01).sum()}")
# Concentration decides whether one-hot encoding is viable and whether a
# "rare" bucket is needed. See chapter 6 on high-cardinality features.

# --- A REUSABLE UNIVARIATE REPORT ---
def univariate(df: pd.DataFrame) -> pd.DataFrame:
    rows = {}
    for col in df.select_dtypes("number").columns:
        s = df[col].dropna()
        q1, q3 = s.quantile([.25, .75])
        rows[col] = {
            "n": len(s), "mean": s.mean(), "median": s.median(),
            "std": s.std(), "iqr": q3 - q1, "skew": s.skew(),
            "p01": s.quantile(.01), "p99": s.quantile(.99),
            "pct_below_mean": (s < s.mean()).mean(),
            "zeros": int((s == 0).sum()),
        }
    return pd.DataFrame(rows).T.round(2)

print("\n", univariate(data).to_string())

  FOUR DISTRIBUTIONS, ONE MEAN

  symmetric            right skew           bimodal            spike at 100
       ####                #                  ##    ##              #
      ######              ###                ####  ####            ##
     ########            #####               ##########           ####
    ##########          #######_____         ##########          ######
  --------------      ----------------     --------------     ------------
       ^mean            ^mean                   ^mean              ^mean
     mean=median    mean > median          mean is between     mean is a
                                           the two peaks       default value
                                           and describes
                                           nobody

  WHAT EACH VIEW SHOWS

  view        centre  spread  skew  modality  outliers  binning choice
  ---------------------------------------------------------------------
  summary     yes     yes     hint  NO        hint      none
  histogram   yes     yes     yes   yes       yes       YES, big effect
  box plot    yes     yes     yes   NO        yes       none
  violin      yes     yes     yes   yes       partly    bandwidth
  ECDF        yes     yes     yes   partly    yes       NONE

  THE ORDER TO LOOK

    1. mean vs median      -> is the mean usable at all?
    2. histogram, 2-3 bin widths -> shape, peaks, gaps
    3. p01 and p99         -> tails, without being dragged by one value
    4. most frequent value -> sentinels and defaults
    5. count of zeros      -> is zero a measurement or a missing marker?

Key Takeaways

• Centre, spread, and shape are independent questions, and summary statistics answer only the first two.
• A large gap between the mean and the median tells you the mean is not a typical value and the median should be reported.
• A second peak almost always means two populations are mixed in one column, which is a structural finding rather than a modelling detail.
• Histogram appearance depends on bin width, so inspect several; box plots hide modality entirely, and the ECDF involves no binning choice.
• An implausibly frequent value is a sentinel or default, not a measurement.

🧪 Practice

1. Build four columns that share a mean but differ in skew, modality, and the presence of a spike, then show that a table of means and standard deviations cannot distinguish them.
2. Plot the same bimodal column with 5, 15, and 60 bins and describe how the apparent number of peaks changes.
3. Write a univariate report function returning centre, spread, tail percentiles, and the fraction of values below the mean for every numeric column, and run it on a real dataset.
4. Interview: A stakeholder reports that average session length rose from 4 to 6 minutes and concludes engagement improved. What do you check before agreeing? (Hint: which single view would tell you whether the typical session changed at all?)

Bivariate Relationships

Once you know what each variable looks like alone, the next question is how they move together. Bivariate analysis is where most findings actually come from, because a relationship between two variables is a candidate explanation, and an explanation is what a stakeholder wants.

The right technique depends entirely on the pair of types involved, and this is the single most useful thing to have memorised.

Pair	Primary view	Primary summary
numeric vs numeric	scatter plot	correlation, fitted line
numeric vs categorical	box or violin per group	group means, effect size
categorical vs categorical	grouped bar, heatmap	contingency table, rates
numeric vs time	line chart	trend, seasonality

Numeric versus numeric belongs in a scatter plot, and the scatter plot is the most information-dense chart there is. It shows the direction and strength of a relationship, its functional form, whether the spread changes across the range, and which individual points do not fit. No summary statistic carries all of that, which is why you plot before you compute.

The specific failure to guard against is trusting a correlation coefficient without looking. Anscombe's quartet is four datasets with identical means, variances, and correlations, and wildly different shapes. It exists to make one point: the number is a projection of the picture, and projections lose information.

For large datasets a scatter plot overplots, meaning points pile on top of each other until the dense region is a solid mass and you can no longer see where the mass is. The fixes are transparency, sampling, or switching to a binned view such as a hexbin or a two-dimensional histogram, which counts points per cell and shows density honestly.

Numeric versus categorical is a comparison of distributions across groups. Box plots side by side are the workhorse. Two habits improve them enormously: sort the groups by their median rather than alphabetically, because ordering is free information, and show the group sizes, because a box drawn from six observations deserves less trust than one drawn from six thousand and the chart does not say which is which.

Categorical versus categorical is a contingency table, and the important move is deciding which way to normalise. Raw counts are dominated by whichever category is common. Row percentages answer "of this group, what fraction did X?" and column percentages answer "of those who did X, what fraction were in this group?" Those are different questions with different answers, and presenting one while the reader assumes the other is a common source of confusion.

Numeric versus time is a line chart, and the thing to remember is that time imposes order, so the points are connected for a reason. Never connect points across a gap without marking it, and never use a line chart for a categorical x-axis, because the line implies an interpolation that does not exist.

import numpy as np
import pandas as pd

rng = np.random.default_rng(3)

# --- ANSCOMBE'S QUARTET: identical statistics, four different shapes. ---
anscombe = pd.DataFrame({
    "x1": [10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5],
    "y1": [8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68],
    "x2": [10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5],
    "y2": [9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74],
    "x3": [10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5],
    "y3": [7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73],
    "x4": [8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8],
    "y4": [6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89],
})
rows = []
for i in (1, 2, 3, 4):
    x, y = anscombe[f"x{i}"], anscombe[f"y{i}"]
    slope, intercept = np.polyfit(x, y, 1)
    rows.append({"set": i, "mean_x": x.mean(), "mean_y": round(y.mean(), 2),
                 "sd_y": round(y.std(ddof=1), 2),
                 "pearson_r": round(x.corr(y), 3),
                 "slope": round(slope, 2), "intercept": round(intercept, 2)})
print(pd.DataFrame(rows).to_string(index=False))
print("\nidentical to two decimals, and the four shapes are:")
print("  1: a genuine linear relationship")
print("  2: a clean curve, so a line is the wrong model")
print("  3: perfectly linear except one outlier that sets the slope")
print("  4: no relationship at all; one leverage point creates the whole fit")

# --- FUNCTIONAL FORM: a correlation near zero does not mean independent. ---
x = rng.uniform(-3, 3, 800)
pairs = {
    "linear": 2 * x + rng.normal(0, 1, 800),
    "quadratic": x ** 2 + rng.normal(0, 1, 800),       # symmetric, so r ~ 0
    "sinusoidal": np.sin(2 * x) + rng.normal(0, .2, 800),
    "heteroscedastic": 2 * x + rng.normal(0, 1, 800) * np.abs(x),
}
print("\nPearson r vs a rank correlation vs mutual information:")
from sklearn.feature_selection import mutual_info_regression
for name, y in pairs.items():
    r = np.corrcoef(x, y)[0, 1]
    rho = pd.Series(x).corr(pd.Series(y), method="spearman")
    mi = mutual_info_regression(x.reshape(-1, 1), y, random_state=0)[0]
    print(f"  {name:16s} r {r:+.3f}   spearman {rho:+.3f}   MI {mi:.3f}")
# The quadratic pair has r near zero and high mutual information: strongly
# related, not linearly. Reporting "no correlation" here would be wrong.

# --- OVERPLOTTING: what a dense scatter hides. ---
big_x = rng.normal(0, 1, 60_000)
big_y = big_x * 0.6 + rng.normal(0, 0.8, 60_000)

def hexbin_text(xs, ys, nx=28, ny=11) -> str:
    """A 2-D histogram rendered as density characters."""
    H, _, _ = np.histogram2d(xs, ys, bins=[nx, ny])
    ramp = " .:-=+*#%@"
    hi = H.max()
    out = []
    for j in range(ny - 1, -1, -1):
        line = "".join(ramp[min(int(H[i, j] / hi * (len(ramp) - 1)), len(ramp) - 1)]
                       for i in range(nx))
        out.append("  " + line)
    return "\n".join(out)

print("\n60,000 points as a density grid (a plain scatter would be a blob):")
print(hexbin_text(big_x, big_y))
# Transparency, sampling, or binning. A solid mass of ink is not a finding.

# --- NUMERIC VS CATEGORICAL: sort by median, and show n. ---
groups = pd.DataFrame({
    "plan": rng.choice(["free", "basic", "pro", "enterprise"], 1200,
                       p=[.55, .25, .15, .05]),
})
base = {"free": 20, "basic": 45, "pro": 90, "enterprise": 300}
groups["spend"] = [base[p] * rng.lognormal(0, .4) for p in groups["plan"]]

by_plan = groups.groupby("plan")["spend"].agg(
    n="size", median="median", q1=lambda s: s.quantile(.25),
    q3=lambda s: s.quantile(.75), mean="mean")
print("\n", by_plan.sort_values("median").round(1).to_string())
print("\nsorted by median, not alphabetically, and n is visible so the reader")
print("knows the enterprise box rests on far fewer observations.")

# A standardised effect size makes group gaps comparable across variables.
a = groups.loc[groups["plan"] == "pro", "spend"]
b = groups.loc[groups["plan"] == "basic", "spend"]
pooled = np.sqrt(((len(a)-1)*a.var(ddof=1) + (len(b)-1)*b.var(ddof=1))
                 / (len(a)+len(b)-2))
print(f"\nCohen's d, pro vs basic: {(a.mean()-b.mean())/pooled:.2f}")

# --- CATEGORICAL VS CATEGORICAL: which way you normalise is the question. ---
ct = pd.DataFrame({
    "channel": rng.choice(["ads", "organic", "referral"], 2000, p=[.5, .35, .15]),
})
p_convert = {"ads": .05, "organic": .12, "referral": .20}
ct["converted"] = [rng.random() < p_convert[c] for c in ct["channel"]]

counts = pd.crosstab(ct["channel"], ct["converted"])
print("\nraw counts (dominated by the biggest channel):\n", counts.to_string())
print("\nrow % -- 'of this channel, what fraction converted?'")
print((pd.crosstab(ct["channel"], ct["converted"], normalize="index") * 100)
      .round(1).to_string())
print("\ncolumn % -- 'of converters, which channel did they come from?'")
print((pd.crosstab(ct["channel"], ct["converted"], normalize="columns") * 100)
      .round(1).to_string())
# Ads have the LOWEST conversion rate and supply the MOST converters. Both
# statements are true and they answer different questions.

# A chi-squared test asks whether the pattern is more than sampling noise.
from scipy.stats import chi2_contingency
chi2, p, dof, _ = chi2_contingency(counts)
print(f"\nchi-squared {chi2:.1f}, dof {dof}, p {p:.2e}")

  PICK THE VIEW FROM THE TYPE PAIR

              | numeric              | categorical
  ------------+----------------------+---------------------------
  numeric     | scatter (hexbin if   | box/violin per group,
              | dense), + fitted     | sorted by median, n shown
              | line, + residuals    |
  ------------+----------------------+---------------------------
  categorical | (transpose of above) | contingency table,
              |                      | grouped bar or heatmap
  ------------+----------------------+---------------------------
  time        | line; never connect across a gap unmarked

  WHY YOU PLOT BEFORE YOU COMPUTE

  set 1          set 2          set 3          set 4
   linear         curved         outlier        leverage point
     .  .           . . .          .               .
   .  .           .     .        .               |
  .  .           .       .      . . . . .        | (all x equal)
  ------------   ------------   ------------   ------------
  r = 0.816      r = 0.816      r = 0.816      r = 0.816

  All four share mean, variance, correlation, and fitted line.
  Only one of them should be summarised by that line.

  ROW % VS COLUMN %

    of this channel, what      of the converters, where
    fraction converted?        did they come from?
    ads       5%               ads       31%
    organic  12%               organic   52%
    referral 20%               referral  17%

    "Ads convert worst" and "ads bring the most conversions" are both true.
    State which normalisation you used, every time.

Key Takeaways

• The type pair decides the view: scatter for numeric pairs, grouped box plots for numeric against categorical, a contingency table for two categoricals.
• A correlation coefficient is a lossy projection of a scatter plot, so plot first; Anscombe's quartet shares every summary statistic and no shape.
• A near-zero linear correlation does not mean independence, as a symmetric curve demonstrates.
• Dense scatter plots overplot, so use transparency, sampling, or a binned density view instead of a solid mass of ink.
• For contingency tables, row and column percentages answer different questions, so always state which one you are showing.

🧪 Practice

1. Reproduce Anscombe's quartet, confirm the four correlations agree to two decimals, and describe what each scatter plot shows.
2. Generate a quadratic relationship, report the Pearson correlation and a rank correlation, and explain why one is near zero.
3. Build a contingency table of channel against conversion, then present it three ways: counts, row percentages, and column percentages, with a one-sentence takeaway for each.
4. Interview: A scatter plot of 200,000 points looks like a solid rectangle. Name three ways to recover the structure, with a trade-off for each. (Hint: you can change the opacity, the number of rows, or the mark.)

Multivariate Exploration

Two variables at a time is not enough, because the relationship between two variables can be created, hidden, or reversed by a third. Multivariate exploration is how you find that third variable, and it is the step that separates a finding that survives scrutiny from one that embarrasses you in a review.

The central phenomenon is confounding. A confounder is a variable that influences both of the others, producing an association between them that is not a direct effect. Ice cream sales and drowning deaths correlate because both rise in summer. Temperature is the confounder, and controlling for it makes the association disappear.

Its dramatic form is Simpson's paradox, where the direction of a relationship reverses when you condition on a group. The classic case is university admissions where each department admitted men and women at similar rates, or favoured women, while the aggregate showed a bias against women, because women applied disproportionately to more competitive departments. Both the aggregate and the per-department numbers were correct. Only one of them answered the question being asked.

The practical defence is to make conditioning a habit. Whenever you report a relationship, look at it again within the levels of your most important categorical variables. If it holds everywhere, you have something robust. If it reverses or vanishes somewhere, you have found the real story.

The techniques for seeing more than two variables at once:

Small multiples, also called faceting, repeat the same chart across levels of a third variable. This is the most reliable multivariate technique because it relies on position comparison, which humans do well, rather than on colour or size, which we do poorly. Faceting is covered in depth later in this chapter.

Encoding extra variables in a scatter plot through colour, size, and shape. Useful and strictly limited: colour handles roughly five to seven categories before adjacent hues blur, size communicates only large differences, and shape handles very few levels. Past that, facet.

The scatter plot matrix, a grid of every pairwise scatter, gives you all two-way relationships at once. It scales badly, becoming unreadable past about eight variables, but for a first pass on a modest dataset it is excellent value.

The correlation heatmap summarises all pairs numerically, which scales further at the cost of hiding every shape. Use it to decide which pairs to plot, never as the final word.

Parallel coordinates draw one line per observation across parallel axes, which reveals clusters and is the rare technique that handles many dimensions at once. It is dense and takes practice to read.

Dimensionality reduction, principal component analysis and neighbour embeddings such as t-SNE and UMAP, projects many dimensions down to two for plotting. These are powerful and easy to over-interpret: in a neighbour embedding, distances between clusters and cluster sizes are not meaningful, only local neighbourhoods are.

import numpy as np
import pandas as pd

rng = np.random.default_rng(7)

# --- SIMPSON'S PARADOX, built deliberately. ---
# Two departments. Each admits women at a HIGHER rate than men. Women apply
# mostly to the competitive department, so the aggregate reverses.
rows = []
for dept, n_m, n_w, rate_m, rate_w in [
        ("easy",        800,  200, 0.62, 0.68),
        ("competitive", 200,  800, 0.14, 0.20)]:
    for sex, n, rate in (("M", n_m, rate_m), ("W", n_w, rate_w)):
        admits = int(round(n * rate))
        rows += [{"dept": dept, "sex": sex, "admitted": 1}] * admits
        rows += [{"dept": dept, "sex": sex, "admitted": 0}] * (n - admits)
apps = pd.DataFrame(rows)

print("WITHIN each department, women are admitted at a higher rate:")
within = (apps.groupby(["dept", "sex"])["admitted"]
              .agg(rate="mean", n="size").round(3))
print(within.to_string())

print("\nAGGREGATED over departments, the direction reverses:")
overall = apps.groupby("sex")["admitted"].agg(rate="mean", n="size").round(3)
print(overall.to_string())

print("\nthe mechanism -- where each group applied:")
print((pd.crosstab(apps["sex"], apps["dept"], normalize="index") * 100)
      .round(0).to_string())
print("\nBoth tables are correct. The aggregate answers 'what happened';")
print("the per-department table answers 'was the process biased'.")

# --- CONFOUNDING: an association that a third variable fully explains. ---
n = 1500
temperature = rng.uniform(0, 35, n)
ice_cream = 20 + 3.0 * temperature + rng.normal(0, 12, n)
drownings = 1 + 0.25 * temperature + rng.normal(0, 1.2, n)
conf = pd.DataFrame({"temp": temperature, "ice_cream": ice_cream,
                     "drownings": drownings})

print(f"\nraw correlation ice cream vs drownings:"
      f" {conf['ice_cream'].corr(conf['drownings']):+.3f}")

# Condition on the confounder: within narrow temperature bands, it vanishes.
conf["band"] = pd.qcut(conf["temp"], 6)
within_band = conf.groupby("band", observed=True).apply(
    lambda g: g["ice_cream"].corr(g["drownings"]), include_groups=False)
print("\ncorrelation within each temperature band:")
print(within_band.round(3).to_string())
print(f"\nmean within-band correlation: {within_band.mean():+.3f}  <- ~0")

# The regression version: add the confounder and the coefficient collapses.
import statsmodels.api as sm
naive = sm.OLS(conf["drownings"],
               sm.add_constant(conf[["ice_cream"]])).fit()
adjusted = sm.OLS(conf["drownings"],
                  sm.add_constant(conf[["ice_cream", "temp"]])).fit()
print(f"\nice_cream coefficient alone:      {naive.params['ice_cream']:+.4f}"
      f"  (p {naive.pvalues['ice_cream']:.1e})")
print(f"ice_cream coefficient + temp:     {adjusted.params['ice_cream']:+.4f}"
      f"  (p {adjusted.pvalues['ice_cream']:.2f})")

# --- HOW MANY VARIABLES CAN ONE CHART CARRY? ---
print("""
  channel   comfortable capacity            degrades by
  -------------------------------------------------------------------
  x, y      continuous, high precision      overplotting
  colour    5-7 categories, or 1 ramp       adjacent hues blur; CVD
  size      3-4 distinguishable steps       area is judged poorly
  shape     4-5 glyphs                      small sizes become identical
  facet     10-20 small panels              panel size, then page space
  Past these, facet rather than adding another channel.""")

# --- SCATTER PLOT MATRIX: the pairwise summary behind the picture. ---
feat = pd.DataFrame({
    "a": rng.normal(0, 1, 600),
})
feat["b"] = 0.8 * feat["a"] + rng.normal(0, .6, 600)     # correlated with a
feat["c"] = feat["a"] ** 2 + rng.normal(0, .5, 600)      # curved in a
feat["d"] = rng.normal(0, 1, 600)                        # unrelated

pearson = feat.corr()
spearman = feat.corr(method="spearman")
print("\nPearson (linear only):\n", pearson.round(2).to_string())
print("\npairs where the rank correlation is much larger than Pearson,")
print("which flags a non-linear relationship worth plotting:")
diff = (spearman.abs() - pearson.abs()).round(2)
for i in feat.columns:
    for j in feat.columns:
        if i < j and diff.loc[i, j] > 0.1:
            print(f"  {i} vs {j}: pearson {pearson.loc[i,j]:+.2f},"
                  f" spearman {spearman.loc[i,j]:+.2f}")

# --- CONDITIONING AS A HABIT: check a relationship inside every segment. ---
def check_across_segments(df, x, y, by, min_n=30):
    out = []
    overall = df[x].corr(df[y])
    for level, g in df.groupby(by, observed=True):
        if len(g) >= min_n:
            out.append({by: level, "n": len(g), "corr": round(g[x].corr(g[y]), 3)})
    res = pd.DataFrame(out)
    res["sign_flips"] = np.sign(res["corr"]) != np.sign(overall)
    print(f"overall corr({x}, {y}) = {overall:+.3f}")
    return res

print("\n", check_across_segments(conf, "ice_cream", "drownings", "band")
        .to_string(index=False))
print("\nevery segment near zero while the overall is strongly positive is the")
print("signature of a confounder, not of a weak effect.")

  SIMPSON'S PARADOX, DRAWN

  WITHIN departments              AGGREGATED
  admit rate                      admit rate
   |  M   W                        |
  .7|  #   #                      .5|  #
  .6|  #   #   easy               .4|  #      M
  .5|                             .3|  #      #
  .2|          W   competitive    .2|  #      #   W
  .1|      M   #                  .1|  #      #   #
    +-------------                  +-------------
    women higher in BOTH            women lower OVERALL

  The mechanism is the application mix, not the admission rule:
    men   -> 80% easy, 20% competitive
    women -> 20% easy, 80% competitive

  CONFOUNDING, AS A DIAGRAM

           temperature
            /        \
           v          v
     ice cream ....  drownings
       (no arrow between them: the dotted line is the OBSERVED
        association, produced entirely by the common cause)

    Condition on temperature -> the dotted association disappears.

  THE MULTIVARIATE LADDER

    2 vars   scatter / box / crosstab
    3 vars   + colour, or facet by the third
    4 vars   + size or shape, or a 2-way facet grid
    5-8      scatter matrix, correlation heatmap, parallel coordinates
    8+       dimensionality reduction (and interpret it cautiously)

Key Takeaways

• A third variable can create, hide, or reverse a two-variable relationship, so conditioning on segments should be a habit rather than a special step.
• Simpson's paradox means both the aggregate and the conditional numbers can be correct while answering different questions.
• Faceting beats extra visual channels because position comparison is the one thing human vision does precisely.
• Colour carries five to seven categories, size three or four steps, and shape four or five glyphs; past that, facet.
• Correlation heatmaps decide which pairs to plot and never substitute for plotting them.

🧪 Practice

1. Construct a dataset exhibiting Simpson's paradox and present both the conditional and aggregate tables with a sentence on which answers which question.
2. Simulate a confounded pair, show the raw correlation, then show it collapsing within narrow bands of the confounder.
3. Write a function that recomputes a correlation within every level of a grouping column and flags any level whose sign differs from the overall.
4. Interview: Your A/B test shows the new checkout is worse overall, yet better on both mobile and desktop. What happened and what do you report? (Hint: what else changed between the two groups besides the variant?)

Correlation Analysis

Correlation measures how much two variables move together. It is the most used and most misused number in data work, so it is worth being precise about what each version measures and what none of them measure.

Pearson's r is the workhorse. It measures the strength of a linear association, ranges from -1 to 1, and is the covariance divided by the product of the standard deviations, which makes it scale-free. Its assumptions are the source of most misuse: it assumes the relationship is linear, it is sensitive to outliers, and it is meaningful for continuous variables.

Spearman's rho is Pearson's r computed on the ranks instead of the values. Because ranks discard the spacing, it measures whether the relationship is monotonic rather than straight, and it is robust to outliers and to any monotonic transformation. If a variable is heavily skewed, Spearman is usually the more honest choice.

Kendall's tau counts concordant and discordant pairs. It is more robust than Spearman on small samples and with ties, at a higher computational cost, and its values run smaller than Spearman's for the same data, so do not compare them directly.

Mutual information measures any statistical dependence, linear or not. It catches the relationships the others miss, at the cost of being unbounded above in its raw form, harder to interpret, and requiring more data.

For categorical associations, Cramér's V derives an interpretable zero-to-one measure from a chi-squared statistic, and the correlation ratio handles a numeric against a categorical variable.

Then the four warnings, each of which accounts for real mistakes.

Correlation is not causation. This is repeated so often that it has stopped landing, so state the alternatives explicitly: a correlation between A and B is consistent with A causing B, B causing A, a third variable causing both, selection in how the data was collected, or coincidence. Establishing which requires an experiment or an identification strategy, covered in the chapters on experimentation and causal inference.

Zero correlation is not independence. Pearson's r is zero for any symmetric non-monotonic relationship, and a U-shape is a strong dependence with r near zero.

Correlation is not comparable across different ranges. Restricting the range of a variable attenuates its correlation, sometimes dramatically. A correlation computed among admitted students will be lower than the same correlation in the full applicant pool, purely because of the selection. This is called range restriction and it makes correlations from different populations non-comparable.

Statistical significance is not strength. With a hundred thousand rows, a correlation of 0.01 has a tiny p-value and no practical meaning. Report the coefficient and a confidence interval, not just the p-value.

One more issue specific to exploration: when you compute a correlation matrix over fifty variables, you have computed 1,225 coefficients, and by chance alone dozens will look notable. This is the multiple comparisons problem, and treating the largest coefficient in a big matrix as a discovery is how spurious findings enter a report. Correct for it, or treat the matrix strictly as a source of hypotheses to test on fresh data.

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

rng = np.random.default_rng(11)
n = 600
x = rng.uniform(-3, 3, n)

shapes = {
    "linear":         2 * x + rng.normal(0, 1.0, n),
    "monotone curve": np.exp(x) + rng.normal(0, 1.0, n),
    "u-shape":        x ** 2 + rng.normal(0, 0.8, n),
    "sinusoid":       np.sin(2.2 * x) + rng.normal(0, 0.25, n),
    "no relation":    rng.normal(0, 1, n),
}

from sklearn.feature_selection import mutual_info_regression
print(f"{'shape':16s} {'pearson':>8s} {'spearman':>9s} {'kendall':>8s} {'MI':>6s}")
for name, y in shapes.items():
    r = stats.pearsonr(x, y).statistic
    rho = stats.spearmanr(x, y).statistic
    tau = stats.kendalltau(x, y).statistic
    mi = mutual_info_regression(x.reshape(-1, 1), y, random_state=0)[0]
    print(f"{name:16s} {r:+8.3f} {rho:+9.3f} {tau:+8.3f} {mi:6.3f}")
print("\nthe u-shape and the sinusoid have Pearson r near zero and high mutual")
print("information: strongly dependent, not linearly related.")

# --- ONE OUTLIER CAN MANUFACTURE OR DESTROY A CORRELATION ---
clean_x = rng.normal(0, 1, 60)
clean_y = rng.normal(0, 1, 60)                        # genuinely unrelated
print(f"\n60 unrelated points: r = {np.corrcoef(clean_x, clean_y)[0,1]:+.3f}")
dirty_x = np.append(clean_x, 12.0)
dirty_y = np.append(clean_y, 12.0)                    # one leverage point
print(f"after adding ONE point at (12, 12): r = "
      f"{np.corrcoef(dirty_x, dirty_y)[0,1]:+.3f}")
print(f"Spearman is unmoved: {stats.spearmanr(dirty_x, dirty_y).statistic:+.3f}"
      "  <- ranks cannot be dragged")

# --- CONFIDENCE INTERVALS: a coefficient without one is half a result. ---
def pearson_ci(xs, ys, conf=0.95):
    r = np.corrcoef(xs, ys)[0, 1]
    k = len(xs)
    z = np.arctanh(r)                       # Fisher transform to stabilise
    se = 1 / np.sqrt(k - 3)
    crit = stats.norm.ppf(1 - (1 - conf) / 2)
    return r, np.tanh(z - crit * se), np.tanh(z + crit * se)

for k in (20, 100, 5000):
    xs = rng.normal(0, 1, k)
    ys = 0.3 * xs + rng.normal(0, 1, k)
    r, lo, hi = pearson_ci(xs, ys)
    print(f"\nn={k:>5}  r={r:+.3f}  95% CI [{lo:+.3f}, {hi:+.3f}]"
          f"  width {hi-lo:.3f}")
# At n=20 the interval spans almost everything. The point estimate alone
# would have looked like a finding.

# --- SIGNIFICANCE IS NOT STRENGTH ---
big = 200_000
xb = rng.normal(0, 1, big)
yb = 0.008 * xb + rng.normal(0, 1, big)
res = stats.pearsonr(xb, yb)
print(f"\nn={big}: r={res.statistic:+.4f}, p={res.pvalue:.2e}")
print(f"r-squared = {res.statistic**2:.6f} -> explains "
      f"{res.statistic**2*100:.4f}% of the variance. Significant, not useful.")

# --- RANGE RESTRICTION: selection attenuates correlation. ---
full_x = rng.normal(0, 1, 5000)
full_y = 0.7 * full_x + rng.normal(0, 0.7, 5000)
print(f"\nfull population r = {np.corrcoef(full_x, full_y)[0,1]:+.3f}")
for pct in (50, 20, 5):
    cut = np.percentile(full_x, 100 - pct)
    keep = full_x >= cut
    print(f"  top {pct:>2}% on x only: r = "
          f"{np.corrcoef(full_x[keep], full_y[keep])[0,1]:+.3f}")
print("the relationship did not change; the observed correlation did.")

# --- MULTIPLE COMPARISONS IN A CORRELATION MATRIX ---
p_vars = 50
noise = pd.DataFrame(rng.normal(size=(200, p_vars)))   # NOTHING is related
corr = noise.corr().values
iu = np.triu_indices(p_vars, k=1)
vals = np.abs(corr[iu])
n_pairs = len(vals)
print(f"\n{p_vars} pure-noise variables -> {n_pairs} pairwise correlations")
print(f"largest |r| by chance: {vals.max():.3f}")
print(f"pairs with |r| > 0.15: {(vals > 0.15).sum()}")
# Bonferroni: the threshold a single test would use, divided by the count.
alpha = 0.05
print(f"uncorrected alpha {alpha} -> expect ~{int(n_pairs*alpha)} false hits")
print(f"Bonferroni threshold: {alpha/n_pairs:.2e}")

# --- CATEGORICAL ASSOCIATION: Cramer's V ---
def cramers_v(a: pd.Series, b: pd.Series) -> float:
    tab = pd.crosstab(a, b)
    chi2 = stats.chi2_contingency(tab).statistic
    nn = tab.values.sum()
    return float(np.sqrt(chi2 / (nn * (min(tab.shape) - 1))))

ca = pd.Series(rng.choice(list("xyz"), 900))
cb = pd.Series(np.where(rng.random(900) < 0.75, ca, rng.choice(list("xyz"), 900)))
print(f"\nCramer's V, related pair:   {cramers_v(ca, cb):.3f}")
print(f"Cramer's V, unrelated pair: "
      f"{cramers_v(ca, pd.Series(rng.choice(list('xyz'), 900))):.3f}")

# --- A CORRELATION SCREEN THAT DOES NOT MISLEAD ---
def screen(df: pd.DataFrame) -> pd.DataFrame:
    cols, out = df.select_dtypes("number").columns, []
    for i, a in enumerate(cols):
        for b in cols[i+1:]:
            s = df[[a, b]].dropna()
            if len(s) < 30:
                continue
            r, lo, hi = pearson_ci(s[a].values, s[b].values)
            rho = stats.spearmanr(s[a], s[b]).statistic
            out.append({"a": a, "b": b, "n": len(s), "pearson": round(r, 3),
                        "ci_lo": round(lo, 3), "ci_hi": round(hi, 3),
                        "spearman": round(rho, 3),
                        # A large gap flags a non-linear or outlier-driven pair.
                        "nonlinear_flag": abs(rho) - abs(r) > 0.15})
    return pd.DataFrame(out).sort_values("pearson", key=abs, ascending=False)

demo = pd.DataFrame({"x": x, **shapes})
print("\n", screen(demo).head(6).to_string(index=False))

  WHICH COEFFICIENT MEASURES WHAT

  measure        detects          robust to    range      needs
  ---------------------------------------------------------------------------
  Pearson r      linear only      nothing      [-1, 1]    continuous, linear
  Spearman rho   any monotonic    outliers     [-1, 1]    ordinal is fine
  Kendall tau    any monotonic    outliers,    [-1, 1]    small n ok
                                  ties
  mutual info    ANY dependence   outliers     [0, inf)   more data
  Cramer's V     categorical      -            [0, 1]     contingency table

  THE SHAPES AND WHAT EACH COEFFICIENT SAYS

  linear          monotone curve     u-shape          sinusoid
      .              .                 .     .         .   .
    .              .                    . .           . . . .
  .              .                        .             .   .
  r  +0.87        r  +0.62            r  ~0.00       r  ~0.00
  rho +0.86       rho +0.95           rho ~0.00      rho ~0.00
  MI  0.55        MI  0.72            MI  0.61       MI  0.84
                  ^ rho > r says      ^ only MI sees these two at all
                    "curved"

  FIVE EXPLANATIONS FOR ANY CORRELATION

    A -> B          the causal story you hoped for
    B -> A          reverse causation
    C -> A and B    confounding (the common one)
    selection       the sample was chosen in a way that links them
    chance          especially in a large correlation matrix

  RANGE RESTRICTION

    full range          top 20% only
    r = +0.71           r = +0.32
    . . . . . . . .                . . .
    . . . . . . . .    ------->    . . .   the slice is nearly round,
    . . . . . . . .                . . .   so the visible tilt is small

Key Takeaways

• Pearson measures linear association only, Spearman measures any monotonic association and resists outliers, and mutual information detects any dependence.
• A Spearman coefficient much larger than Pearson's is a reliable flag that the relationship is curved or outlier-driven.
• Report a confidence interval with every coefficient, since a small sample produces an interval spanning nearly every value.
• Statistical significance is not strength; at large sample sizes a meaningless correlation is highly significant.
• A correlation matrix over many variables produces notable-looking coefficients by chance, so treat it as a hypothesis generator and correct for multiplicity.

🧪 Practice

1. Compute Pearson, Spearman, and mutual information for a linear pair, a monotone curve, and a U-shape, and state which coefficient each pair needs.
2. Add a single leverage point to two unrelated variables and report how far Pearson moves and how far Spearman does.
3. Generate fifty unrelated columns, compute the full correlation matrix, and report the largest absolute coefficient found by chance alone.
4. Take a correlated pair, keep only the top 20% by one variable, and show the correlation shrink without the underlying relationship changing.
5. Interview: A colleague reports that hours of training correlate with performance at r equals 0.4, p below 0.001, in a sample of 12,000, and recommends mandatory training. What are your two main objections? (Hint: one is about what r-squared implies, the other about the five explanations for any correlation.)

Hypothesis Generation from Data

Exploration produces candidate explanations. Testing produces conclusions. Keeping those two activities separate is the most important methodological discipline in this chapter, and the reason is statistical rather than stylistic.

When you explore a dataset you implicitly conduct hundreds of tests. You look at a dozen segments, several time windows, and many variable pairs, and you stop when something looks interesting. That stopping rule is the problem: you have selected the most extreme pattern from many, and the most extreme pattern in pure noise still looks striking. Calculating a p-value for it afterwards is invalid, because the p-value assumes you specified the test before seeing the data.

The names for this are worth knowing because you will hear them in review. HARKing is hypothesising after the results are known, presenting a post-hoc pattern as a prediction. p-hacking is trying analyses until one reaches significance. The garden of forking paths is the subtler version: no dishonesty, just a long sequence of reasonable choices each nudged by the data, arriving at a result that would not replicate.

The standard example is the Texas sharpshooter: fire at a barn wall, then draw the target around the tightest cluster of holes and claim accuracy. Nothing was falsified. The target was just chosen afterwards.

So what do you actually do? You separate the modes explicitly.

In exploration mode, generate freely and record everything. Look at segments, plot distributions, follow hunches. Write down each pattern as a candidate hypothesis with the evidence and the effect size. Do not compute p-values, and do not phrase anything as a conclusion.

In confirmation mode, state the hypothesis and the test before touching the confirmation data. Then test on data the exploration never saw: a held-out split, a later time period, a different region, or ideally an experiment.

The mechanism that makes this practical is a holdout. Split the data, or better, split it by time, and explore only the first part. This costs you some sample size and buys you the ability to make a claim that survives. For work that matters, pre-registering the analysis plan, even informally in a document, is a strong discipline.

A good hypothesis from exploration has four parts: a mechanism that could plausibly produce the pattern, a direction so it is falsifiable, a magnitude you expect, and a test that would distinguish it from the obvious alternatives. A pattern without a mechanism is usually noise or a data artefact, and the first alternatives to rule out are always the boring ones: a pipeline change, a logging bug, a definition change, or a shift in population mix.

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

rng = np.random.default_rng(42)

# --- THE PROBLEM: exploration finds "findings" in pure noise. ---
n = 4000
noise = pd.DataFrame({
    "outcome": rng.normal(0, 1, n),                 # unrelated to everything
    "country": rng.choice(["US", "CA", "DE", "FR", "JP", "BR"], n),
    "device":  rng.choice(["ios", "android", "web"], n),
    "plan":    rng.choice(["free", "pro"], n),
    "weekday": rng.choice(range(7), n),
})

# A realistic exploration: slice by every segment and look for the biggest gap.
findings = []
for dim in ["country", "device", "plan", "weekday"]:
    g = noise.groupby(dim)["outcome"]
    means = g.mean()
    hi, lo = means.idxmax(), means.idxmin()
    a = noise.loc[noise[dim] == hi, "outcome"]
    b = noise.loc[noise[dim] == lo, "outcome"]
    t = stats.ttest_ind(a, b, equal_var=False)
    findings.append({"dimension": dim, "best": hi, "worst": lo,
                     "gap": round(means.max() - means.min(), 3),
                     "p_value": round(t.pvalue, 4)})
res = pd.DataFrame(findings).sort_values("p_value")
print("'findings' from data with NO real signal:\n", res.to_string(index=False))
print("\nThe outcome is pure noise. Every gap above is a chance fluctuation,")
print("and the p-values are invalid because the comparison was CHOSEN after")
print("looking at the means.")

# --- HOW MANY COMPARISONS DID THAT ACTUALLY RUN? ---
combos = sum(noise[d].nunique() * (noise[d].nunique() - 1) // 2
             for d in ["country", "device", "plan", "weekday"])
print(f"\nimplicit pairwise comparisons across four dimensions: {combos}")
print(f"expected 'significant at 0.05' results by chance: {combos * 0.05:.1f}")
print(f"probability at least one pair hits p<0.05: "
      f"{1 - 0.95 ** combos:.1%}")

# --- THE GARDEN OF FORKING PATHS, quantified. ---
paths = 0
for outlier_rule in ["keep", "drop_3sd"]:
    for transform in ["raw", "log"]:
        for window in ["all", "last_90d"]:
            for adjust in ["none", "by_country"]:
                paths += 1
print(f"\nfour innocuous analysis choices -> {paths} distinct analyses")
print("Each is defensible. Reporting only the one that 'worked' is p-hacking")
print("even when no single step felt dishonest.")

# --- THE FIX: SPLIT FIRST, EXPLORE ON ONE HALF ONLY. ---
# A dataset with ONE real effect (device) plus many noise dimensions.
m = 8000
real = pd.DataFrame({
    "country": rng.choice(["US", "CA", "DE", "FR"], m),
    "device":  rng.choice(["ios", "android", "web"], m),
    "plan":    rng.choice(["free", "pro"], m),
})
true_lift = {"ios": 0.30, "android": 0.0, "web": 0.0}   # the only real signal
real["outcome"] = [rng.normal(true_lift[d], 1) for d in real["device"]]

# Time-based split is preferred; here a random one for brevity.
explore = real.iloc[:4000]
confirm = real.iloc[4000:]

print("\n--- EXPLORATION HALF: generate candidates, no p-values ---")
candidates = []
for dim in ["country", "device", "plan"]:
    means = explore.groupby(dim)["outcome"].agg(["mean", "size"])
    gap = means["mean"].max() - means["mean"].min()
    candidates.append({"dimension": dim, "gap": round(gap, 3),
                       "top": means["mean"].idxmax()})
cand = pd.DataFrame(candidates).sort_values("gap", ascending=False)
print(cand.to_string(index=False))
print("Recorded as CANDIDATES. No conclusions drawn yet.")

print("\n--- CONFIRMATION HALF: test only the pre-stated hypotheses ---")
# State the hypotheses BEFORE touching `confirm`, then test exactly those.
HYPOTHESES = [("device", "ios"), ("country", cand.iloc[0]["top"])]
for dim, level in HYPOTHESES:
    a = confirm.loc[confirm[dim] == level, "outcome"]
    b = confirm.loc[confirm[dim] != level, "outcome"]
    t = stats.ttest_ind(a, b, equal_var=False)
    d = (a.mean() - b.mean()) / np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
    ci = stats.t.interval(0.95, len(a)+len(b)-2,
                          loc=a.mean()-b.mean(),
                          scale=np.sqrt(a.var(ddof=1)/len(a) + b.var(ddof=1)/len(b)))
    print(f"  {dim}={level:<8} diff {a.mean()-b.mean():+.3f}"
          f"  95% CI [{ci[0]:+.3f}, {ci[1]:+.3f}]"
          f"  d {d:+.2f}  p {t.pvalue:.4f}"
          f"  {'REPLICATES' if t.pvalue < 0.05 else 'does not replicate'}")
# Only the real effect survives. The chance pattern from exploration does not.

# --- STRUCTURING A CANDIDATE HYPOTHESIS ---
print("""
  A USABLE HYPOTHESIS HAS FOUR PARTS

  pattern     "iOS users convert 4 points higher than Android users"
  mechanism   "iOS checkout skips the address form via autofill"
  direction   "removing the form should raise Android conversion"
  magnitude   "we expect 2-4 points, not 0.2 and not 20"
  test        "ship autofill to 50% of Android traffic for two weeks"

  BORING ALTERNATIVES TO RULE OUT FIRST

    a logging or tagging change    a definition change in the metric
    a pipeline backfill           a shift in traffic mix / population
    a bot or internal-traffic      seasonality or a known campaign
      filter that changed
  Check these before any interesting explanation. They are more likely.""")

  THE TWO MODES, AND WHY THEY MUST NOT MIX

  EXPLORATION                        CONFIRMATION
  ------------------------------     ------------------------------
  goal: generate candidates          goal: test a stated claim
  many comparisons, informal         one pre-stated comparison
  p-values MEANINGLESS here          p-values valid here
  output: a ranked hypothesis list   output: an estimate + interval
  data: the exploration split        data: untouched holdout / later
                                           period / an experiment

  THE TEXAS SHARPSHOOTER

  fire first...              ...then draw the target
    .  . .   .                    .  . .   .
     . ..  .                       .(..). .     <- 'a tight cluster!'
   .   .  . .                    .   .  . .
  Nothing was falsified. The target was chosen after the shots.

  WHY THE P-VALUE IS INVALID AFTER EXPLORATION

    4 dimensions, 27 implicit pairwise comparisons
    P(at least one p < 0.05 | no real effect) = 1 - 0.95^27 = 75%

    So a "significant" segment gap found by browsing is the EXPECTED
    outcome under the null, not evidence against it.

  THE MINIMUM VIABLE DISCIPLINE

    1. split by time before you start
    2. explore only the earlier part
    3. write candidates in a document, with effect sizes, before testing
    4. test on the later part, or run an experiment
    5. report what failed to replicate as well as what did

Key Takeaways

• Exploration implicitly runs hundreds of comparisons, so a p-value computed on a pattern you found by browsing is invalid.
• The garden of forking paths produces unreplicable results through a chain of individually defensible choices, with no dishonesty required.
• Separate the modes: generate candidates on an exploration split, then test pre-stated hypotheses on a holdout, a later period, or an experiment.
• A usable hypothesis states a mechanism, a direction, an expected magnitude, and a test that distinguishes it from alternatives.
• Rule out the boring explanations first, since a logging change or a shift in traffic mix is more likely than an interesting effect.

🧪 Practice

1. Generate an outcome unrelated to four categorical dimensions, search for the largest segment gap, and report its p-value to show how easily noise produces a finding.
2. Count the implicit comparisons in that search and compute the probability that at least one reaches significance by chance.
3. Split a dataset containing exactly one real effect, generate candidates on the first half, and test them on the second, showing which replicate.
4. Write a candidate hypothesis with all four parts for a pattern in any dataset you have, plus three boring alternatives to rule out.
5. Interview: You find that users who use feature X retain at twice the rate of users who do not, and product wants to promote X heavily. What is your concern and what would you propose? (Hint: who chooses to use a feature, and what would distinguish the effect of X from the kind of user who adopts it?)

Segment and Cohort Analysis

An aggregate metric is an average over a population, and an average can be stable while everything underneath it changes. Segmentation and cohort analysis are the two standard ways of looking underneath, and they answer different questions.

Segmentation splits users by an attribute: plan, country, device, acquisition channel, company size. It answers "who behaves differently right now?" It is the tool for finding where a problem lives, and its main risk is that slicing enough ways will always produce an extreme slice, which is the multiple comparisons problem from the previous topic.

Cohort analysis groups users by when they started and follows each group over time. It answers "is the product getting better?", and it is the tool that separates two things an aggregate metric confuses: changes in behaviour, and changes in the mix of users.

That distinction is the whole reason cohort analysis exists, and it is worth spelling out. Suppose total retention is flat month over month. Two very different worlds produce that number. In one, every cohort retains identically, and nothing is changing. In the other, newer cohorts retain much better while the shrinking older cohorts retain worse, and the two effects cancel. The aggregate is identical; the businesses are not. Only a cohort view distinguishes them.

The analogy is a school's average test score. If it is flat year over year, you cannot tell whether teaching is unchanged, or whether teaching improved while the intake got weaker. Following each entering class separately answers it.

The standard construction is a cohort table: rows are the cohort, usually the signup month, columns are the periods since signup, and cells hold the retained fraction. Reading it has two directions and both matter. Reading across a row shows how one cohort decays, and the shape of that decay is a property of the product. Reading down a column compares cohorts at the same age, which is the only fair comparison and the one that shows whether things are improving.

Three practical cautions.

The triangle is not missing data. Recent cohorts have not existed long enough to have a period-12 number. The empty cells are unmeasured, not zero, and averaging down a column while ignoring that the cell counts differ produces a number that drifts for purely structural reasons.

Small cohorts are noisy. A cohort of forty users produces retention percentages that jump around by several points for no reason. Show the counts next to the rates, or your reader will interpret noise as trend.

Define the event precisely. Retention depends entirely on what counts as active, and whether periods are calendar months or rolling 30-day windows from each user's own signup date. Those two choices produce materially different tables, and mixing them across analyses is a common source of numbers that do not reconcile.

Beyond retention, the same structure supports revenue cohorts, which reveal whether newer customers are worth more, and survival analysis, which handles the censoring properly when you care about time-to-event rather than a fixed checkpoint.

import numpy as np
import pandas as pd

rng = np.random.default_rng(5)

# --- BUILD A SYNTHETIC EVENT LOG WITH IMPROVING COHORTS ---
signup_months = pd.period_range("2024-01", "2024-08", freq="M")
records = []
for i, m in enumerate(signup_months):
    n_users = [500, 520, 600, 640, 700, 900, 1100, 1200][i]
    # Later cohorts retain better: the product is genuinely improving.
    base_retention = 0.55 + 0.025 * i
    for u in range(n_users):
        uid = f"{m}-{u}"
        records.append({"user": uid, "cohort": m, "period": 0})
        alive = True
        for p in range(1, 9):
            if m + p > signup_months[-1]:       # the future has not happened
                break
            # Monthly survival decays with age, improves with cohort.
            p_survive = base_retention * np.exp(-0.12 * (p - 1))
            alive = alive and (rng.random() < p_survive)
            if alive:
                records.append({"user": uid, "cohort": m, "period": p})
log = pd.DataFrame(records)

# --- THE COHORT TABLE ---
sizes = log[log["period"] == 0].groupby("cohort")["user"].nunique()
counts = log.pivot_table(index="cohort", columns="period",
                         values="user", aggfunc="nunique")
rates = counts.div(sizes, axis=0) * 100

print("cohort sizes:\n", sizes.to_string())
print("\nretention %, rows = cohort, columns = months since signup:")
print(rates.round(1).to_string())
print("\nNaN cells are UNMEASURED, not zero: those months have not happened.")

# --- READ IT TWO WAYS ---
print("\nACROSS a row (cohort 2024-01 decaying over its own lifetime):")
print(rates.loc[signup_months[0]].dropna().round(1).to_string())

print("\nDOWN column 3 (every cohort at the same age -- the FAIR comparison):")
col = pd.DataFrame({"retention_m3": rates[3].round(1), "cohort_size": sizes})
print(col.dropna().to_string())
print("month-3 retention is rising across cohorts: the product is improving.")

# --- WHY THE AGGREGATE HIDES THIS ---
# The blended 'active users / total users' number mixes cohorts of every age.
active_by_month = {}
for m in signup_months:
    total = alive = 0
    for c in signup_months:
        if c > m:
            continue
        age = (m - c).n
        if age in counts.columns and not pd.isna(counts.loc[c, age]):
            alive += counts.loc[c, age]
            total += sizes[c]
    active_by_month[str(m)] = 100 * alive / total
blended = pd.Series(active_by_month).round(1)
print("\nblended retention across all cohorts:\n", blended.to_string())
print("\nThis single series moves for two reasons at once: real changes in")
print("behaviour, and a changing mix of cohort ages. It cannot separate them.")

# --- THE CANCELLATION CASE: flat aggregate, two opposite trends. ---
demo = pd.DataFrame({
    "cohort": ["old"] * 3 + ["new"] * 3,
    "period": [1, 2, 3] * 2,
    "retention": [60, 50, 42, 75, 66, 58],       # new cohorts clearly better
    "users": [8000, 8000, 8000, 1000, 1000, 1000],  # but old dominates volume
})
weighted = (demo.assign(w=demo["retention"] * demo["users"])
                .groupby("period").apply(
                    lambda g: g["w"].sum() / g["users"].sum(),
                    include_groups=False))
print("\nper-cohort retention:\n",
      demo.pivot(index="cohort", columns="period",
                 values="retention").to_string())
print("\nsize-weighted blend:\n", weighted.round(1).to_string())
print("The blend tracks the large old cohort and understates the improvement.")

# --- SEGMENTATION: where does a problem live, and is the slice trustworthy? ---
users = pd.DataFrame({
    "plan": rng.choice(["free", "pro", "enterprise"], 3000, p=[.7, .25, .05]),
    "device": rng.choice(["ios", "android", "web"], 3000),
})
p_churn = {"free": .30, "pro": .12, "enterprise": .05}
users["churned"] = [rng.random() < p_churn[p] for p in users["plan"]]

def segment_report(df, dim, outcome, min_n=100):
    g = df.groupby(dim)[outcome].agg(n="size", rate="mean")
    overall = df[outcome].mean()
    g["rate"] = (g["rate"] * 100).round(1)
    g["lift_vs_overall"] = (g["rate"] - overall * 100).round(1)
    # A standard error on a proportion, so the reader sees the noise floor.
    g["std_err_pts"] = (100 * np.sqrt(
        (g["rate"]/100) * (1 - g["rate"]/100) / g["n"])).round(2)
    g["trustworthy"] = g["n"] >= min_n
    return g.sort_values("rate", ascending=False)

print("\nchurn by plan:\n", segment_report(users, "plan", "churned").to_string())
print("\nchurn by device (no real effect -- note the gaps sit within ~2 SE):")
print(segment_report(users, "device", "churned").to_string())
# Printing the standard error next to the rate is what stops a reader from
# reading a 1.5-point noise gap as a finding.

# --- CROSS-SEGMENT: where slicing gets dangerous. ---
cross = (users.groupby(["plan", "device"])["churned"]
              .agg(n="size", rate="mean"))
cross["rate"] = (cross["rate"] * 100).round(1)
print(f"\n{len(cross)} cells from two dimensions;"
      f" {(cross['n'] < 100).sum()} of them have n < 100")
print(cross.to_string())
print("\nThe most extreme cell in a grid of small cells is almost always")
print("noise. Require a minimum n before reporting a cell at all.")

  THE COHORT TABLE AND ITS TWO READINGS

                   months since signup
  cohort    n     0     1     2     3     4     5
  2024-01  500  100   55.2  48.1  42.0  37.4  33.1   --> READ ACROSS:
  2024-02  520  100   57.9  50.4  44.6  39.8   .          one cohort's decay
  2024-03  600  100   60.1  53.0  47.2   .     .          (a product property)
  2024-04  640  100   62.8  55.6   .    .     .
  2024-05  700  100   65.0   .    .     .     .
  2024-06  900  100    .     .    .     .     .
                       |
                       v
                  READ DOWN: every cohort at the SAME age.
                  The only fair comparison, and the one that
                  answers "are we getting better?"

  The empty triangle is UNMEASURED time, not churn and not zero.

  WHAT THE BLENDED NUMBER CANNOT TELL YOU

  aggregate retention: 52%, 52%, 52%   <- flat, three months running

  world A                          world B
  every cohort retains 52%         new cohorts 75%, old cohorts 40%,
  nothing is changing              and the mix shifts to cancel out
                                   the product is improving fast
  Identical aggregate. Only the cohort table separates them.

  SEGMENT VS COHORT

             segmentation                cohort analysis
  splits by  an attribute (plan, geo)    a start date
  answers    where is the problem now?   are we improving over time?
  risk       multiple comparisons        misreading the empty triangle
  guardrail  minimum n + standard error  show cohort sizes beside rates

Key Takeaways

• Segmentation finds where a problem lives right now; cohort analysis shows whether the product is improving over time.
• A flat aggregate metric is consistent with both "nothing is changing" and "behaviour improved while the mix shifted", and only a cohort view separates them.
• Read a cohort table across a row for one cohort's decay and down a column to compare cohorts at the same age, which is the only fair comparison.
• The empty triangle is unmeasured time rather than zero, so never average a column without accounting for which cohorts are present.
• Show cohort sizes and a standard error beside every rate, because the most extreme small segment is almost always noise.

🧪 Practice

1. Build a cohort retention table from an event log with signup dates, and print the cohort sizes next to it.
2. Construct a case where per-cohort retention improves every month while the size-weighted blended metric stays flat or falls.
3. Write a segment report that outputs the rate, the row count, and the standard error per level, and use it to argue that one apparent gap is noise.
4. Cross two segmentation dimensions and count how many resulting cells fall below a minimum sample size you consider reportable.
5. Interview: Monthly active users are flat for six months and leadership concludes the product has stalled. Which two views would you build first, and what could each reveal? (Hint: one separates new from returning users, the other follows each signup month separately.)

<a id="72-visualization-principles"></a>

7.2 Visualization Principles

A chart is a claim about data that the reader will believe, which makes chart design a correctness concern rather than a matter of taste; this section covers the grammar that generates charts, the rules for choosing among them, and the specific constructions that mislead.

Grammar of Graphics

Most people learn charts as a catalogue: here is a bar chart, here is a scatter plot, here is a pie chart, pick one. That model breaks down as soon as you need something the catalogue does not contain, and it gives you no way to reason about why one chart works and another does not.

The grammar of graphics replaces the catalogue with a small set of components that combine. Leland Wilkinson formalised it in 1999 and Hadley Wickham's ggplot2 made it the dominant mental model; Vega-Lite, plotnine, and Altair all implement it. The insight is that a chart is not a type, it is a mapping from data columns to visual properties, plus a choice of geometry.

The components:

Data is the tidy table you start from, one row per observation. The grammar assumes tidy input, which is why chapter 6's tidy-data discipline pays off here.

Aesthetic mappings connect columns to visual properties: x position, y position, colour, size, shape, opacity. This is the heart of the grammar. You do not "colour the bars blue", you map the region column to the colour channel, and the library derives the legend, the scale, and the assignment.

Geometries are the marks that represent each row: point, line, bar, area, text, tile. Changing only the geometry turns a scatter plot into a line chart, which reveals that these are not different kinds of chart so much as different marks over the same mapping.

Statistical transformations are computations the chart performs before drawing: binning for a histogram, counting for a bar chart of frequencies, smoothing for a trend line, computing quartiles for a box plot. Recognising that a histogram is "bar geometry plus a binning stat" is what makes the grammar generative.

Scales control how data values become visual values: the axis range, whether it is linear or logarithmic, which colours a categorical variable receives.

Coordinate systems position the marks. Cartesian is the default; flipping it gives horizontal bars; polar coordinates turn a stacked bar into a pie chart, which is a genuinely useful thing to understand, because it explains why pie charts are hard to read: they encode with angle instead of length.

Facets split the data into panels by a variable, producing small multiples.

The payoff is composability. A histogram, a density plot, a box plot, and a violin plot are all the same data with different stats and geometries. Once you think this way you stop looking for the right chart type in a list and start asking which channels your variables should map to.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")                     # headless: render to file, not screen
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
df = pd.DataFrame({
    "month": pd.date_range("2024-01-01", periods=24, freq="MS").tolist() * 2,
    "region": ["north"] * 24 + ["south"] * 24,
})
df["revenue"] = (100 + np.arange(48) % 24 * 4
                 + np.where(df["region"] == "north", 30, 0)
                 + rng.normal(0, 8, 48)).round(1)

# --- THE SAME MAPPING, THREE GEOMETRIES ---
# Only the mark changes. The data, the aesthetics, and the scales are fixed.
fig, axes = plt.subplots(1, 3, figsize=(13, 3.2), sharey=True)
north = df[df["region"] == "north"]
for ax, geom in zip(axes, ["point", "line", "bar"]):
    if geom == "point":
        ax.scatter(north["month"], north["revenue"], s=18, color="#2a78d6")
    elif geom == "line":
        ax.plot(north["month"], north["revenue"], lw=2, color="#2a78d6")
    else:
        ax.bar(north["month"], north["revenue"], width=20, color="#2a78d6")
    ax.set_title(f"geom = {geom}", loc="left", fontsize=10)
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/geoms.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- AESTHETIC MAPPING: colour is a MAPPING, not a decoration. ---
# Wrong mental model: "make the north line blue".
# Right mental model: map `region` -> colour; the legend follows automatically.
PALETTE = {"north": "#2a78d6", "south": "#eb6834"}    # fixed, per entity

fig, ax = plt.subplots(figsize=(7, 3.4))
for region, g in df.groupby("region"):
    ax.plot(g["month"], g["revenue"], lw=2, color=PALETTE[region], label=region)
    # Direct-label the endpoint so identity is not colour-alone.
    last = g.iloc[-1]
    ax.annotate(region, (last["month"], last["revenue"]),
                xytext=(6, 0), textcoords="offset points",
                color=PALETTE[region], fontsize=9, va="center")
ax.legend(frameon=False, loc="upper left", fontsize=9)
ax.set_ylabel("revenue")
ax.spines[["top", "right"]].set_visible(False)
ax.grid(axis="y", lw=0.5, color="#e6e6e4")            # hairline, solid, recessive
ax.set_axisbelow(True)
fig.savefig("/tmp/claude-0/mapping.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- STATISTICAL TRANSFORMATIONS: four charts, one column. ---
# Each is "a stat + a geom". This is what makes the grammar generative.
values = rng.lognormal(3.2, 0.5, 900)
recipes = {
    "histogram":   "stat=bin,      geom=bar",
    "density":     "stat=kde,      geom=area",
    "box plot":    "stat=quantile, geom=box",
    "violin":      "stat=kde,      geom=violin (mirrored)",
    "ECDF":        "stat=ecdf,     geom=step",
}
for name, recipe in recipes.items():
    print(f"  {name:11s} = {recipe}")

fig, axes = plt.subplots(1, 4, figsize=(13, 3))
axes[0].hist(values, bins=30, color="#2a78d6")
axes[0].set_title("bin + bar", loc="left", fontsize=10)
from scipy.stats import gaussian_kde
grid = np.linspace(values.min(), values.max(), 200)
axes[1].fill_between(grid, gaussian_kde(values)(grid), color="#2a78d6", alpha=0.25)
axes[1].plot(grid, gaussian_kde(values)(grid), color="#2a78d6", lw=2)
axes[1].set_title("kde + area", loc="left", fontsize=10)
axes[2].boxplot(values, vert=True, widths=0.4)
axes[2].set_title("quantile + box", loc="left", fontsize=10)
xs = np.sort(values)
axes[3].step(xs, np.arange(1, len(xs) + 1) / len(xs), color="#2a78d6", lw=2)
axes[3].set_title("ecdf + step", loc="left", fontsize=10)
for ax in axes:
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/stats.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- COORDINATE SYSTEMS: a pie chart IS a stacked bar in polar coordinates. ---
parts = pd.Series({"a": 40, "b": 30, "c": 18, "d": 12})
fig, axes = plt.subplots(1, 2, figsize=(9, 3.2))
left = 0
for i, (k, v) in enumerate(parts.items()):
    axes[0].barh([0], [v], left=left, height=0.4,
                 color=["#2a78d6", "#eb6834", "#1baf7a", "#eda100"][i],
                 edgecolor="white", linewidth=2)   # 2px surface gap, not a border
    left += v
axes[0].set_title("cartesian: stacked bar (length)", loc="left", fontsize=10)
axes[0].set_yticks([])
axes[1].pie(parts.values, labels=parts.index,
            colors=["#2a78d6", "#eb6834", "#1baf7a", "#eda100"],
            wedgeprops={"edgecolor": "white", "linewidth": 2})
axes[1].set_title("polar: the same bar (angle)", loc="left", fontsize=10)
fig.savefig("/tmp/claude-0/coords.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("\nSame data, same stat, same geom -- only the coordinate system differs.")
print("The pie encodes with ANGLE, which people judge far worse than LENGTH,")
print("which is the whole reason to prefer the bar.")

# --- FACETS: split into panels by a variable. ---
fig, axes = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
for ax, (region, g) in zip(axes, df.groupby("region")):
    ax.plot(g["month"], g["revenue"], lw=2, color=PALETTE[region])
    ax.set_title(region, loc="left", fontsize=10)
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/facets.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- THE SAME CHART, DECLARED GRAMMATICALLY (plotnine / ggplot2 syntax) ---
print("""
  ggplot2 / plotnine, where the grammar is explicit:

    (ggplot(df, aes(x="month", y="revenue", color="region"))
     + geom_line(size=1)
     + geom_smooth(method="lm", se=False)    # a second layer, same mapping
     + scale_y_continuous(limits=(0, None))  # a bar-safe zero baseline
     + facet_wrap("~region")                 # small multiples
     + labs(title="Revenue by region", y="Revenue (USD)"))

  Read it as: data + mapping + geometry + scale + facet. Adding a layer or
  swapping a geom is a one-line change, because nothing about the mapping
  was baked into the chart type.
""")

  THE SEVEN COMPONENTS

    DATA           a tidy table, one row per observation
      |
    AESTHETICS     column -> channel   (x, y, colour, size, shape, alpha)
      |
    STATISTICS     a computation before drawing (bin, count, smooth, quantile)
      |
    GEOMETRY       the mark (point, line, bar, area, tile, text)
      |
    SCALES         data value -> visual value (range, log, palette)
      |
    COORDINATES    cartesian, flipped, polar
      |
    FACETS         split into panels by a variable
      |
      v
    the chart

  ONE COLUMN, FIVE CHARTS -- ONLY stat AND geom CHANGE

    stat=bin      + geom=bar     -> histogram
    stat=kde      + geom=area    -> density plot
    stat=quantile + geom=box     -> box plot
    stat=kde      + geom=violin  -> violin plot
    stat=ecdf     + geom=step    -> ECDF

  WHY THE PIE CHART IS HARD TO READ

    stacked bar  = stat=identity, geom=bar, coord=cartesian  -> LENGTH
    pie chart    = stat=identity, geom=bar, coord=POLAR      -> ANGLE

    Same numbers. Length is judged accurately; angle is not.
    The grammar makes the pie's weakness a property of one component.

Key Takeaways

• A chart is a mapping from columns to visual channels plus a geometry, not an entry in a catalogue of chart types.
• Aesthetic mappings mean you map a column to colour rather than colouring a mark, which is what makes legends and scales derive automatically.
• A histogram, density plot, box plot, and violin plot are the same column under different statistical transformations and geometries.
• A pie chart is a stacked bar in polar coordinates, which is exactly why it encodes with angle and reads worse than a bar.
• The grammar assumes tidy input, so the reshaping discipline from the previous chapter is a prerequisite.

🧪 Practice

1. Take one numeric column and produce a histogram, a density plot, a box plot, and an ECDF, labelling each with its stat and geom.
2. Build the same line chart twice, once by hard-coding colours per series and once by mapping a column to colour, and describe what the second version gives you for free.
3. Render the same part-to-whole data as a stacked bar and as a pie, then ask someone to rank the segments from each and compare their accuracy.
4. Interview: A stakeholder asks for "a chart showing revenue by region over time, split by product line". Decompose the request into grammar components. (Hint: four variables, and you have x, y, colour, and facet available.)

Choosing the Right Chart Type

The most common charting mistake is not an ugly chart, it is a chart that answers a different question than the reader has. Choosing well starts from one question: what does the reader need to do with this data? The job determines the form, and the form determines the encoding.

Before picking a chart, ask whether you need one at all. A single current value is a stat tile, a number with a label and maybe a delta, and a one-bar bar chart is a worse version of it. A handful of headline numbers is a row of stat tiles. The number a report leads with is a hero figure, set large. More than about seven meaningful classes is usually a table, possibly beside a chart, because past that point adjacent colours blur and the chart stops functioning.

When a chart is right, map the job to the form:

The reader must	Default form	Why
Compare magnitudes	bar or column	length from a shared baseline is read precisely
See a trend over time	line	connection implies ordered continuity
See a distribution	histogram, ECDF, violin	shape is the subject
Compare distributions	box or violin per group	compact, sortable
See a relationship	scatter, hexbin if dense	two positions at once
See part of a whole	stacked bar	length, not angle
See an ordered share	diverging stacked bar	centred on the neutral point
Compare before and after	dumbbell or slope chart	two points, one connection
See a value on a grid	heatmap	many cells, sequential colour
See a value on a map	choropleth or symbol map	geography is the question
Notice one series	emphasis: highlight one, gray the rest	the story is one line

A few rules behind the table are worth stating, because they resolve most real decisions.

Bars need a zero baseline; lines do not. A bar encodes value as length from the axis, so truncating the axis lies about the ratio between bars. A line encodes change through slope, so a truncated axis is legitimate and often necessary to see the variation at all. This is the single most useful asymmetry in charting.

Never use a dual axis. Two y-scales on one plot invent a relationship, because the alignment of the two scales is arbitrary: slide one and the series appear to cross wherever you like. The honest alternatives are two charts stacked, small multiples, or indexing both series to 100 at a common starting point so they share one axis.

Emphasis is the most underused form. When the story is "this one segment moved", colouring eight series in eight hues buries it. One series in an accent colour with the rest in gray is usually the honest answer to "make this chart clearer".

Horizontal bars for long labels and many categories. Rotated text is hard to read and there is more vertical room on a page than horizontal. Sort bars by value unless the categories have an inherent order such as time or size tiers.

Pie charts are acceptable only for part-to-whole at a glance with few segments. For comparing close values, use a bar or just print the numbers.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

BLUE, ORANGE, GRAY, INK = "#2a78d6", "#eb6834", "#c9c9c5", "#0b0b0b"

# --- IS IT EVEN A CHART? A single value is a stat tile, not a bar. ---
fig, axes = plt.subplots(1, 2, figsize=(9, 2.6))
axes[0].bar(["conversion"], [3.4], width=0.3, color=BLUE)
axes[0].set_title("one-bar bar chart (wrong)", loc="left", fontsize=10)
axes[0].set_ylim(0, 5)
axes[1].axis("off")
axes[1].text(0, .62, "3.4%", fontsize=38, color=INK, va="center")   # hero figure
axes[1].text(0, .28, "Conversion rate", fontsize=10, color="#52514e")
axes[1].text(0, .10, "+0.5 pts vs last month", fontsize=9, color="#008300")
axes[1].set_title("stat tile (right)", loc="left", fontsize=10)
fig.savefig("/tmp/claude-0/stat_tile.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- THE BASELINE ASYMMETRY: bars must start at zero, lines need not. ---
cats = ["A", "B", "C", "D"]
vals = [102, 105, 100, 108]
months = pd.date_range("2024-01-01", periods=12, freq="MS")
series = 100 + np.cumsum(np.random.default_rng(1).normal(0.4, 1.2, 12))

fig, axes = plt.subplots(1, 4, figsize=(14, 3))
axes[0].bar(cats, vals, color=BLUE, width=0.55)
axes[0].set_ylim(0, 120)
axes[0].set_title("bar, zero baseline: HONEST", loc="left", fontsize=9)

axes[1].bar(cats, vals, color=ORANGE, width=0.55)
axes[1].set_ylim(98, 110)          # D now looks ~4x B. It is 3% larger.
axes[1].set_title("bar, truncated: LIES", loc="left", fontsize=9)

axes[2].plot(months, series, lw=2, color=BLUE)
axes[2].set_ylim(0, series.max() * 1.1)
axes[2].set_title("line, zero baseline: hides the trend", loc="left", fontsize=9)

axes[3].plot(months, series, lw=2, color=BLUE)
axes[3].set_title("line, fitted range: FINE", loc="left", fontsize=9)
for ax in axes:
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/baseline.png", dpi=110, bbox_inches="tight")
plt.close(fig)

print("bar D is 8% above bar C, and the truncated panel makes it look 4x.")
print("A bar encodes LENGTH from the axis; a line encodes SLOPE. Hence the")
print("asymmetry: zero is mandatory for bars and optional for lines.")

# --- NEVER A DUAL AXIS. Index to a common base instead. ---
users = np.array([12, 14, 15, 18, 21, 24, 26, 29, 31, 34, 36, 39]) * 1000
sessions = np.array([190, 250, 230, 300, 360, 420,
                     390, 470, 520, 560, 600, 640]) * 1000

fig, axes = plt.subplots(1, 3, figsize=(14, 3.2))
ax1 = axes[0]
ax2 = ax1.twinx()                                  # the anti-pattern
ax1.plot(months, users, color=BLUE, lw=2)
ax2.plot(months, sessions, color=ORANGE, lw=2)
ax1.set_title("dual axis: WRONG", loc="left", fontsize=9)
# Slide the right-hand scale and the crossing point moves anywhere you want.
ax2.set_ylim(0, sessions.max() * 2.4)

axes[1].plot(months, users / users[0] * 100, color=BLUE, lw=2, label="users")
axes[1].plot(months, sessions / sessions[0] * 100, color=ORANGE, lw=2,
             label="sessions")
axes[1].axhline(100, lw=0.8, color="#c9c9c5")
axes[1].set_title("indexed to 100: one axis, honest", loc="left", fontsize=9)
axes[1].legend(frameon=False, fontsize=9)

axes[2].plot(months, sessions / users, color=BLUE, lw=2)
axes[2].set_title("the ratio, if that IS the question", loc="left", fontsize=9)
for ax in axes:
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/dual_axis.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- EMPHASIS BEATS EIGHT HUES WHEN THE STORY IS ONE SERIES ---
rng = np.random.default_rng(4)
many = pd.DataFrame(
    {f"region_{i}": 100 + np.cumsum(rng.normal(0.2, 1.1, 12)) for i in range(7)},
    index=months)
many["region_hero"] = 100 + np.cumsum(rng.normal(1.6, 1.1, 12))   # the story

fig, axes = plt.subplots(1, 2, figsize=(11, 3.4), sharey=True)
for i, col in enumerate(many.columns):
    axes[0].plot(many.index, many[col], lw=2, label=col)
axes[0].set_title("8 categorical hues: the story is buried", loc="left", fontsize=9)
axes[0].legend(frameon=False, fontsize=6, ncol=2)

for col in many.columns:
    is_hero = col == "region_hero"
    axes[1].plot(many.index, many[col], lw=2 if is_hero else 1.4,
                 color=BLUE if is_hero else GRAY, zorder=3 if is_hero else 1)
axes[1].annotate("region_hero", (many.index[-1], many["region_hero"].iloc[-1]),
                 xytext=(-78, 6), textcoords="offset points",
                 color=BLUE, fontsize=9, fontweight="bold")
axes[1].set_title("emphasis: one accent, rest gray", loc="left", fontsize=9)
for ax in axes:
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/emphasis.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- HORIZONTAL AND SORTED for long labels and many categories ---
labels = ["Enterprise self-serve", "Mid-market direct", "SMB partner channel",
          "Consumer mobile app", "Education licensing", "Government contracts"]
amounts = [412, 388, 265, 190, 96, 74]
fig, axes = plt.subplots(1, 2, figsize=(12, 3.4))
axes[0].bar(labels, amounts, color=BLUE, width=0.6)
axes[0].tick_params(axis="x", rotation=45, labelsize=7)
axes[0].set_title("vertical, rotated labels: hard to read", loc="left", fontsize=9)
order = np.argsort(amounts)                            # ascending -> biggest on top
axes[1].barh([labels[i] for i in order], [amounts[i] for i in order],
             color=BLUE, height=0.6)
axes[1].tick_params(labelsize=8)
axes[1].set_title("horizontal and sorted: readable", loc="left", fontsize=9)
for ax in axes:
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/horizontal.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- A DECISION FUNCTION, to make the reasoning explicit ---
def suggest_form(n_values: int, x_type: str, n_series: int,
                 story: str = "compare") -> str:
    if n_values == 1:
        return "stat tile or hero figure (never a one-bar chart)"
    if story == "one_series_moved":
        return "emphasis: accent the one series, gray the rest"
    if x_type == "time":
        extra = " (or small multiples above 5 series)" if n_series > 5 else ""
        return "line" + extra
    if x_type == "numeric":
        dense = " -> hexbin or 2-D histogram if dense" if n_values > 5000 else ""
        return "scatter" + dense
    if x_type == "categorical":
        if n_series > 1:
            return "grouped bar (few series) or stacked bar (part-to-whole)"
        if n_values > 12:
            return "horizontal bar, sorted by value"
        return "bar, sorted by value, zero baseline"
    return "table"

for case in [(1, "none", 1, "compare"), (48, "time", 8, "one_series_moved"),
             (40000, "numeric", 1, "compare"), (25, "categorical", 1, "compare"),
             (12, "time", 3, "compare")]:
    print(f"  {str(case):32s} -> {suggest_form(*case)}")

  THE JOB DECIDES THE FORM

  job                        form                     encoding channel
  ---------------------------------------------------------------------------
  one current value          stat tile / hero number  text
  a few headline numbers     KPI row of stat tiles    text
  compare magnitudes         bar (horizontal if long) length from baseline
  trend over time            line                     position + slope
  one series is the point    emphasis (1 hue + gray)  position + salience
  distribution shape         histogram / ECDF         position + area
  compare distributions      box / violin per group   position
  relationship               scatter (hexbin if dense) two positions
  part-to-whole              stacked bar              length
  ordered share (Likert)     diverging stacked bar    length from centre
  value on a grid            heatmap                  sequential colour
  more than ~7 classes       a table                  text

  THE BASELINE RULE

  BAR: encodes LENGTH            LINE: encodes SLOPE
   108|    #  <- looks 4x B       truncating is FINE, and often required
   105| #  #                      to see the variation at all
   102|#   #
      +------- truncated at 98    zero-baselining a line chart flattens
      A B C D                     the trend into a useless straight line
  Truncating a bar axis lies about the RATIO. Zero is mandatory.

  DUAL AXIS: WHY IT IS ALWAYS WRONG

    users     0 -> 40k   on the left
    sessions  0 -> 800k  on the right
    The crossing point is set by YOUR choice of the two ranges, not by the
    data. Slide one scale and the story reverses. Fixes: two charts, small
    multiples, or index both to 100 at t0 on ONE axis.

  SERIES-COUNT LADDER (categorical colour)

    1-3    colour alone is comfortable; direct-label
    4      direct labels become mandatory
    5-6    soft cap: legend or small multiples
    7-8    ceiling; past it fold the tail into "Other" or facet
    9+     never generate another hue -- it is indistinguishable under CVD

Key Takeaways

• Ask what the reader must do with the data, and whether a stat tile, hero number, or table serves better than any chart.
• Bars encode length so they require a zero baseline; lines encode slope so a truncated axis is legitimate.
• Dual-axis charts invent a relationship because the alignment of the two scales is arbitrary; index both series to a common base instead.
• When one series is the story, emphasis with a single accent and gray context beats eight categorical hues.
• Use horizontal sorted bars for long labels and many categories, and stop adding hues past about eight series.

🧪 Practice

1. Render the same four-category data as a zero-baselined bar chart and a truncated one, and state the ratio a reader would infer from each.
2. Take two series of very different scale, build the dual-axis version, then rebuild it indexed to 100 and explain what the first version implied.
3. Convert an eight-series line chart into an emphasis chart and note how many colours you needed.
4. Write a function mapping value count, x-axis type, series count, and story to a recommended chart form, and test it on five real requests.
5. Interview: A manager asks you to "make the difference clearer" in a bar chart where the bars are 102 and 108. What do you do? (Hint: what is the honest way to show a 6% difference, and what would truncating the axis imply?)

Encoding Data with Visual Channels

A chart works by converting numbers into visual properties, and human vision reads some of those properties far more accurately than others. That ranking is not a style preference, it is a measured perceptual fact, and knowing it turns chart design into an ordering problem: put your most important variable on your most accurate channel.

The foundational work is Cleveland and McGill's 1984 experiments, which asked people to judge magnitudes from different encodings and measured the errors. The resulting ranking, from most to least accurate for quantitative data:

1. Position along a common scale. Two bars on the same axis, two points in the same scatter. This is by a wide margin the most accurate channel.
2. Position along identical non-aligned scales. The same variable in two panels of a small-multiple grid.
3. Length. A bar not sharing a baseline with the one you compare it to, such as segments inside different stacked bars.
4. Angle and slope. Pie wedges, line steepness.
5. Area. Bubble sizes, treemap rectangles.
6. Volume. Three-dimensional representations.
7. Colour saturation and lightness. Heatmap intensity, choropleth shading.
8. Colour hue. Which category something is.

Two consequences follow directly.

Use position for what matters most. If the reader's main question is "which is bigger", put that variable on an axis. This is why bars beat pies: the pie moves the comparison from position and length down to angle.

Area is judged poorly, and the error is systematic. People underestimate area differences. A circle with twice the radius has four times the area and is usually read as roughly twice as large. This means bubble charts communicate only coarse differences, and it means any chart sizing marks by value must scale by area, not radius, or it exaggerates by the square.

The channels also differ in type suitability, which is a separate question from accuracy. Hue is poor for magnitude and excellent for identity, because hues have no natural order: nobody thinks green is more than purple. Lightness has a natural order, which is exactly why sequential colour scales vary lightness rather than hue. Shape carries identity for a handful of categories and no magnitude at all.

A final rule that removes a surprising amount of clutter: do not encode the same variable twice. Colouring bars darker where they are taller adds no information, because bar length already showed it, and it spends the colour channel, which was the one channel still free to carry something else. Reserve each channel for one variable.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

BLUE, ORANGE, GRAY = "#2a78d6", "#eb6834", "#c9c9c5"
values = np.array([28.0, 34.0, 31.0, 40.0, 36.0])
labels = list("ABCDE")

# --- THE SAME FIVE NUMBERS ON FIVE CHANNELS, best to worst ---
fig, axes = plt.subplots(1, 5, figsize=(15, 2.9))

axes[0].bar(labels, values, color=BLUE, width=0.6)     # POSITION + LENGTH
axes[0].set_ylim(0, 45)
axes[0].set_title("1. position, common scale", loc="left", fontsize=9)

# LENGTH without a shared baseline: segments inside separate stacks.
offsets = np.array([0, 12, 5, 18, 9])
axes[1].bar(labels, values, bottom=offsets, color=BLUE, width=0.6)
axes[1].set_title("3. length, no shared baseline", loc="left", fontsize=9)

axes[2].pie(values, labels=labels, wedgeprops={"edgecolor": "white",
                                               "linewidth": 2})
axes[2].set_title("4. angle", loc="left", fontsize=9)

# AREA: the sizes must be proportional to VALUE, so radius ~ sqrt(value).
axes[3].scatter(range(5), [0] * 5, s=values * 40, color=BLUE)  # s IS area
axes[3].set_xticks(range(5), labels)
axes[3].set_yticks([])
axes[3].set_title("5. area", loc="left", fontsize=9)

axes[4].imshow(values.reshape(1, -1), cmap="Blues", aspect="auto")
axes[4].set_xticks(range(5), labels)
axes[4].set_yticks([])
axes[4].set_title("7. colour lightness", loc="left", fontsize=9)

for ax in axes[:2]:
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/channels.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("Rank A..E from each panel. Accuracy falls sharply left to right,")
print("and the underlying numbers never changed.")

# --- THE AREA TRAP: scaling radius by value exaggerates by the square. ---
a, b = 10.0, 20.0                       # b is exactly twice a
r_wrong = np.array([a, b])              # radius proportional to value
r_right = np.sqrt(np.array([a, b]))     # AREA proportional to value
print(f"\nvalue ratio                     {b/a:.1f}x")
print(f"area ratio, radius ~ value     {(r_wrong[1]/r_wrong[0])**2:.1f}x"
      "   <- looks 4x bigger")
print(f"area ratio, radius ~ sqrt(value) {(r_right[1]/r_right[0])**2:.1f}x"
      "   <- correct")

fig, axes = plt.subplots(1, 2, figsize=(8, 2.8))
for ax, radii, title in [(axes[0], r_wrong, "radius ~ value (WRONG)"),
                         (axes[1], r_right * 4.5, "area ~ value (right)")]:
    ax.scatter([0, 1], [0, 0], s=(radii * 30) ** 2 / 30, color=BLUE)
    ax.set_title(title, loc="left", fontsize=9)
    ax.set_xlim(-.6, 1.8); ax.set_ylim(-1, 1); ax.axis("off")
fig.savefig("/tmp/claude-0/area_trap.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("\nmatplotlib's `s` parameter is already an AREA in points squared,")
print("which is why passing the raw value is correct there and passing a")
print("radius is the bug. Check the convention in whatever library you use.")

# --- CHANNEL SUITABILITY BY DATA TYPE ---
suitability = pd.DataFrame([
    ("position",   "excellent", "excellent", "excellent", "unlimited"),
    ("length",     "excellent", "good",      "poor",      "unlimited"),
    ("angle",      "fair",      "fair",      "poor",      "unlimited"),
    ("area",       "fair",      "fair",      "poor",      "unlimited"),
    ("lightness",  "good",      "excellent", "poor",      "~7 bins"),
    ("hue",        "poor",      "poor",      "excellent", "5-8 levels"),
    ("shape",      "none",      "none",      "good",      "4-5 levels"),
    ("texture",    "none",      "fair",      "good",      "3-4 levels"),
], columns=["channel", "quantitative", "ordinal", "nominal", "capacity"])
print("\n", suitability.to_string(index=False))

# --- DO NOT DOUBLE-ENCODE: one variable, one channel. ---
fig, axes = plt.subplots(1, 2, figsize=(9.5, 3))
import matplotlib.cm as cm
norm = plt.Normalize(values.min(), values.max())
axes[0].bar(labels, values, color=cm.Blues(norm(values)), width=0.6)
axes[0].set_title("value on BOTH length and colour (wasteful)",
                  loc="left", fontsize=9)
axes[0].set_ylim(0, 45)

# The colour channel is now free to carry a SECOND variable.
is_new = [False, True, False, True, False]
axes[1].bar(labels, values, width=0.6,
            color=[ORANGE if n else BLUE for n in is_new])
axes[1].set_title("length = value, colour = launched this year",
                  loc="left", fontsize=9)
axes[1].set_ylim(0, 45)
handles = [plt.Rectangle((0, 0), 1, 1, color=c) for c in (BLUE, ORANGE)]
axes[1].legend(handles, ["existing", "new"], frameon=False, fontsize=8)
for ax in axes:
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/double_encode.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- CHANNEL BUDGET: how many variables can one chart carry honestly? ---
print("""
  A SCATTER PLOT'S CHANNEL BUDGET

    x            variable 1   continuous, precise
    y            variable 2   continuous, precise
    colour hue   variable 3   nominal, <= 5-7 levels
    size (area)  variable 4   continuous, COARSE differences only
    shape        variable 5   nominal, <= 4-5 levels
    facet        variable 6   nominal, 10-20 panels

  Six is the theoretical ceiling and four is the practical one. Past that,
  facet: position comparison across panels beats a fifth channel inside one.
""")

  THE ACCURACY RANKING (Cleveland & McGill, 1984)

  most accurate
    1  position, common scale     |###  #####  ####   two bars, one axis
    2  position, identical scales  small multiples
    3  length                      stacked-bar segments
    4  angle / slope               pie wedges, line steepness
    5  area                        bubble size, treemap
    6  volume                      3-D anything
    7  colour saturation/lightness heatmap intensity
    8  colour hue                  which category
  least accurate

  Put the variable the reader cares about most on channel 1.

  THE AREA TRAP

    value doubles: 10 -> 20

    radius ~ value              area ~ value
      o        (  )               o       ( )
      r=10     r=20               r=3.2   r=4.5
      area 1x  area 4x            area 1x area 2x
      READS AS 4x                 READS AS ~2x (correct)

    Scale by sqrt(value) when you control the radius. Check your library:
    matplotlib's `s` is already an area.

  ONE VARIABLE, ONE CHANNEL

    BAD                          GOOD
    length = revenue             length = revenue
    colour = revenue             colour = is_new_product
    (colour says nothing new)    (colour carries a second variable)

Key Takeaways

• Human accuracy at reading channels is ranked and measured, so put the most important variable on position along a common scale.
• Area is judged poorly and systematically underestimated, so bubble charts convey only coarse differences.
• Scale marks by area rather than radius, or a doubling looks like a quadrupling.
• Hue carries identity and not magnitude; lightness carries magnitude, which is why sequential scales vary lightness.
• Never encode one variable on two channels, because it spends a channel that could have carried something new.

🧪 Practice

1. Render the same five values as bars, pie wedges, bubbles, and heatmap cells, then ask someone to rank the values from each and record the errors.
2. Draw two circles whose radii are proportional to values of 10 and 20, then redraw with area proportional, and state the perceived ratio in each.
3. Take a chart that colours bars by their own height and rework it so colour carries a second variable.
4. Interview: You need to show revenue, growth rate, region, and customer segment in one view. How do you allocate the channels? (Hint: which two variables deserve the axes, and what do you do with the fourth?)

Color Theory and Accessibility

Colour is the channel people reach for first and reason about least. It is also the one with a hard accessibility constraint: roughly 8% of men and 0.5% of women have some form of colour vision deficiency, most commonly difficulty distinguishing red from green. A chart that encodes meaning in hue alone is unreadable to a meaningful share of any audience.

The discipline that fixes most colour problems is to ask what job the colour is doing. There are only four, and each has one correct structure.

Categorical colour encodes identity: which series is which. The structure is a fixed set of distinguishable hues, assigned in a fixed order. The two rules that matter: assign hues in a fixed sequence so the same entity always gets the same colour, and never generate a ninth hue, because at that point adjacent colours are indistinguishable under colour vision deficiency.

Sequential colour encodes magnitude: how much. The structure is one hue varying from light to dark. Lightness carries the order, which is why a single-hue ramp is readable in greyscale and by every viewer. Never use a rainbow for magnitude: hues have no perceptual order, so a rainbow scale introduces boundaries where the data is smooth and hides differences where it is not. The classic offender is the old jet colormap.

Diverging colour encodes polarity: which side of a meaningful midpoint. The structure is two hues that read as opposites, separated by a neutral grey midpoint. Two requirements: the midpoint must be neutral, because a hue there implies a category that does not exist, and the two poles must be genuinely opposite, so blue against orange or blue against red works while blue against teal does not, both being cool. Diverging scales must also be centred on the real midpoint. Centring a profit-and-loss scale anywhere but zero misstates who is profitable.

Status colour encodes state: good, warning, serious, critical. These meanings are reserved, so a status colour must never double as "series four", and status must always ship with an icon or a label, never colour alone.

Three practices make colour reliable rather than hopeful.

Use perceptually uniform colormaps. Viridis, magma, and cividis are constructed so equal steps in data are equal steps in perceived lightness. The older rainbow maps are not, which means they create visual features that are not in the data.

Never encode by colour alone. Colour should be redundant with something else: a direct label, a shape, a position, a texture, or a line style. This is the single most effective accessibility practice, and it also helps when a chart is printed in greyscale or photocopied.

Check colour choices with a simulator rather than by eye. Colour vision deficiency simulation is a computation, not an aesthetic judgement, and tools exist for it in every language. The test is whether every pair of colours that can appear side by side remains distinguishable under simulated protanopia and deuteranopia. Anything you can compute, you should compute rather than guess.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import to_rgb

# --- A VALIDATED CATEGORICAL PALETTE, in fixed slot order. ---
# Fixed order is the colour-vision-deficiency safety mechanism: the sequence
# is chosen so that ADJACENT slots stay separable. Never shuffle it, and
# never append a 9th hue.
CATEGORICAL = ["#2a78d6",  # 1 blue
               "#eb6834",  # 2 orange
               "#1baf7a",  # 3 aqua
               "#eda100",  # 4 yellow
               "#e87ba4",  # 5 magenta
               "#008300",  # 6 green
               "#4a3aa7",  # 7 violet
               "#e34948"]  # 8 red

# --- SIMULATE COLOUR VISION DEFICIENCY: compute it, do not eyeball it. ---
# Machado, Oliveira & Fernandes (2009) linear transforms at severity 1.0.
CVD = {
    "protanopia": np.array([[0.1522, 0.1178, -0.028],
                            [0.1052, 0.8933,  0.0018],
                            [-0.0038, -0.0291, 1.0329]]),
    "deuteranopia": np.array([[0.3674, 0.8607, -0.2281],
                              [0.2801, 0.6725,  0.0474],
                              [-0.0118, 0.0428, 0.9690]]),
}

def srgb_to_linear(c):
    c = np.asarray(c, dtype=float)
    return np.where(c <= 0.04045, c / 12.92, ((c + 0.055) / 1.055) ** 2.4)

def linear_to_srgb(c):
    c = np.clip(c, 0, 1)
    return np.where(c <= 0.0031308, c * 12.92, 1.055 * c ** (1/2.4) - 0.055)

def simulate(hex_color: str, kind: str) -> tuple:
    lin = srgb_to_linear(to_rgb(hex_color))
    return tuple(linear_to_srgb(CVD[kind] @ lin))

def oklab(hex_color: str) -> np.ndarray:
    """OKLab is roughly perceptually uniform, so Euclidean distance in it
    approximates perceived difference."""
    r, g, b = srgb_to_linear(to_rgb(hex_color))
    l = 0.4122214708*r + 0.5363325363*g + 0.0514459929*b
    m = 0.2119034982*r + 0.6806995451*g + 0.1073969566*b
    s = 0.0883024619*r + 0.2817188376*g + 0.6299787005*b
    l_, m_, s_ = np.cbrt([l, m, s])
    return np.array([
        0.2104542553*l_ + 0.7936177850*m_ - 0.0040720468*s_,
        1.9779984951*l_ - 2.4285922050*m_ + 0.4505937099*s_,
        0.0259040371*l_ + 0.7827717662*m_ - 0.8086757660*s_,
    ])

def delta_e(c1, c2) -> float:
    """Perceptual distance, OKLab scaled by 100. Target >= 8 for adjacent
    series; below ~6 the pair is effectively the same colour."""
    return float(np.linalg.norm(oklab(c1) - oklab(c2)) * 100)

# Check adjacent slot pairs, simulated. This is the test that matters for
# stacked bars and multi-line charts, where only neighbours touch.
print("adjacent slot separation under simulated CVD (target >= 8):")
worst = 1e9
for i in range(len(CATEGORICAL) - 1):
    a, b = CATEGORICAL[i], CATEGORICAL[i+1]
    row = {"pair": f"{i+1}-{i+2}", "normal": round(delta_e(a, b), 1)}
    for kind in CVD:
        sa = "#%02x%02x%02x" % tuple(int(v*255) for v in simulate(a, kind))
        sb = "#%02x%02x%02x" % tuple(int(v*255) for v in simulate(b, kind))
        row[kind[:5]] = round(delta_e(sa, sb), 1)
        worst = min(worst, row[kind[:5]])
    print(f"  {row}")
print(f"worst adjacent pair under CVD: {worst:.1f}")

# A NAIVE palette that looks fine and collapses under simulation.
naive = ["#d62728", "#2ca02c"]            # the classic red/green pair
print(f"\nred vs green, normal vision:  {delta_e(*naive):.1f}")
for kind in CVD:
    sa = "#%02x%02x%02x" % tuple(int(v*255) for v in simulate(naive[0], kind))
    sb = "#%02x%02x%02x" % tuple(int(v*255) for v in simulate(naive[1], kind))
    print(f"red vs green, {kind:12s}: {delta_e(sa, sb):.1f}"
          f"  {'FAILS' if delta_e(sa, sb) < 8 else 'ok'}")

# --- SEQUENTIAL: lightness must increase monotonically. ---
def lightness(hex_color: str) -> float:
    return oklab(hex_color)[0]

print("\nis the ramp monotone in lightness?")
for cmap_name in ("viridis", "Blues", "jet"):
    cmap = plt.get_cmap(cmap_name)
    hexes = ["#%02x%02x%02x" % tuple(int(v*255) for v in cmap(t)[:3])
             for t in np.linspace(0, 1, 9)]
    ls = [lightness(h) for h in hexes]
    diffs = np.diff(ls)
    mono = bool(np.all(diffs > 0) or np.all(diffs < 0))
    print(f"  {cmap_name:8s} monotone: {str(mono):5s}"
          f"  lightness range {min(ls):.2f}-{max(ls):.2f}"
          f"  {'' if mono else '<- creates bands the data does not have'}")

fig, axes = plt.subplots(3, 2, figsize=(10, 4.2))
grid = np.linspace(0, 1, 256).reshape(1, -1)
for row, name in enumerate(["viridis", "Blues", "jet"]):
    axes[row][0].imshow(grid, cmap=name, aspect="auto")
    axes[row][0].set_title(f"{name}", loc="left", fontsize=9)
    cmap = plt.get_cmap(name)
    ls = [lightness("#%02x%02x%02x" % tuple(int(v*255) for v in cmap(t)[:3]))
          for t in np.linspace(0, 1, 60)]
    axes[row][1].plot(ls, color="#0b0b0b", lw=2)
    axes[row][1].set_title("perceived lightness along the ramp",
                           loc="left", fontsize=8)
    for ax in axes[row]:
        ax.set_xticks([]); ax.set_yticks([])
fig.savefig("/tmp/claude-0/ramps.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- DIVERGING: neutral midpoint, and CENTRED ON THE REAL MIDPOINT. ---
profit = np.array([-42.0, -18.0, -3.0, 5.0, 12.0, 30.0])
fig, axes = plt.subplots(1, 2, figsize=(10, 2.8))
from matplotlib.colors import TwoSlopeNorm
axes[0].bar(range(6), profit,
            color=plt.get_cmap("RdBu")(
                plt.Normalize(profit.min(), profit.max())(profit)))
axes[0].axhline(0, color="#0b0b0b", lw=0.8)
axes[0].set_title("diverging NOT centred at 0: WRONG", loc="left", fontsize=9)

norm = TwoSlopeNorm(vmin=profit.min(), vcenter=0.0, vmax=profit.max())
axes[1].bar(range(6), profit, color=plt.get_cmap("RdBu")(norm(profit)))
axes[1].axhline(0, color="#0b0b0b", lw=0.8)
axes[1].set_title("centred at 0: colour agrees with sign", loc="left", fontsize=9)
for ax in axes:
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/diverging.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("\nuncentred, a small LOSS can be shaded the same as a small PROFIT.")
print("Diverging scales must be centred on the value that means 'neither'.")

# --- REDUNDANT ENCODING: never colour alone. ---
rng = np.random.default_rng(2)
x = np.arange(12)
fig, axes = plt.subplots(1, 2, figsize=(11, 3.2), sharey=True)
styles = [("-", "o"), ("--", "s"), (":", "^")]
for i, name in enumerate(["alpha", "beta", "gamma"]):
    y = 10 + i * 3 + np.cumsum(rng.normal(0, .8, 12))
    axes[0].plot(x, y, color=CATEGORICAL[i], lw=2)          # colour only
    ls, mk = styles[i]
    axes[1].plot(x, y, color=CATEGORICAL[i], lw=2, linestyle=ls,
                 marker=mk, markersize=5, label=name)       # colour + style
    axes[1].annotate(name, (x[-1], y[-1]), xytext=(5, 0),
                     textcoords="offset points", fontsize=8,
                     color=CATEGORICAL[i], va="center")
axes[0].set_title("colour alone: unreadable in greyscale or with CVD",
                  loc="left", fontsize=9)
axes[1].set_title("colour + line style + marker + direct label",
                  loc="left", fontsize=9)
axes[1].legend(frameon=False, fontsize=8)
for ax in axes:
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/redundant.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- CONTRAST: marks need >= 3:1 against the surface; text 4.5:1. ---
def rel_luminance(hex_color: str) -> float:
    r, g, b = srgb_to_linear(to_rgb(hex_color))
    return 0.2126 * r + 0.7152 * g + 0.0722 * b

def contrast(a: str, b: str) -> float:
    la, lb = rel_luminance(a), rel_luminance(b)
    hi, lo = max(la, lb), min(la, lb)
    return (hi + 0.05) / (lo + 0.05)

SURFACE = "#fcfcfb"
print(f"\ncontrast against the light surface (marks need >= 3:1):")
for i, c in enumerate(CATEGORICAL, 1):
    ratio = contrast(c, SURFACE)
    print(f"  slot {i} {c}  {ratio:5.2f}:1"
          f"  {'ok' if ratio >= 3 else 'needs a visible label or table view'}")

  THE FOUR JOBS OF COLOUR

  job          encodes            structure                    capacity
  ---------------------------------------------------------------------------
  categorical  identity           fixed distinct hues, in a     5-8 levels
                                  fixed assignment order
  ordinal      position in order  ONE hue, monotone lightness   ~7 steps
  sequential   magnitude          ONE hue, light -> dark        continuous
  diverging    polarity           two opposite hues + NEUTRAL   continuous
                                  grey midpoint
  status       state              reserved scale + icon/label   4 states

  WHY A RAINBOW FAILS FOR MAGNITUDE

  jet:      blue cyan green yellow red
  lightness  .3   .8    .6    .95  .4     <- NOT monotone
                  ^           ^
            bright bands the data does not contain, and dark regions where
            real differences disappear.

  viridis:  dark purple -> teal -> yellow
  lightness  .15  .35  .55  .75  .95      <- monotone, so equal data steps
                                             look like equal steps

  DIVERGING MUST BE CENTRED ON THE REAL MIDPOINT

    values: -42 ... +30, scale centred on the data's mean (-2.7)
      -3 (a LOSS) shades nearly the same as +5 (a PROFIT)   WRONG

    scale centred on 0
      every negative is one hue, every positive the other    RIGHT

  THE ACCESSIBILITY RULES, IN PRIORITY ORDER

    1. never encode by colour ALONE (add label, shape, style, or position)
    2. simulate protanopia and deuteranopia; compute, do not eyeball
    3. one-hue ramps for magnitude, so greyscale printing still works
    4. marks >= 3:1 contrast on the surface, text >= 4.5:1
    5. a legend for 2+ series, and direct labels for <= 4
    6. a table view exists for every chart

Key Takeaways

• Decide what job the colour is doing: identity, order, magnitude, polarity, or status, and each has one correct structure.
• Categorical hues are assigned in a fixed order so an entity keeps its colour, and a ninth generated hue is indistinguishable under colour vision deficiency.
• Sequential scales vary lightness within one hue; rainbow scales are not monotone in lightness and invent bands the data does not have.
• Diverging scales need a neutral midpoint, two genuinely opposite hues, and centring on the value that means "neither".
• Never encode by colour alone, and verify palettes with a simulation rather than by eye, because it is a computation.

🧪 Practice

1. Implement a colour vision deficiency simulation and report the perceptual distance between red and green under normal vision and under deuteranopia.
2. Plot perceived lightness along the viridis, Blues, and jet ramps and identify which are monotone.
3. Render a profit-and-loss bar chart with a diverging scale centred on the data mean and again centred on zero, and describe the misreading the first invites.
4. Take a three-series line chart encoded by colour alone and add two redundant channels, then convert it to greyscale to confirm it still reads.
5. Interview: A stakeholder wants your dashboard to use the company's red and green brand colours for "bad" and "good". How do you respond? (Hint: what share of viewers is affected, and what can carry the meaning alongside hue?)

Scales, Axes, and Annotation

Scales and axes are the chart's units. They tell the reader what the marks mean, and almost every misleading chart is misleading through its axes rather than through its data. Annotation is what converts a chart from a display of numbers into a statement someone can act on.

Scale type is the first decision. A linear scale puts equal data differences at equal distances and is the default. A logarithmic scale puts equal ratios at equal distances, which is correct in three situations: the data spans several orders of magnitude, the meaningful comparison is proportional rather than absolute, or the growth is exponential and you want it to appear straight. A log scale must always be labelled, because a reader who misses it will read growth completely wrong. Note that log scales cannot show zero or negative values, which is a hard constraint rather than an inconvenience.

Axis range carries the asymmetry from the earlier topic: zero is mandatory for bars because they encode length, and optional for lines because they encode slope. Two more range rules. When comparing panels in a small-multiple grid, share the axis range, or the comparison the layout invites is invalid. And when a range is truncated deliberately, say so.

Tick marks and labels should be few and readable. Four to seven ticks is usually right; more becomes visual noise. Use round numbers, format them for humans with thousands separators or abbreviations, and put units in the axis label rather than repeating them on every tick.

Gridlines should be hairline, solid, and one step off the background. They exist to support reading a value, not to be seen. Dashed gridlines add noise and falsely suggest a threshold or a projection.

Then annotation, which is the highest-value and most neglected part of a chart. A chart without a written conclusion makes the reader do the interpretation, and they may not do it the way you would.

The elements worth including:

• A title that states the finding, not the contents. "Checkout conversion fell 1.8 points after the March release" beats "Conversion over time". The title is the one piece of text everyone reads.
• Direct labels on the marks that matter. A label on the endpoint, the extreme, or the one series being discussed. Never a number on every point, which is chaos and goes unread.
• Reference lines for a target, a baseline, an average, or a threshold, so the reader can see whether a value is acceptable rather than just what it is.
• Event markers for a release, a campaign, an outage, or a definition change. A step in a metric with no annotation invites the reader to invent an explanation.
• A source and an as-of timestamp. Which table, which filters, and when it was computed. This is what makes a chart auditable and reconcilable.
• Honest treatment of gaps and incomplete periods. A partial final period looks like a collapse. Either exclude it or mark it explicitly.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter

BLUE, GRAY, INK, MUTED = "#2a78d6", "#c9c9c5", "#0b0b0b", "#52514e"

# --- LINEAR VS LOG: the same data, two different readings. ---
months = pd.date_range("2022-01-01", periods=36, freq="MS")
users = 1000 * 1.14 ** np.arange(36)          # steady 14% monthly growth

fig, axes = plt.subplots(1, 2, figsize=(11, 3.2))
axes[0].plot(months, users, lw=2, color=BLUE)
axes[0].set_title("linear: growth looks like it began in year 3",
                  loc="left", fontsize=9)
axes[1].plot(months, users, lw=2, color=BLUE)
axes[1].set_yscale("log")
axes[1].set_title("log: a straight line = CONSTANT growth rate",
                  loc="left", fontsize=9)
axes[1].set_ylabel("users (log scale)")        # ALWAYS label a log axis
for ax in axes:
    ax.spines[["top", "right"]].set_visible(False)
    ax.grid(axis="y", lw=0.5, color="#e6e6e4", linestyle="-")  # solid hairline
    ax.set_axisbelow(True)
fig.savefig("/tmp/claude-0/log_scale.png", dpi=110, bbox_inches="tight")
plt.close(fig)

growth = pd.Series(users).pct_change().dropna()
print(f"actual monthly growth: min {growth.min():.1%},"
      f" max {growth.max():.1%} -- perfectly constant")
print("The linear panel implies acceleration that does not exist. On a log")
print("axis a constant growth RATE is a straight line, which is the point.")

# --- LOG CANNOT SHOW ZERO OR NEGATIVES ---
with_zero = np.array([0.0, 5.0, 50.0, 500.0])
safe = np.where(with_zero > 0, with_zero, np.nan)
print(f"\nlog10 of the series: {np.log10(safe)}")
print("the zero becomes NaN and silently vanishes from the plot.")
print("Use a symlog scale, or plot the zeros separately and say so.")

# --- SHARED AXES IN SMALL MULTIPLES, or the comparison is invalid. ---
rng = np.random.default_rng(3)
panels = {"north": 100 + np.cumsum(rng.normal(1.0, 2, 24)),
          "south": 40 + np.cumsum(rng.normal(0.3, 1, 24)),
          "east":  180 + np.cumsum(rng.normal(-0.4, 3, 24))}
idx = pd.date_range("2024-01-01", periods=24, freq="MS")

fig, axes = plt.subplots(2, 3, figsize=(12, 4.6))
for col, (name, vals) in enumerate(panels.items()):
    axes[0][col].plot(idx, vals, lw=2, color=BLUE)        # independent scales
    axes[0][col].set_title(f"{name} (own scale)", loc="left", fontsize=8)
    axes[1][col].plot(idx, vals, lw=2, color=BLUE)
    axes[1][col].set_ylim(0, 260)                         # SHARED scale
    axes[1][col].set_title(f"{name} (shared scale)", loc="left", fontsize=8)
for ax in axes.ravel():
    ax.spines[["top", "right"]].set_visible(False)
    ax.tick_params(labelsize=7)
fig.savefig("/tmp/claude-0/shared_axes.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("\nTop row: all three look similar in size. Bottom row: south is a")
print("quarter of east. Independent scales make a grid layout lie.")

# --- TICK FORMATTING: humans, not machines ---
def human(v, _=None):
    for lim, suf in ((1e9, "B"), (1e6, "M"), (1e3, "k")):
        if abs(v) >= lim:
            return f"{v/lim:,.1f}{suf}".replace(".0", "")
    return f"{v:,.0f}"

fig, axes = plt.subplots(1, 2, figsize=(10, 2.9))
big = np.array([1_240_000, 3_800_000, 2_150_000, 5_600_000])
for ax, fmt in zip(axes, [None, human]):
    ax.bar(list("ABCD"), big, color=BLUE, width=0.55)
    if fmt:
        ax.yaxis.set_major_formatter(FuncFormatter(fmt))
        ax.set_yticks(np.linspace(0, 6e6, 4))       # 4 ticks, round numbers
        ax.set_title("4 round ticks, abbreviated", loc="left", fontsize=9)
    else:
        ax.set_title("default: too many ticks, raw digits",
                     loc="left", fontsize=9)
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/ticks.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- ANNOTATION: turn a chart into a statement. ---
conv = np.concatenate([
    np.random.default_rng(9).normal(4.2, .12, 14),     # before the release
    np.random.default_rng(9).normal(2.4, .12, 10),     # after
])
days = pd.date_range("2024-03-01", periods=24, freq="D")
RELEASE = days[14]
TARGET = 4.0

fig, ax = plt.subplots(figsize=(9, 3.8))
ax.plot(days, conv, lw=2, color=BLUE, zorder=3)

# 1. A reference line for the target, labelled in place.
ax.axhline(TARGET, lw=1, color=GRAY, zorder=1)
ax.annotate(f"target {TARGET}%", (days[0], TARGET), xytext=(2, 4),
            textcoords="offset points", fontsize=8, color=MUTED)

# 2. An event marker, so the reader does not invent an explanation.
ax.axvline(RELEASE, lw=1, color="#eb6834", zorder=2)
ax.annotate("v2.4 checkout release", (RELEASE, ax.get_ylim()[1]),
            xytext=(4, -10), textcoords="offset points",
            fontsize=8, color="#eb6834")

# 3. Selective direct labels: the endpoint and the extreme, not every point.
last = conv[-1]
ax.annotate(f"{last:.1f}%", (days[-1], last), xytext=(6, 0),
            textcoords="offset points", fontsize=9, fontweight="bold",
            color=BLUE, va="center")

# 4. A title that states the FINDING, with a subtitle for the mechanism.
ax.set_title("Checkout conversion fell 1.8 points after the v2.4 release",
             loc="left", fontsize=12, color=INK, pad=16)
ax.text(0, 1.04, "Daily conversion, all traffic. The drop begins the day of "
        "the release and has not recovered.",
        transform=ax.transAxes, fontsize=8.5, color=MUTED)

# 5. Source and as-of stamp: what makes the chart auditable.
ax.text(0, -0.22, "Source: analytics.fct_checkout_sessions  |  "
        "filters: bots excluded, web only  |  as of 2024-03-25 06:00 UTC",
        transform=ax.transAxes, fontsize=7, color=MUTED)

ax.set_ylabel("conversion rate (%)")
ax.set_ylim(0, 5.5)
ax.spines[["top", "right"]].set_visible(False)
ax.grid(axis="y", lw=0.5, color="#e6e6e4")
ax.set_axisbelow(True)
fig.savefig("/tmp/claude-0/annotated.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- THE PARTIAL-PERIOD TRAP ---
full = pd.Series([980, 1010, 1040, 1080, 1120],
                 index=pd.period_range("2024-01", periods=5, freq="M"))
partial = full.copy()
partial.iloc[-1] = 380           # month 5 is only 10 days in
fig, axes = plt.subplots(1, 2, figsize=(10, 2.9))
axes[0].bar([str(p) for p in partial.index], partial.values, color=BLUE, width=.55)
axes[0].set_title("partial final month: looks like a collapse",
                  loc="left", fontsize=9)
axes[1].bar([str(p) for p in partial.index[:-1]], partial.values[:-1],
            color=BLUE, width=.55)
axes[1].bar([str(partial.index[-1])], [partial.values[-1]], width=.55,
            color=BLUE, alpha=0.35, hatch="//")
axes[1].annotate("partial\n(10 of 31 days)", (4, 420), fontsize=7,
                 ha="center", color=MUTED)
axes[1].set_title("marked as incomplete", loc="left", fontsize=9)
for ax in axes:
    ax.tick_params(labelsize=7)
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/partial.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("\nAn unmarked partial period is the most common accidental lie in")
print("routine reporting. Exclude it or hatch it, every time.")

  WHEN TO USE A LOG SCALE

  use it when                          do NOT when
  --------------------------------     ------------------------------------
  data spans >= 2 orders of magnitude   values include 0 or negatives
  ratios are the meaningful comparison  the audience will not notice the axis
  growth is exponential                 differences, not ratios, are the point
  Always label the axis "(log scale)".

  CONSTANT GROWTH, TWO AXES

  linear                          log
        |                  /            |        ____----
        |                 /             |  __----
        |               /               |--
        |         _____/                |
        |___-----                       |
        +------------------             +------------------
    "growth started in year 3"     "growth rate never changed"
    Both are the same 14%/month. Only one reads correctly.

  SHARED SCALES IN A GRID

  own scales                      shared scale
  north  south  east              north  south  east
   /\     /\     /\                /\
   ||     ||     ||                ||           /\
   ||     ||     ||                ||    __     ||
  all "the same size"             south is 1/4 of east
  A grid invites comparison, so unshared scales make the layout lie.

  THE ANNOTATION CHECKLIST

    [ ] title states the FINDING, not the contents
    [ ] subtitle gives the population, filter, and mechanism
    [ ] units in the axis label, not on every tick
    [ ] 4-7 round ticks, human-formatted (1.2M, not 1240000)
    [ ] reference line for target / baseline / threshold
    [ ] event markers for releases, campaigns, outages, definition changes
    [ ] direct labels on the endpoint and the extreme only
    [ ] gridlines hairline and SOLID, one step off the surface
    [ ] partial periods excluded or visibly marked
    [ ] source table, filters, and an as-of timestamp

Key Takeaways

• A log scale makes equal ratios equal distances, so constant exponential growth becomes a straight line; always label it, and remember it cannot show zero.
• Zero baselines are mandatory for bars and optional for lines, and panels in a grid must share an axis range or the comparison is invalid.
• Use four to seven round, human-formatted ticks, and keep gridlines hairline and solid so they support reading without being seen.
• A title should state the finding rather than the contents, because it is the one piece of text every reader reads.
• Mark events, targets, and partial periods explicitly, and stamp every chart with its source and as-of time.

🧪 Practice

1. Plot a series growing at a constant percentage on linear and log axes and describe what a reader would conclude from each.
2. Build a three-panel small-multiple grid with independent axes, then with a shared axis, and state which comparison each supports.
3. Take an unannotated line chart and add a finding-based title, a target reference line, an event marker, one direct label, and a source stamp.
4. Render a bar chart whose final period is incomplete, then mark it as partial, and describe the wrong conclusion the first version invites.
5. Interview: A chart shows a metric dropping sharply on one date and a stakeholder asks what caused it. You do not know yet. What do you add to the chart before sending it? (Hint: what would stop the reader from inventing a cause, and what would let someone else reproduce the number?)

Avoiding Misleading Visuals

A chart can be factually accurate and still cause the reader to believe something false. That gap between accurate and honest is where most misleading visuals live, and almost none of them involve fabricated data. They involve a defensible-sounding choice that changes the conclusion.

This matters for a practical reason beyond ethics: you will make most of these mistakes accidentally, because every one of them is a default in some tool or a natural-seeming response to a stakeholder request. Knowing the catalogue is the only reliable defence.

Truncated bar axes are the most common. Bars encode value as length, so starting the axis at 98 instead of 0 makes a 6% difference look like a fourfold one. The chart is accurate and the reader's takeaway is wrong.

Dual axes invent correlations. Because the alignment of two independent scales is arbitrary, you can make two series appear to move together, diverge, or cross wherever you choose. There is no honest version of this chart.

Inconsistent bin widths or intervals distort. A histogram with uneven bins must plot density rather than count, or the wide bins look artificially tall. A time axis with uneven spacing, treated as evenly spaced categories, compresses and stretches the trend arbitrarily.

Cherry-picked ranges are the subtlest, because nothing about the chart is wrong. Choose a start date at a local minimum and any series shows growth. Always ask whether the window was chosen before or after seeing the result.

Area and volume exaggeration multiplies. Scaling an icon's width and height by a value squares the perceived difference; a three-dimensional shape cubes it. Scale by area, and avoid three-dimensional charts entirely, since perspective also makes near items look larger than far ones of the same value.

Reversed or inverted axes invert the message. The notorious example is a chart of firearm deaths with an inverted y-axis, which made a sharp rise look like a decline to anyone reading the shape rather than the numbers.

Aggregation that hides the distribution. A bar of group means conceals whether the groups overlap almost completely. Show the distribution, or the interval, or the sample size.

Missing uncertainty is arguably the biggest one in routine work. A point estimate drawn as a solid line implies precision it does not have. Two bars whose confidence intervals overlap heavily should not be described as different, and a chart without intervals gives the reader no way to know.

Unequal denominators. Comparing raw counts across groups of different sizes is the classic error: a region with more crimes may simply have more people. Rates, not counts, unless the count itself is the question.

The defence is a checklist, applied to your own charts before sending them, because the question is never whether you meant to mislead.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

BLUE, ORANGE, GRAY, MUTED = "#2a78d6", "#eb6834", "#c9c9c5", "#52514e"
rng = np.random.default_rng(8)

# --- 1. TRUNCATED BAR AXIS: same data, two conclusions. ---
vals = np.array([102.0, 108.0])
print("TRUNCATION MATH")
ratio = vals[1] / vals[0]
print(f"  true ratio: {ratio:.2f}x  ({(ratio-1)*100:.0f}% larger)")
for floor in (0, 98, 101):
    drawn = (vals - floor) / (vals[0] - floor) if vals[0] > floor else vals
    print(f"  axis starting at {floor:>3}: bar B looks {drawn[1]:.1f}x bar A")

fig, axes = plt.subplots(1, 3, figsize=(11, 2.9))
for ax, floor in zip(axes, (0, 98, 101)):
    ax.bar(["A", "B"], vals, color=BLUE if floor == 0 else ORANGE, width=.5)
    ax.set_ylim(floor, 110)
    ax.set_title(f"axis from {floor}", loc="left", fontsize=9)
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/truncate.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- 2. CHERRY-PICKED RANGE: one series, four honest-looking stories. ---
full = pd.Series(100 + np.cumsum(rng.normal(0.05, 2.2, 240)),
                 index=pd.date_range("2023-01-01", periods=240, freq="D"))
windows = {"full series": (0, 240), "growth story": None,
           "decline story": None, "flat story": (100, 160)}
# Find the windows that produce each story, which is exactly what a motivated
# analyst does without ever falsifying a number.
best_up = max(((i, j) for i in range(0, 180, 10) for j in range(i+40, 240, 10)),
              key=lambda ij: full.iloc[ij[1]-1] - full.iloc[ij[0]])
best_dn = min(((i, j) for i in range(0, 180, 10) for j in range(i+40, 240, 10)),
              key=lambda ij: full.iloc[ij[1]-1] - full.iloc[ij[0]])
windows["growth story"], windows["decline story"] = best_up, best_dn

fig, axes = plt.subplots(1, 4, figsize=(14, 2.8))
for ax, (name, (i, j)) in zip(axes, windows.items()):
    seg = full.iloc[i:j]
    ax.plot(seg.index, seg.values, lw=2, color=BLUE)
    change = (seg.iloc[-1] / seg.iloc[0] - 1) * 100
    ax.set_title(f"{name}: {change:+.0f}%", loc="left", fontsize=9)
    ax.tick_params(labelsize=6)
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/cherry_pick.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print(f"\nCHERRY-PICKING: one random-walk series yields "
      f"{(full.iloc[best_up[1]-1]/full.iloc[best_up[0]]-1)*100:+.0f}% and "
      f"{(full.iloc[best_dn[1]-1]/full.iloc[best_dn[0]]-1)*100:+.0f}% windows.")
print("Ask whether the window was chosen before or after seeing the answer.")

# --- 3. UNEQUAL BINS: a histogram of counts becomes a lie. ---
x = rng.exponential(10, 4000)
edges = [0, 2, 4, 6, 8, 40]                     # the last bin is 16x wider
counts, _ = np.histogram(x, bins=edges)
widths = np.diff(edges)
fig, axes = plt.subplots(1, 2, figsize=(9.5, 2.9))
axes[0].bar(edges[:-1], counts, width=widths, align="edge", color=ORANGE,
            edgecolor="white")
axes[0].set_title("unequal bins, plotting COUNT: wrong", loc="left", fontsize=9)
axes[1].bar(edges[:-1], counts / widths, width=widths, align="edge",
            color=BLUE, edgecolor="white")
axes[1].set_title("same bins, plotting DENSITY: right", loc="left", fontsize=9)
for ax in axes:
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/bins.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print(f"\nUNEQUAL BINS: the widest bin holds {counts[-1]} of {len(x)} values")
print(f"over {widths[-1]} units. As a count it towers; as a density it is the")
print("flattest bin in the chart, which is the truth.")

# --- 4. MISSING UNCERTAINTY: bars of means hide overlapping distributions. ---
a = rng.normal(52, 14, 60)
b = rng.normal(56, 14, 60)
from scipy import stats
t = stats.ttest_ind(a, b)
fig, axes = plt.subplots(1, 3, figsize=(12, 3.1))
axes[0].bar(["A", "B"], [a.mean(), b.mean()], color=BLUE, width=.5)
axes[0].set_title("means only: 'B is higher'", loc="left", fontsize=9)

se = [a.std(ddof=1)/np.sqrt(len(a)), b.std(ddof=1)/np.sqrt(len(b))]
axes[1].bar(["A", "B"], [a.mean(), b.mean()], color=BLUE, width=.5,
            yerr=[1.96*s for s in se], capsize=6, ecolor="#0b0b0b")
axes[1].set_title("with 95% intervals: they overlap", loc="left", fontsize=9)

for i, (name, d) in enumerate([("A", a), ("B", b)]):
    axes[2].scatter(np.full(len(d), i) + rng.normal(0, .06, len(d)), d,
                    s=10, alpha=.45, color=BLUE)
axes[2].set_xticks([0, 1], ["A", "B"])
axes[2].set_title("the distributions: almost identical", loc="left", fontsize=9)
for ax in axes:
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig("/tmp/claude-0/uncertainty.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print(f"\nMISSING UNCERTAINTY: means {a.mean():.1f} vs {b.mean():.1f},"
      f" p = {t.pvalue:.2f}. The bar chart implies a difference the data")
print("does not support. Error bars or a dot plot would have shown it.")

# --- 5. UNEQUAL DENOMINATORS: counts versus rates. ---
regions = pd.DataFrame({
    "region": ["Metro", "Suburb", "Rural"],
    "incidents": [4200, 1800, 300],
    "population": [2_100_000, 600_000, 60_000],
})
regions["per_100k"] = (regions["incidents"] / regions["population"] * 1e5).round(0)
print("\nUNEQUAL DENOMINATORS:\n", regions.to_string(index=False))
print("By count, Metro is worst by far. By rate, Rural is the highest.")
print("Counts answer 'where is the most volume'; rates answer 'where is the")
print("highest risk'. State which question you are answering.")

# --- 6. A SELF-AUDIT FUNCTION ---
def audit_chart(spec: dict) -> list[str]:
    """Run before sending any chart. Every rule here comes from a real
    misleading-chart failure mode."""
    issues = []
    if spec.get("geom") == "bar" and spec.get("y_min", 0) != 0:
        issues.append("CRITICAL: truncated axis on a length encoding")
    if spec.get("dual_axis"):
        issues.append("CRITICAL: dual axis invents a correlation; index instead")
    if spec.get("y_inverted"):
        issues.append("CRITICAL: inverted axis reverses the visual message")
    if spec.get("uneven_bins") and spec.get("stat") == "count":
        issues.append("plot density, not count, with unequal bin widths")
    if spec.get("aggregate") and not spec.get("shows_uncertainty"):
        issues.append("means shown without intervals or sample sizes")
    if spec.get("compares_groups") and not spec.get("normalised"):
        issues.append("comparing counts across groups of different size")
    if spec.get("partial_final_period") and not spec.get("partial_marked"):
        issues.append("incomplete final period will read as a collapse")
    if spec.get("faceted") and not spec.get("shared_scale"):
        issues.append("grid layout invites comparison but scales differ")
    if spec.get("n_series", 1) > 8:
        issues.append("more than 8 series: fold the tail or facet")
    if spec.get("colour_only_encoding"):
        issues.append("encoded by colour alone: add label, shape, or style")
    if spec.get("log_scale") and not spec.get("axis_labelled_log"):
        issues.append("log axis not labelled as such")
    if not spec.get("source"):
        issues.append("no source or as-of stamp: not auditable")
    return issues

bad = {"geom": "bar", "y_min": 98, "dual_axis": True, "aggregate": True,
       "compares_groups": True, "n_series": 11, "colour_only_encoding": True}
print("\naudit of a deliberately bad chart spec:")
for issue in audit_chart(bad):
    print("  -", issue)

good = {"geom": "bar", "y_min": 0, "aggregate": True, "shows_uncertainty": True,
        "compares_groups": True, "normalised": True, "n_series": 3,
        "source": "analytics.fct_orders", "partial_final_period": False}
print(f"\naudit of a careful spec: {audit_chart(good) or 'no issues'}")

  THE CATALOGUE

  technique                what it does                    severity
  ---------------------------------------------------------------------------
  truncated bar axis       exaggerates ratios              critical
  dual axis                invents a correlation           critical
  inverted axis            reverses the message            critical
  3-D / volume encoding    exaggerates by the cube         critical
  cherry-picked window     manufactures any trend          high
  count with uneven bins   inflates the wide bins          high
  unequal denominators     confuses volume with rate       high
  no uncertainty shown     implies precision               high
  unmarked partial period  reads as a collapse             high
  unshared facet scales    makes the grid layout lie       medium
  aggregation only         hides overlapping distributions medium
  colour-only encoding     excludes ~8% of male readers    medium

  TRUNCATION, MEASURED

    values 102 and 108 -- a 6% difference

    axis from 0    A|######  B|#######        reads as 6%     honest
    axis from 98   A|##      B|########       reads as 4x     misleading
    axis from 101  A|#       B|#######        reads as 7x     egregious

  THE SAME RANDOM WALK, THREE STORIES

    full:   ---\/--\/-\/--\/---   +1%   "flat"
    window: ------/               +18%  "strong growth"
    window:        \-----         -14%  "worrying decline"
  No number was altered. Only the start and end dates were chosen.

  THE PRE-SEND CHECK

    1. is any length encoding truncated?
    2. is there a second y-axis?
    3. was the date range chosen before or after seeing the result?
    4. are the denominators comparable?
    5. is uncertainty visible?
    6. is the final period complete?
    7. would the opposite conclusion also fit this chart?

Key Takeaways

• Most misleading charts contain accurate data and one defensible-sounding choice that changes the conclusion.
• Truncating a bar axis exaggerates ratios, and a dual axis invents a correlation whose crossing point you control.
• With unequal bin widths, plot density rather than count, and compare rates rather than counts across groups of different size.
• A point estimate without uncertainty implies precision it does not have, and overlapping intervals mean you cannot call the values different.
• Audit your own charts against the catalogue before sending, since the question is never whether you intended to mislead.

🧪 Practice

1. Render two bars of 102 and 108 with the axis starting at 0, 98, and 101, and compute the ratio a reader would infer from each.
2. Generate a random walk with no real trend, then find the windows producing the largest apparent growth and largest apparent decline.
3. Build a histogram with one bin sixteen times wider than the rest, plot it as counts and as density, and explain the difference.
4. Take a bar chart of two group means whose confidence intervals overlap and produce two alternative versions that show the overlap.
5. Implement the chart audit function and run it against three charts from a dashboard you use.
6. Interview: A stakeholder asks you to "zoom in so the improvement is visible" on a bar chart. How do you handle it? (Hint: what encoding do bars use, and what alternative form shows a small change honestly?)

<a id="73-visualization-practice"></a>

7.3 Visualization Practice

Principles need tools, and the tools differ sharply in what they make easy; this section covers the library landscape, the standard statistical plots, and the specialised forms for interaction, geography, networks, and faceted comparison.

Static Plotting Libraries

The Python plotting landscape looks chaotic until you notice that the libraries sit at different levels of abstraction over the same job, and that almost all of them are built on one foundation.

Matplotlib is that foundation. It is imperative and low-level: you create a figure, add axes, and draw marks one call at a time. Everything is controllable, which means everything must be specified, and a publication-quality figure takes real effort. Two API styles coexist and mixing them causes confusion, so pick one. The pyplot state-machine style, where functions act on an implicit "current axes", is fine for quick exploration. The object-oriented style, where you hold a figure and axes explicitly, is what you should use for anything reusable, because it composes into functions and subplots without hidden state.

Seaborn sits on matplotlib and adds statistical defaults. It understands tidy dataframes, so you pass column names rather than arrays, and it computes the statistics a plot implies: confidence intervals, regression fits, kernel density estimates, and category aggregation. Its figure-level functions produce faceted grids in one call. Because the result is still matplotlib objects, you can drop down to matplotlib for fine control, which is the combination most Python analysts settle on.

plotnine implements the grammar of graphics directly, giving Python ggplot2's API. If you think in the grammar, this is the most direct expression of it.

Plotly and Altair produce interactive output, covered in the next topic but worth placing here: Plotly is imperative and JavaScript-backed, Altair is declarative and compiles to Vega-Lite.

Two practical matters matter more than library choice.

Backends and headless rendering. Matplotlib needs a backend, and in a script, container, or CI job there is no display. Setting the Agg backend before importing pyplot makes rendering to a file work anywhere. Forgetting this is the most common cause of a plotting script that works in a notebook and crashes in production. Also close figures explicitly in loops, or memory grows until the process dies.

Output format. Vector formats such as SVG and PDF scale to any size and keep text selectable, which is what you want for documents and slides. Raster formats such as PNG are right for the web and for scatter plots with hundreds of thousands of points, where a vector file would contain a path per point and become unopenable. Set an explicit resolution for raster output, because the default is usually too low for print.

A note on the R side, since statistical work often crosses over: ggplot2 is the reference implementation of the grammar, and its consistency is the reason the grammar spread. If you learn one of ggplot2 and plotnine you have learned both.

import numpy as np
import pandas as pd

# --- BACKEND FIRST, before pyplot is imported. This is the headless fix. ---
import matplotlib
matplotlib.use("Agg")            # render to file; no display needed
import matplotlib.pyplot as plt

BLUE, ORANGE, GRAY, INK, MUTED = "#2a78d6", "#eb6834", "#c9c9c5", "#0b0b0b", "#52514e"
rng = np.random.default_rng(0)

df = pd.DataFrame({
    "month": pd.date_range("2024-01-01", periods=12, freq="MS").tolist() * 3,
    "region": np.repeat(["north", "south", "east"], 12),
})
df["revenue"] = (100 + np.tile(np.arange(12), 3) * 3
                 + np.repeat([30, 0, 15], 12) + rng.normal(0, 6, 36)).round(1)

# --- TWO MATPLOTLIB API STYLES. Use the second for anything reusable. ---
# 1. pyplot state machine: acts on an implicit "current axes".
plt.figure(figsize=(5, 2.5))
plt.plot(df.loc[df.region == "north", "month"],
         df.loc[df.region == "north", "revenue"])
plt.title("pyplot style")
plt.savefig("/tmp/claude-0/style_pyplot.png", dpi=110, bbox_inches="tight")
plt.close()

# 2. Object-oriented: explicit figure and axes, no hidden state.
fig, ax = plt.subplots(figsize=(5, 2.5))
ax.plot(df.loc[df.region == "north", "month"],
        df.loc[df.region == "north", "revenue"], lw=2, color=BLUE)
ax.set_title("object-oriented style", loc="left")
fig.savefig("/tmp/claude-0/style_oo.png", dpi=110, bbox_inches="tight")
plt.close(fig)
# The OO style is what lets you write a function that returns an Axes and
# drop it into any grid position. Mixing the two is the usual source of
# "why did my plot land on the wrong subplot".

# --- A REUSABLE STYLE, set once rather than repeated per chart. ---
plt.rcParams.update({
    "figure.dpi": 110,
    "font.size": 10,
    "axes.titlesize": 11,
    "axes.titlelocation": "left",
    "axes.spines.top": False,
    "axes.spines.right": False,
    "axes.grid": True,
    "axes.axisbelow": True,
    "grid.color": "#e6e6e4",
    "grid.linewidth": 0.5,
    "grid.linestyle": "-",           # solid: dashed grids add noise
    "legend.frameon": False,
    "figure.facecolor": "#fcfcfb",
    "axes.facecolor": "#fcfcfb",
    "axes.prop_cycle": plt.cycler(color=[BLUE, ORANGE, "#1baf7a", "#eda100"]),
})

# --- A CHART FUNCTION: the unit of reuse. Takes an Axes, returns it. ---
def revenue_lines(data: pd.DataFrame, ax=None, highlight: str | None = None):
    """Plot revenue by region. If `highlight` is given, use emphasis."""
    if ax is None:
        _, ax = plt.subplots(figsize=(7, 3.2))
    palette = {"north": BLUE, "south": ORANGE, "east": "#1baf7a"}
    for region, g in data.groupby("region"):
        is_focus = highlight is None or region == highlight
        ax.plot(g["month"], g["revenue"], lw=2 if is_focus else 1.3,
                color=palette[region] if is_focus else GRAY,
                zorder=3 if is_focus else 1)
        if is_focus:
            last = g.iloc[-1]
            ax.annotate(region, (last["month"], last["revenue"]),
                        xytext=(5, 0), textcoords="offset points",
                        fontsize=8, color=palette[region], va="center")
    ax.set_ylabel("revenue (USD k)")
    return ax

fig, axes = plt.subplots(1, 2, figsize=(12, 3.2), sharey=True)
revenue_lines(df, axes[0])
revenue_lines(df, axes[1], highlight="north")
axes[0].set_title("all series")
axes[1].set_title("emphasis on north")
fig.savefig("/tmp/claude-0/reusable.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- SEABORN: statistical defaults over tidy data. ---
import seaborn as sns
sns.set_theme(style="whitegrid", rc={"grid.linewidth": 0.5})

# One call computes the aggregation AND a bootstrap confidence band.
fig, axes = plt.subplots(1, 3, figsize=(14, 3.2))
sns.lineplot(data=df, x="month", y="revenue", hue="region",
             errorbar=("ci", 95), ax=axes[0])
axes[0].set_title("lineplot: aggregation + CI computed for you")

sns.boxplot(data=df, x="region", y="revenue", ax=axes[1])
sns.stripplot(data=df, x="region", y="revenue", ax=axes[1],
              size=3, color=INK, alpha=.5)     # overlay the raw points
axes[1].set_title("box + strip: summary and raw data together")

sns.regplot(data=df, x=df["month"].map(pd.Timestamp.toordinal), y="revenue",
            scatter_kws={"s": 12, "alpha": .6}, ax=axes[2])
axes[2].set_title("regplot: fit + interval in one call")
fig.savefig("/tmp/claude-0/seaborn.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# FIGURE-LEVEL functions build the whole grid and return a FacetGrid, not Axes.
g = sns.relplot(data=df, x="month", y="revenue", col="region",
                kind="line", height=2.6, aspect=1.2)
g.set_titles("{col_name}")
g.figure.savefig("/tmp/claude-0/seaborn_facet.png", dpi=110, bbox_inches="tight")
plt.close(g.figure)
# Axes-level (boxplot, lineplot) take `ax=`; figure-level (relplot, catplot,
# displot) own the whole figure. Passing `ax=` to a figure-level function is
# a common error, because it does not accept one.

# --- OUTPUT FORMAT: vector for documents, raster for dense scatter. ---
import os
dense = rng.normal(size=(120_000, 2))
fig, ax = plt.subplots(figsize=(4, 4))
ax.scatter(dense[:, 0], dense[:, 1], s=1, alpha=.15, color=BLUE,
           rasterized=True)     # rasterize the DATA, keep the text vector
fig.savefig("/tmp/claude-0/dense.pdf", bbox_inches="tight", dpi=150)
fig.savefig("/tmp/claude-0/dense_vector.svg", bbox_inches="tight")
plt.close(fig)
for path in ("/tmp/claude-0/dense.pdf", "/tmp/claude-0/dense_vector.svg"):
    print(f"  {os.path.basename(path):20s} {os.path.getsize(path)/1e6:6.2f} MB")
print("`rasterized=True` on the data layer keeps a PDF small while leaving")
print("axis text selectable. Without it, 120k points become 120k vector paths.")

# --- MEMORY: close figures in a loop or the process grows without bound. ---
for i in range(3):
    fig, ax = plt.subplots()
    ax.plot([0, 1], [0, i])
    fig.savefig(f"/tmp/claude-0/loop_{i}.png")
    plt.close(fig)                    # <- omit this and figures accumulate
print(f"\nopen figures after the loop: {len(plt.get_fignums())}")

# --- THE GRAMMAR, IN PYTHON, via plotnine ---
print("""
  plotnine (ggplot2's grammar in Python):

    from plotnine import ggplot, aes, geom_line, facet_wrap, labs, theme_minimal
    (ggplot(df, aes("month", "revenue", color="region"))
     + geom_line(size=1)
     + facet_wrap("~region")
     + labs(title="Revenue by region", y="Revenue (USD k)")
     + theme_minimal())
""")

  THE ABSTRACTION STACK

  level        library      you specify              it computes
  ---------------------------------------------------------------------------
  grammar      plotnine     mapping + geom + facet   stats, scales, legends
  statistical  seaborn      columns + plot kind      aggregation, CI, KDE, fit
  imperative   matplotlib   every mark, by hand      nothing
                    |
                    v
  seaborn and plotnine both RENDER through matplotlib, so you can always
  drop down for fine control. That is why the usual answer is "seaborn for
  the plot, matplotlib for the polish".

  CHOOSING

  need                                    reach for
  ---------------------------------------------------------------------------
  full control, publication figure        matplotlib, object-oriented
  statistical plot with CI or a fit       seaborn
  faceted grid in one call                seaborn figure-level (relplot/catplot)
  grammar-of-graphics thinking            plotnine (Python) / ggplot2 (R)
  interactive HTML                        plotly / altair
  100k+ points                            datashader, or rasterized=True

  THE TWO MATPLOTLIB GOTCHAS

  1. BACKEND        matplotlib.use("Agg") BEFORE importing pyplot,
                    or a headless script raises or hangs.
  2. FIGURE LEAK    plt.close(fig) inside every loop, or memory grows
                    until the process is killed.

  AXES-LEVEL VS FIGURE-LEVEL (seaborn)

    axes-level    lineplot, boxplot, scatterplot   accept ax=, compose freely
    figure-level  relplot, catplot, displot        own the figure, do faceting
    Passing ax= to a figure-level function is an error, not an oversight.

Key Takeaways

• Matplotlib is the foundation and seaborn and plotnine render through it, so the usual workflow is a high-level call plus matplotlib polish.
• Use matplotlib's object-oriented API for anything reusable, since the pyplot state machine hides which axes you are drawing on.
• Set the Agg backend before importing pyplot for scripts and containers, and close figures inside loops to avoid unbounded memory growth.
• Seaborn computes the statistics a plot implies, and its figure-level functions own the whole figure rather than accepting an axes.
• Vector output for documents, raster or rasterized data layers for dense scatter plots.

🧪 Practice

1. Write a chart function that accepts an optional axes argument and returns it, then place it into both a single figure and a two-panel grid.
2. Configure a style dictionary once with hairline solid gridlines, no top or right spines, and a fixed colour cycle, and apply it to three charts.
3. Save a 100,000-point scatter as an SVG and as a PDF with the data layer rasterized, and compare the file sizes.
4. Interview: A colleague's plotting script works in a notebook and fails in the nightly job with a display error. What is your first fix? (Hint: what does matplotlib need before it can draw, and what does a container lack?)

Statistical Plots: Distributions and Relationships

Beyond the basic chart types there is a set of plots designed to answer specific statistical questions. Knowing them is largely a matter of knowing which question each one answers, because reaching for the wrong one wastes the reader's attention on a view that cannot show what you need.

For a single distribution, the ladder runs from crude to precise. A histogram bins and counts, and its appearance depends on bin width. A kernel density estimate smooths instead of binning, which trades the bin-width choice for a bandwidth choice and can invent smoothness that is not there, including mass below zero for a strictly positive variable. A rug plot or strip plot shows every observation and works only for small samples. The ECDF involves no smoothing parameter at all and is the most honest single view, at the cost of being less immediately intuitive.

For comparing distributions across groups, the box plot is compact and sortable and hides modality. The violin plot shows the full shape at the cost of the bandwidth choice. The strip and swarm plots show every point, which is ideal below a few hundred observations per group. A ridgeline plot stacks densities vertically and handles many groups gracefully. The reliable pattern in practice is box or violin plus the raw points overlaid, which gives you summary and detail together, with the sample size visible.

For a relationship, the scatter plot is the base, and the additions that matter are a fitted line with an interval, and a residual plot. The residual plot deserves emphasis because it is the diagnostic that reveals what the scatter hides: plotting residuals against fitted values should show a formless cloud, and any pattern in it means the model is wrong in a specific way. A curve means a missing non-linear term; a funnel means the variance changes with the level.

For a matrix of relationships, the pair plot gives every pairwise scatter plus each variable's own distribution on the diagonal.

Two more views worth having in reach. The QQ plot compares a sample against a theoretical distribution by plotting quantile against quantile, so deviation from the diagonal shows exactly where and how a distribution departs from normality, which a histogram shows only vaguely. The bootstrap is the general-purpose way to get an interval for any statistic, including ones with no formula, by resampling with replacement.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats

BLUE, ORANGE, GRAY, INK = "#2a78d6", "#eb6834", "#c9c9c5", "#0b0b0b"
sns.set_theme(style="whitegrid", rc={"grid.linewidth": .5})
rng = np.random.default_rng(1)

# --- ONE DISTRIBUTION, FIVE VIEWS, one bimodal sample. ---
x = np.concatenate([rng.normal(30, 6, 400), rng.normal(58, 7, 300)])

fig, axes = plt.subplots(1, 5, figsize=(16, 2.9))
axes[0].hist(x, bins=30, color=BLUE)
axes[0].set_title("histogram (bin width)", fontsize=9)

sns.kdeplot(x=x, ax=axes[1], color=BLUE, fill=True)
axes[1].set_title("KDE (bandwidth)", fontsize=9)

axes[2].boxplot(x, widths=.4)
axes[2].set_title("box (hides both peaks)", fontsize=9)

sns.violinplot(y=x, ax=axes[3], color=BLUE)
axes[3].set_title("violin (shows them)", fontsize=9)

xs = np.sort(x)
axes[4].step(xs, np.arange(1, len(xs)+1)/len(xs), color=BLUE, lw=2)
axes[4].set_title("ECDF (no parameter)", fontsize=9)
fig.savefig("/tmp/claude-0/one_dist.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# THE BANDWIDTH PROBLEM mirrors the bin-width problem.
print("KDE bandwidth changes the apparent number of modes:")
from scipy.stats import gaussian_kde
grid = np.linspace(x.min(), x.max(), 400)
for bw in (0.15, 0.4, 1.2):
    dens = gaussian_kde(x, bw_method=bw)(grid)
    peaks = sum(1 for i in range(1, len(dens)-1)
                if dens[i] > dens[i-1] and dens[i] > dens[i+1])
    print(f"  bw_method={bw:<5} apparent modes: {peaks}")

# A KDE ON BOUNDED DATA leaks mass past the boundary.
positive = rng.exponential(5, 800)               # strictly positive
kde = gaussian_kde(positive)
below_zero = kde.integrate_box_1d(-np.inf, 0)
print(f"\nKDE of a strictly positive variable puts {below_zero:.1%} of its"
      " mass below zero")
print("A histogram cannot do this. For bounded data prefer a histogram or ECDF.")

# --- COMPARING GROUPS: box/violin PLUS the points, with n visible. ---
groups = pd.DataFrame({
    "team": np.repeat(["alpha", "beta", "gamma", "delta"], [120, 90, 25, 140]),
})
offsets = {"alpha": 50, "beta": 62, "gamma": 58, "delta": 45}
groups["score"] = [rng.normal(offsets[t], 9) for t in groups["team"]]

order = groups.groupby("team")["score"].median().sort_values().index  # sort!
fig, axes = plt.subplots(1, 3, figsize=(14, 3.4))
sns.boxplot(data=groups, x="team", y="score", order=order, ax=axes[0],
            color=GRAY, width=.5)
axes[0].set_title("box, sorted by median", fontsize=9)

sns.violinplot(data=groups, x="team", y="score", order=order, ax=axes[1],
               color=GRAY, inner="quartile")
axes[1].set_title("violin: full shape", fontsize=9)

sns.boxplot(data=groups, x="team", y="score", order=order, ax=axes[2],
            color="#f0f0ee", width=.5, showfliers=False)
sns.stripplot(data=groups, x="team", y="score", order=order, ax=axes[2],
              size=2.5, alpha=.5, color=BLUE, jitter=.25)
# Print n under each category: a box from 25 points deserves less trust.
counts = groups["team"].value_counts()
axes[2].set_xticks(range(len(order)),
                   [f"{t}\nn={counts[t]}" for t in order], fontsize=8)
axes[2].set_title("box + points + n: what to actually ship", fontsize=9)
fig.savefig("/tmp/claude-0/groups.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- RELATIONSHIP: the scatter, and the RESIDUAL PLOT that audits it. ---
xr = rng.uniform(1, 10, 400)
y_curved = 3 + 0.4 * xr ** 2 + rng.normal(0, 2, 400)          # missing x^2
y_funnel = 3 + 2.2 * xr + rng.normal(0, 1, 400) * xr          # variance grows
y_good = 3 + 2.2 * xr + rng.normal(0, 2.5, 400)               # well specified

fig, axes = plt.subplots(2, 3, figsize=(13, 5.4))
for col, (name, yy) in enumerate([("well specified", y_good),
                                  ("missing curvature", y_curved),
                                  ("changing variance", y_funnel)]):
    slope, intercept = np.polyfit(xr, yy, 1)
    fitted = intercept + slope * xr
    resid = yy - fitted
    axes[0][col].scatter(xr, yy, s=9, alpha=.5, color=BLUE)
    axes[0][col].plot(np.sort(xr), intercept + slope * np.sort(xr),
                      color=ORANGE, lw=2)
    axes[0][col].set_title(f"{name}  (R2={np.corrcoef(xr, yy)[0,1]**2:.2f})",
                           fontsize=9)
    axes[1][col].scatter(fitted, resid, s=9, alpha=.5, color=BLUE)
    axes[1][col].axhline(0, color=INK, lw=.8)
    axes[1][col].set_title("residuals vs fitted", fontsize=9)
    axes[1][col].set_xlabel("fitted")
fig.savefig("/tmp/claude-0/residuals.png", dpi=110, bbox_inches="tight")
plt.close(fig)

print("\nR-squared for all three fits is respectable. The residual plots")
print("distinguish them: a formless cloud, a clear arc, and a widening funnel.")
# Quantify what the eye sees, so it goes in a report rather than a hunch.
for name, yy in [("well specified", y_good), ("missing curvature", y_curved),
                 ("changing variance", y_funnel)]:
    slope, intercept = np.polyfit(xr, yy, 1)
    resid = yy - (intercept + slope * xr)
    lo, hi = resid[xr < np.median(xr)], resid[xr >= np.median(xr)]
    # Correlation of residuals with x^2 catches curvature; the SD ratio
    # catches heteroscedasticity.
    print(f"  {name:20s} corr(resid, x^2) {np.corrcoef(resid, xr**2)[0,1]:+.2f}"
          f"   sd ratio hi/lo {hi.std()/lo.std():.2f}")

# --- QQ PLOT: where and how a distribution departs from normal. ---
samples = {"normal": rng.normal(0, 1, 500),
           "right-skewed": rng.lognormal(0, .7, 500),
           "heavy-tailed": rng.standard_t(3, 500)}
fig, axes = plt.subplots(1, 3, figsize=(12, 3.1))
for ax, (name, s) in zip(axes, samples.items()):
    stats.probplot(s, dist="norm", plot=ax)
    ax.set_title(f"{name}  (Shapiro p={stats.shapiro(s[:400]).pvalue:.1e})",
                 fontsize=9)
    ax.get_lines()[0].set(markersize=3, color=BLUE, alpha=.6)
    ax.get_lines()[1].set(color=ORANGE, lw=1.5)
fig.savefig("/tmp/claude-0/qq.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("\nA QQ plot says WHERE the departure is: skew bends one end, heavy")
print("tails lift both. A histogram only says 'not quite normal'.")

# --- PAIR PLOT: every pairwise relationship at once. ---
feats = pd.DataFrame({"a": rng.normal(0, 1, 400)})
feats["b"] = .8 * feats["a"] + rng.normal(0, .6, 400)
feats["c"] = feats["a"] ** 2 + rng.normal(0, .4, 400)
feats["group"] = rng.choice(["x", "y"], 400)
pg = sns.pairplot(feats, hue="group", diag_kind="kde", height=1.9,
                  plot_kws={"s": 10, "alpha": .5})
pg.figure.savefig("/tmp/claude-0/pairplot.png", dpi=100, bbox_inches="tight")
plt.close(pg.figure)

# --- BOOTSTRAP: an interval for ANY statistic, no formula needed. ---
def bootstrap_ci(data, stat_fn, n_boot=4000, conf=.95, seed=0):
    r = np.random.default_rng(seed)
    boots = np.array([stat_fn(r.choice(data, len(data), replace=True))
                      for _ in range(n_boot)])
    lo, hi = np.percentile(boots, [(1-conf)/2*100, (1+conf)/2*100])
    return stat_fn(data), lo, hi

skewed = rng.lognormal(3, .8, 300)
print("\nbootstrap 95% intervals (no closed-form needed):")
for label, fn in [("median", np.median),
                  ("90th percentile", lambda d: np.percentile(d, 90)),
                  ("IQR", lambda d: np.percentile(d, 75) - np.percentile(d, 25)),
                  ("Gini-ish ratio", lambda d: np.percentile(d, 90)/np.median(d))]:
    est, lo, hi = bootstrap_ci(skewed, fn)
    print(f"  {label:16s} {est:8.2f}   [{lo:7.2f}, {hi:7.2f}]")

  WHICH VIEW ANSWERS WHICH QUESTION

  question                              plot            free parameter
  ---------------------------------------------------------------------------
  what shape is this variable?          histogram       bin width
  ... smoothly?                         KDE             bandwidth
  ... with no arbitrary choice?         ECDF            NONE
  where is the mass, exactly?           rug / strip     none (small n only)
  how do groups compare, compactly?     box             none
  ... including shape?                  violin          bandwidth
  ... with every point visible?         strip / swarm   none (small n)
  ... across many groups?               ridgeline       bandwidth
  is there a relationship?              scatter         none
  is my MODEL right?                    residuals vs    none
                                        fitted
  is this normal, and how not?          QQ plot         none
  all pairs at once?                    pair plot       none
  what is the interval on <statistic>?  bootstrap       n_boot

  THE RESIDUAL PLOT IS THE DIAGNOSTIC

  well specified        missing curvature      changing variance
    .  . . .  .            .                      .        .
   . .. . .. ..           . ..            ..     . .   .  . . .
  ---------------      ---....----....--    ---.-.-.--.-.-.---
   . . .. . . .           ..    ..              .  . .    .  .
    .  .  . .                ..                             .
  formless cloud       an ARC: add a term   a FUNNEL: transform y
                                            or use robust errors
  All three can have a respectable R-squared. Only this plot separates them.

  BOX PLOT VS VIOLIN: WHAT THE BOX HIDES

    bimodal data          box plot of it
     ##    ##             +------+
    ####  ####            |      |    identical to the box of a
    ##########            |  --  |    uniform distribution over
                          +------+    the same range

Key Takeaways

• Histograms trade a bin-width choice, kernel density estimates trade a bandwidth choice, and the ECDF has no free parameter at all.
• A kernel density estimate of a strictly positive variable puts mass below zero, so prefer a histogram or ECDF for bounded data.
• Ship box or violin plots with the raw points overlaid and the sample size printed, and sort groups by median rather than alphabetically.
• The residual plot is the diagnostic a scatter plot hides: an arc means a missing non-linear term, a funnel means changing variance.
• A QQ plot shows where and how a distribution departs from normal, and the bootstrap gives an interval for any statistic without a formula.

🧪 Practice

1. Plot one bimodal sample as a histogram, kernel density estimate, box plot, violin, and ECDF, and note which views hide the second peak.
2. Vary a kernel density bandwidth across three values and report the apparent number of modes at each.
3. Fit a linear model to data with a genuine quadratic term, then show that the R-squared looks acceptable while the residual plot reveals the arc.
4. Write a bootstrap function and use it to put a confidence interval on a 90th percentile and on an interquartile range.
5. Interview: A colleague reports R-squared of 0.88 and concludes the model fits well. Which single plot do you ask for, and what three things could it reveal? (Hint: plot what the model got wrong, against what it predicted.)

Interactive Visualization

An interactive chart lets the reader ask follow-up questions without you. That is its entire value, and it is also the source of its main failure: interaction often substitutes for a clear default view, producing a chart that shows nothing until the reader works for it.

The useful rule is that the default view must answer the main question on its own, and interaction handles the follow-ups. If the reader must hover to learn anything, the chart has failed regardless of how good the tooltip is.

The interaction patterns worth knowing, roughly in order of value:

Tooltips on hover give per-mark detail without cluttering the chart. The critical constraint is that tooltips enhance and never gate: every value a tooltip reveals must also be reachable another way, through a direct label, an axis, or a table view. A value available only on hover is invisible to keyboard users, to anyone on touch, and in a printout or screenshot, which is how most charts actually travel.

Two tooltip mechanics matter. On line and area charts use a crosshair that snaps to the nearest x position and shows every series at once, so the reader aims at a date rather than at a two-pixel line. On bars, cells, and dots the mark itself is the target, and the hit area must be larger than the mark, since an eight-pixel dot is a pinpoint nobody reliably hits.

Zoom and pan matter when the data has detail at multiple scales. Always provide a reset control, because a reader who zooms into an empty region has no way back.

Brushing and linking is the most genuinely powerful pattern: selecting points in one chart highlights the same records in every other chart. This is how you explore a multivariate dataset, because it lets the reader ask "who are those outliers on this axis, and where do they sit on the others?"

Filtering controls belong in one row above everything they scope, never inside a chart card and never per-chart. Every chart below re-renders against the same slice, so the numbers always agree. Per-chart filters are how dashboards start disagreeing with themselves.

Cross-filtering makes clicking a mark filter the rest of the view, which is powerful and easy to make confusing without a visible indication of active filters and a way to clear them.

On the tooling side: Plotly is imperative and produces standalone HTML, which is the quickest route from a dataframe to a shareable interactive chart. Altair is declarative over Vega-Lite and expresses linked selections particularly cleanly. Bokeh targets larger data and server-backed apps. D3 offers total control at a much higher cost. For very large datasets, datashader aggregates server-side into an image, which is the only approach that works at tens of millions of points.

One caution that bites in practice: an interactive chart typically embeds the data in the page. A million rows becomes an unusably large HTML file, and it also means any sensitive column you plotted is now readable in the page source. Aggregate before rendering.

import numpy as np
import pandas as pd

rng = np.random.default_rng(2)
n = 900
df = pd.DataFrame({
    "region": rng.choice(["north", "south", "east", "west"], n),
    "segment": rng.choice(["smb", "mid", "enterprise"], n),
    "spend": rng.gamma(3, 120, n).round(2),
    "seats": rng.integers(1, 400, n),
    "month": rng.choice(pd.date_range("2024-01-01", periods=12, freq="MS"), n),
})
df["revenue"] = (df["spend"] * rng.uniform(.8, 1.4, n)).round(2)

# --- PLOTLY: dataframe to standalone interactive HTML. ---
import plotly.express as px
import plotly.graph_objects as go

PALETTE = {"north": "#2a78d6", "south": "#eb6834",
           "east": "#1baf7a", "west": "#eda100"}

fig = px.scatter(
    df, x="seats", y="revenue", color="region",
    color_discrete_map=PALETTE,          # fixed per entity, never by rank
    # hover_data controls the tooltip contents explicitly.
    hover_data={"segment": True, "spend": ":.2f", "month": "|%b %Y"},
    labels={"seats": "Seats", "revenue": "Revenue (USD)"},
    title="Revenue vs seats by region",
)
fig.update_traces(marker={"size": 8, "opacity": .7,
                          "line": {"width": 1, "color": "white"}})  # 2px ring
fig.update_layout(
    template="simple_white", hovermode="closest",
    legend={"orientation": "h", "y": 1.02, "title": None},
)
fig.write_html("/tmp/claude-0/scatter.html", include_plotlyjs="cdn")
# include_plotlyjs="cdn" keeps the file small; "inline" makes it work offline
# at the cost of ~3 MB per file.

# --- A CROSSHAIR TOOLTIP for time series: aim at a date, not at a line. ---
monthly = (df.groupby(["month", "region"])["revenue"].sum()
             .reset_index())
fig2 = go.Figure()
for region, g in monthly.groupby("region"):
    fig2.add_trace(go.Scatter(
        x=g["month"], y=g["revenue"], name=region, mode="lines+markers",
        line={"width": 2, "color": PALETTE[region]}, marker={"size": 7},
    ))
fig2.update_layout(
    template="simple_white",
    # "x unified" shows EVERY series at the hovered date in one tooltip.
    hovermode="x unified",
    xaxis={"showspikes": True, "spikemode": "across", "spikethickness": 1},
    title="Monthly revenue by region",
)
fig2.write_html("/tmp/claude-0/lines.html", include_plotlyjs="cdn")

# --- ALTAIR: declarative, and linked selections are one line. ---
import altair as alt
alt.data_transformers.disable_max_rows()

brush = alt.selection_interval(encodings=["x", "y"])      # drag to select
legend_click = alt.selection_point(fields=["region"], bind="legend")

scatter = (alt.Chart(df).mark_circle(size=60, opacity=.65)
           .encode(
               x=alt.X("seats:Q", title="Seats"),
               y=alt.Y("revenue:Q", title="Revenue (USD)"),
               # Selected points keep their hue; the rest go gray.
               color=alt.condition(brush, "region:N", alt.value("#d8d8d4"),
                                   scale=alt.Scale(
                                       domain=list(PALETTE),
                                       range=list(PALETTE.values()))),
               opacity=alt.condition(legend_click, alt.value(.7), alt.value(.06)),
               tooltip=["region", "segment", "seats", "revenue"])
           .add_params(brush, legend_click)
           .properties(width=340, height=280))

# BRUSHING AND LINKING: this bar chart recomputes from the brushed selection.
bars = (alt.Chart(df).mark_bar()
        .encode(x=alt.X("sum(revenue):Q", title="Revenue in selection"),
                y=alt.Y("segment:N", sort="-x"),
                color=alt.Color("region:N", scale=alt.Scale(
                    domain=list(PALETTE), range=list(PALETTE.values()))))
        .transform_filter(brush)              # <- the link
        .properties(width=340, height=280))

(scatter | bars).properties(
    title="Drag on the scatter to filter the bars"
).save("/tmp/claude-0/linked.html")
print("wrote linked.html: dragging a region of the scatter recomputes the bars")

# --- THE GATING RULE: never let a tooltip be the only path to a value. ---
print("""
  TOOLTIPS ENHANCE, THEY NEVER GATE

  BAD                              GOOD
  bare marks; every number         axis + selective direct labels carry the
  lives in the tooltip             main values; the tooltip adds detail
                                   AND a table view exists

  Why it matters: a tooltip is unavailable to keyboard-only users, awkward
  on touch, and absent from the screenshot that will be pasted into the
  deck someone actually presents.
""")

# --- A TABLE VIEW is the accessible twin of every chart. ---
table = (df.groupby(["region", "segment"])
           .agg(customers=("revenue", "size"), revenue=("revenue", "sum"))
           .round(0).reset_index())
print("the table view that must accompany the chart:")
print(table.head(6).to_string(index=False))

# --- FILE SIZE: an interactive chart embeds the data. ---
import os
big = pd.DataFrame({"x": rng.normal(size=120_000), "y": rng.normal(size=120_000)})
px.scatter(big, x="x", y="y").write_html("/tmp/claude-0/big.html",
                                         include_plotlyjs="cdn")
agg = big.assign(xb=pd.cut(big.x, 60), yb=pd.cut(big.y, 60)).groupby(
    ["xb", "yb"], observed=True).size().reset_index(name="n")
print(f"\n120k raw points as HTML: "
      f"{os.path.getsize('/tmp/claude-0/big.html')/1e6:.1f} MB")
print(f"aggregated to a {len(agg)}-cell density grid instead: "
      f"{len(agg)} rows to embed")
print("Aggregate before rendering. The browser has to parse whatever you ship,")
print("and every plotted column becomes readable in the page source.")

  INTERACTION PATTERNS, BY VALUE

  pattern            gives the reader              cost / caveat
  ---------------------------------------------------------------------------
  tooltip on hover   per-mark detail               must not GATE any value
  crosshair (time)   all series at one x           only for ordered x
  zoom + pan         multi-scale detail            must ship a reset control
  brush + link       "who are those points?"       the most valuable pattern
  filter row         one consistent slice          ONE row, above everything
  cross-filter       drill by clicking a mark      needs visible active state
  legend toggle      hide/show a series            colour must not re-map

  THE DEFAULT-VIEW RULE

    the default view answers the MAIN question alone
              |
              +-- interaction answers FOLLOW-UP questions

    If the chart is blank or meaningless until hovered, no amount of
    tooltip polish fixes it.

  HIT TARGETS

    an 8px dot              the same dot with a hit area
        o                        +--------+
                                 |   o    |   >= 24px, includes the ring;
                                 +--------+   dense scatter uses a
                                              nearest-point layer

  LIBRARY CHOICE

  need                                    reach for
  ---------------------------------------------------------------------------
  fast dataframe -> shareable HTML         plotly express
  declarative spec, linked selections      altair (Vega-Lite)
  server-backed app, bigger data           bokeh
  total control, bespoke design            d3
  10M+ points                              datashader (aggregate to an image)

  THE EMBEDDED-DATA TRAP
    an interactive HTML file CONTAINS its data. 120k rows is a heavy page
    and every plotted column is readable in the source. Aggregate first.

Key Takeaways

• The default view must answer the main question by itself, with interaction reserved for follow-up questions.
• Tooltips enhance and never gate, because a hover-only value is unreachable by keyboard, awkward on touch, and absent from a screenshot.
• Use a snapping crosshair showing all series on time-series charts, and make hit areas larger than the marks they belong to.
• Brushing and linking is the highest-value pattern for multivariate exploration, and filters belong in one row scoping everything below.
• Interactive files embed their data, so aggregate first for both file size and disclosure reasons.

🧪 Practice

1. Build an interactive scatter plot with an explicit tooltip field list and a fixed colour map keyed by entity rather than by rank.
2. Create a linked pair of charts where brushing the first recomputes the second.
3. Produce a time-series chart with a unified crosshair tooltip showing every series at the hovered date.
4. Render 100,000 points to interactive HTML, record the file size, then aggregate to a binned grid and compare.
5. Interview: A stakeholder says your interactive dashboard is "hard to read" although every value is in a tooltip. What is the likely problem? (Hint: what does the chart show before anyone touches it?)

Geospatial Mapping

When location is part of the question, a map is often the only view that answers it, because spatial patterns such as clustering, corridors, and borders are invisible in a table. Maps also carry failure modes no other chart has, and most of them come from the fact that geography is not neutral.

The first is projection. The Earth is a sphere and a screen is flat, and no flat map can preserve area, shape, distance, and direction simultaneously. Every projection sacrifices something. Web Mercator, the default for nearly every online map, preserves local shape and direction, which is what you want for navigation, and grossly inflates area away from the equator: Greenland appears comparable to Africa and is about a fourteenth of its size. If your map is about area or about comparing regions at different latitudes, use an equal-area projection instead.

The second is the choropleth's denominator problem, which is really the unequal-denominator error from earlier with a geographic face. Shading regions by a raw count produces a map of where people live. Every such map looks the same, because population dominates everything. Always map a rate, a density, or a per-capita figure unless the count itself is genuinely the question.

Related and less obvious: a choropleth gives equal visual weight to unequal areas. A large rural region with few people dominates the map while a dense urban area with far more people is a dot. That is an inversion of importance, and the standard alternatives are a proportional symbol map, which sizes a mark by value at each location, a cartogram, which distorts the geography so that area encodes the value, or a hexbin map, which imposes equal-area cells.

The third is binning. Choropleths need the continuous value cut into classes, and the class boundaries change the map's appearance dramatically. Equal intervals, quantiles, natural breaks, and standard deviations give four different maps from one dataset. State which you used, and prefer a scheme whose boundaries are meaningful in the domain over one chosen because it looks good.

The main map types and their jobs:

Type	Encodes	Good for
Choropleth	value per region, by colour	rates over administrative areas
Proportional symbol	value per point, by area	counts at locations
Dot density	one dot per n units	distribution and texture
Heatmap / kernel density	smoothed point intensity	clustering of events
Hexbin	value per equal-area cell	dense points, unbiased by area size
Flow map	movement between places	migration, trade, traffic
Cartogram	value as area	when population should dominate

On tooling, the key concepts are the geometry column, holding points, lines, and polygons, and the coordinate reference system, which says what the coordinates mean. GeoPandas adds both to pandas. The most common practical bug in all of geospatial work is mismatched reference systems, where two layers use different coordinate definitions and simply do not overlap, or a distance is computed in degrees and reported as though it were metres.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

BLUE, ORANGE, INK, MUTED = "#2a78d6", "#eb6834", "#0b0b0b", "#52514e"
rng = np.random.default_rng(4)

# --- THE DENOMINATOR PROBLEM: a count map is a population map. ---
regions = pd.DataFrame({
    "region": ["Metro", "Inner", "Suburb", "Coastal", "Rural", "Highland"],
    "population": [2_400_000, 900_000, 620_000, 180_000, 74_000, 21_000],
    "area_km2": [320, 540, 1_450, 2_100, 9_800, 14_200],
})
# Underlying risk is genuinely HIGHEST in the least populated places.
true_rate = np.array([180, 210, 240, 300, 420, 510])          # per 100k
regions["cases"] = (regions["population"] / 1e5 * true_rate).round(0).astype(int)
regions["per_100k"] = (regions["cases"] / regions["population"] * 1e5).round(0)
regions["pop_density"] = (regions["population"] / regions["area_km2"]).round(0)

print(regions.to_string(index=False))
print("\nranked by CASE COUNT: ", regions.sort_values("cases", ascending=False)
      ["region"].tolist())
print("ranked by RATE:       ", regions.sort_values("per_100k", ascending=False)
      ["region"].tolist())
print("\nThe two orderings are exact opposites. A count choropleth would show")
print("Metro as the problem; the rate shows Highland. Only one answers 'where")
print("is the risk highest'.")
print(f"\ncorrelation of cases with population: "
      f"{regions['cases'].corr(regions['population']):+.3f}"
      "  <- a count map IS a population map")

# --- BINNING SCHEME CHANGES THE MAP ---
values = regions["per_100k"].values
schemes = {}
k = 4
schemes["equal interval"] = np.linspace(values.min(), values.max(), k + 1)
schemes["quantile"] = np.quantile(values, np.linspace(0, 1, k + 1))
schemes["std deviation"] = np.array([
    values.min(), values.mean() - values.std(), values.mean(),
    values.mean() + values.std(), values.max()])
print("\nthe same six values, three class schemes:")
for name, edges in schemes.items():
    labels = pd.cut(values, bins=np.unique(edges), include_lowest=True,
                    labels=False)
    print(f"  {name:15s} edges {np.round(np.unique(edges), 0).tolist()}")
    print(f"  {'':15s} class per region {labels.tolist()}")
print("\nHighland and Rural share a class under one scheme and split under")
print("another. State the scheme on the map.")

# --- PROJECTION: Web Mercator inflates area away from the equator. ---
def mercator_area_factor(lat_deg: float) -> float:
    """Web Mercator scales area by 1/cos^2(latitude)."""
    return 1.0 / np.cos(np.radians(lat_deg)) ** 2

places = pd.DataFrame({
    "place": ["Equator", "Cairo", "Paris", "Stockholm", "Greenland (mid)"],
    "lat": [0, 30, 48.9, 59.3, 72],
})
places["area_inflation"] = places["lat"].map(mercator_area_factor).round(2)
print("\n", places.to_string(index=False))
print("\nA Greenland-latitude polygon is drawn ~10x its true area relative to")
print("an equatorial one. Never compare region SIZES on a Mercator map, and")
print("use an equal-area projection when area is the point.")

# --- A SCHEMATIC CHOROPLETH VS A PROPORTIONAL SYMBOL MAP ---
# Grid positions stand in for geography; the point is the encoding, not the
# map accuracy. Cell sizes deliberately differ, as real regions do.
layout = {                     # (x, y, width, height)
    "Metro":    (4.2, 2.2, 0.8, 0.8),
    "Inner":    (3.4, 2.0, 0.8, 1.2),
    "Suburb":   (3.0, 0.8, 2.0, 1.2),
    "Coastal":  (0.2, 0.8, 2.6, 2.2),
    "Rural":    (0.2, 3.2, 4.8, 1.6),
    "Highland": (0.2, 4.8, 4.8, 1.8),
}
fig, axes = plt.subplots(1, 3, figsize=(14, 3.8))
import matplotlib.patches as mpatches
from matplotlib.cm import ScalarMappable

for ax, (col, title, cmap) in zip(axes, [
        ("cases", "choropleth of COUNTS (a population map)", "Oranges"),
        ("per_100k", "choropleth of RATES (the real pattern)", "Blues"),
        (None, "proportional symbols: area encodes count", None)]):
    if col:
        norm = plt.Normalize(regions[col].min(), regions[col].max())
        cm = plt.get_cmap(cmap)
        for _, r in regions.iterrows():
            x, y, w, h = layout[r["region"]]
            ax.add_patch(mpatches.Rectangle(
                (x, y), w, h, facecolor=cm(norm(r[col])),
                edgecolor="white", linewidth=2))     # 2px surface gap
            ax.text(x + w/2, y + h/2, r["region"][:4], ha="center",
                    va="center", fontsize=7, color=INK)
        fig.colorbar(ScalarMappable(norm=norm, cmap=cm), ax=ax,
                     fraction=.04, pad=.02)
    else:
        for _, r in regions.iterrows():
            x, y, w, h = layout[r["region"]]
            # s is an AREA in matplotlib, so pass the value, not its root.
            ax.scatter([x + w/2], [y + h/2], s=r["cases"] / 6,
                       color=BLUE, alpha=.65, edgecolor="white", linewidth=1.5)
            ax.text(x + w/2, y + h/2 - .45, r["region"][:4], ha="center",
                    fontsize=7, color=MUTED)
        for _, r in regions.iterrows():             # faint region outlines
            x, y, w, h = layout[r["region"]]
            ax.add_patch(mpatches.Rectangle((x, y), w, h, facecolor="none",
                                            edgecolor="#e6e6e4", linewidth=1))
    ax.set_title(title, loc="left", fontsize=9)
    ax.set_xlim(0, 5.6); ax.set_ylim(6.8, 0); ax.axis("off")
fig.savefig("/tmp/claude-0/maps.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("\nNote in panel 1 and 2 how the two largest cells (Rural, Highland)")
print("dominate the ink while holding 2% of the population. That visual")
print("over-weighting is inherent to choropleths.")

# --- HEXBIN: equal-area cells, so region size cannot bias the reading. ---
lon = np.concatenate([rng.normal(-0.12, .06, 4000),     # a dense cluster
                      rng.uniform(-0.45, .25, 1500)])   # diffuse background
lat = np.concatenate([rng.normal(51.51, .04, 4000),
                      rng.uniform(51.35, 51.65, 1500)])
fig, axes = plt.subplots(1, 2, figsize=(10, 3.8))
axes[0].scatter(lon, lat, s=2, alpha=.12, color=BLUE)
axes[0].set_title("raw points: overplotted", loc="left", fontsize=9)
hb = axes[1].hexbin(lon, lat, gridsize=28, cmap="Blues", mincnt=1)
axes[1].set_title("hexbin: equal-area cells, honest density",
                  loc="left", fontsize=9)
fig.colorbar(hb, ax=axes[1], fraction=.04, label="points per cell")
for ax in axes:
    ax.set_xlabel("longitude"); ax.set_ylabel("latitude")
fig.savefig("/tmp/claude-0/hexbin_map.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- COORDINATE REFERENCE SYSTEMS: the most common geospatial bug. ---
print("""
  THE CRS BUG, IN ONE EXAMPLE

    EPSG:4326  degrees of latitude/longitude   (what GPS and GeoJSON give you)
    EPSG:3857  metres, Web Mercator            (what tile maps use)
    EPSG:326xx metres, UTM zone                (what you want for distances)

    gdf.crs                       # ALWAYS check before any spatial operation
    gdf.to_crs(epsg=3857)         # reproject, do not relabel
    gdf.set_crs(epsg=4326)        # only when the CRS was MISSING, never to convert

  Two symptoms of getting it wrong:
    1. two layers plot in completely different places, or one is a dot
    2. a "distance" of 0.03 is degrees, not metres -- off by ~100,000x

  A degree of longitude is ~111 km at the equator and ~0 km at the poles,
  so a buffer or distance in degrees is not a fixed distance at all.
""")

# Rough great-circle distance, for when you must work in degrees.
def haversine_km(lat1, lon1, lat2, lon2):
    R = 6371.0
    p1, p2 = np.radians(lat1), np.radians(lat2)
    dp, dl = p2 - p1, np.radians(lon2 - lon1)
    a = np.sin(dp/2)**2 + np.cos(p1) * np.cos(p2) * np.sin(dl/2)**2
    return 2 * R * np.arcsin(np.sqrt(a))

print(f"  London to Paris: {haversine_km(51.5074, -0.1278, 48.8566, 2.3522):.0f} km")
print(f"  naive euclidean on raw degrees: "
      f"{np.hypot(51.5074-48.8566, -0.1278-2.3522):.2f} -- a meaningless number")

  WHY A COUNT CHOROPLETH IS ALWAYS THE SAME MAP

  cases by region          rate per 100k
  Metro     4,320  ####    Metro      180  ##
  Inner     1,890  ##      Inner      210  ##
  Suburb    1,488  ##      Suburb     240  ###
  Coastal     540  #       Coastal    300  ###
  Rural       311  .       Rural      420  #####
  Highland    107  .       Highland   510  ######

  Exactly reversed. corr(cases, population) = +0.99, so a count map is a
  population map with extra steps.

  THE CHOROPLETH AREA BIAS

    +--------------------------------+
    |                                |   Highland: 14,200 km2, 21k people
    |          HIGHLAND              |   -> dominates the ink
    +--------------------------------+
    |           RURAL                |   Rural: 9,800 km2, 74k people
    +----------+---------+-----------+
    | COASTAL  | SUBURB  |In|Metro|      Metro: 320 km2, 2.4M people
    +----------+---------+--+-----+      -> nearly invisible

    Fixes: proportional symbols, a cartogram, or hexbins.

  PROJECTION TRADE-OFFS

  projection        preserves            inflates            use for
  ---------------------------------------------------------------------------
  Web Mercator      local shape, angle   area, badly, by     web tiles,
                                         1/cos^2(lat)        navigation
  Albers equal-area area                 shape at edges      comparing regions
  Robinson          a compromise         nothing exactly     world overviews
  UTM zone          distance, locally    fails across zones  metric distances

    area inflation on Web Mercator: 1.0x at the equator, 2.3x at Paris,
    10.5x at Greenland's latitude.

Key Takeaways

• Shade a choropleth by a rate or density, never a raw count, because a count map is a population map.
• Choropleths give equal visual weight to unequal areas, so large sparse regions dominate; proportional symbols, cartograms, and hexbins fix that.
• Class boundaries change a choropleth's appearance substantially, so choose a scheme deliberately and state it.
• Web Mercator inflates area by the inverse square of the cosine of latitude, so never compare region sizes on it.
• Always check the coordinate reference system before a spatial operation, and never treat degrees as a distance.

🧪 Practice

1. Build a region table where the ranking by count is the exact reverse of the ranking by rate, and state which question each answers.
2. Classify one set of regional rates with equal-interval, quantile, and standard-deviation schemes and list which regions change class.
3. Compute the Web Mercator area inflation factor at the equator, at 45 degrees, and at 70 degrees latitude.
4. Render dense point data as a raw scatter and as a hexbin, and explain what the second recovers.
5. Interview: A public-health map shades counties by case count and leadership wants to send resources to the darkest counties. What is your objection? (Hint: what other variable is almost perfectly correlated with any raw count by region?)

Network and Hierarchy Visualization

Some data is about relationships rather than values. A social graph, a dependency tree, a set of money transfers, an organisational chart: in all of these the connections are the subject, and a table of edges conveys almost nothing about the structure while a drawing can convey it instantly.

Networks, or graphs, consist of nodes and edges. The standard drawing is a node-link diagram, and its central problem is that there is no natural position for a node: unlike a scatter plot, where position means something, a graph layout must be invented. That is what layout algorithms do.

Force-directed layouts are the common default. They simulate physics: edges act as springs pulling connected nodes together, and all nodes repel each other. The result places densely connected groups near each other, which makes community structure visible. Two things about them matter. They are non-deterministic, so running twice gives different pictures of the same graph, and you must fix the random seed for a reproducible figure. And distance in the layout is not a quantity: two nodes drawn close together are probably related, but the distance does not measure anything, and readers routinely over-interpret it.

The hard limit is density. Past a few hundred nodes with moderate connectivity, a node-link diagram becomes a hairball with no readable structure. The alternatives when you hit that point:

• An adjacency matrix, with nodes on both axes and cells shaded by edge presence or weight. It scales much further, it makes density obvious, and clusters appear as blocks if you order the rows sensibly. It loses the intuitive sense of paths.
• Aggregate the graph. Draw communities as single nodes with weighted edges between them. Usually the right move, because the question is normally about groups rather than individuals.
• Filter. Show the largest component, or the strongest edges, or the neighbourhood of one node of interest.

Before drawing, compute the centrality measures, because they tell you which nodes matter and what to label. Degree centrality counts connections. Betweenness counts how often a node lies on shortest paths, identifying brokers whose removal disconnects the graph. Closeness measures average distance to everyone else. Eigenvector centrality and PageRank weight connections by the importance of who they connect to. These often disagree, and the disagreement is informative: a node with low degree and high betweenness is a bridge.

Hierarchies are the special case of a graph with no cycles and a single parent per node, and their extra structure permits better layouts. A node-link tree shows the structure clearly and wastes space. Treemaps nest rectangles whose area encodes a value, using space efficiently and encoding on area, which is judged poorly, so they show coarse proportions rather than precise ones. Sunburst charts are treemaps in polar coordinates, with the same trade-off plus angle. Icicle plots keep the rectangles but lay them out linearly, which is easier to read than a sunburst. Dendrograms show the merge order from hierarchical clustering, and Sankey diagrams show flow through stages, which is the natural form for funnels and for any conserved quantity moving between categories.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import networkx as nx

BLUE, ORANGE, GRAY, INK = "#2a78d6", "#eb6834", "#c9c9c5", "#0b0b0b"
rng = np.random.default_rng(3)

# --- A GRAPH WITH REAL COMMUNITY STRUCTURE AND ONE BRIDGE NODE ---
G = nx.Graph()
communities = {"a": range(0, 9), "b": range(9, 18), "c": range(18, 26)}
for name, members in communities.items():
    members = list(members)
    for i, u in enumerate(members):
        G.add_node(u, community=name)
        for v in members[i+1:]:
            if rng.random() < 0.55:                # dense WITHIN a community
                G.add_edge(u, v)
# A single low-degree node that connects all three groups: the broker.
G.add_node(99, community="bridge")
for anchor in (2, 11, 20):
    G.add_edge(99, anchor)
G.add_edge(4, 13)                                  # one extra cross edge

print(f"{G.number_of_nodes()} nodes, {G.number_of_edges()} edges, "
      f"density {nx.density(G):.3f}")

# --- CENTRALITY FIRST: it tells you what to label. ---
cent = pd.DataFrame({
    "degree": pd.Series(dict(G.degree())),
    "betweenness": pd.Series(nx.betweenness_centrality(G)).round(3),
    "closeness": pd.Series(nx.closeness_centrality(G)).round(3),
    "pagerank": pd.Series(nx.pagerank(G)).round(3),
})
print("\ntop 5 by betweenness:")
print(cent.sort_values("betweenness", ascending=False).head(5).to_string())
print("\nnode 99 has LOW degree and the HIGHEST betweenness: it is a bridge.")
print("Degree alone would have ranked it near the bottom.")

# Removing the bridge splits the graph, which is what betweenness predicted.
H = G.copy()
H.remove_node(99)
print(f"\ncomponents with node 99: {nx.number_connected_components(G)}")
print(f"components without it:   {nx.number_connected_components(H)}")

# --- LAYOUT: fix the seed, or the same graph looks different every run. ---
pos_a = nx.spring_layout(G, seed=1)
pos_b = nx.spring_layout(G, seed=2)
shift = np.mean([np.linalg.norm(pos_a[n] - pos_b[n]) for n in G])
print(f"\nmean node displacement between two unseeded-equivalent layouts:"
      f" {shift:.2f}")
print("Force-directed layouts are non-deterministic. Pin the seed for any")
print("figure you will publish or compare against.")

PALETTE = {"a": BLUE, "b": ORANGE, "c": "#1baf7a", "bridge": "#e34948"}
fig, axes = plt.subplots(1, 3, figsize=(15, 4.4))

pos = nx.spring_layout(G, seed=7, k=0.5)
sizes = [120 + 900 * cent.loc[n, "betweenness"] for n in G]   # area ~ value
colors = [PALETTE[G.nodes[n]["community"]] for n in G]
nx.draw_networkx_edges(G, pos, ax=axes[0], alpha=.25, width=1)
nx.draw_networkx_nodes(G, pos, ax=axes[0], node_size=sizes, node_color=colors,
                       edgecolors="white", linewidths=1.5)
# Label only the nodes that matter, never all of them.
top = cent.nlargest(3, "betweenness").index
nx.draw_networkx_labels(G, pos, ax=axes[0],
                        labels={n: str(n) for n in top}, font_size=9)
axes[0].set_title("node-link: size = betweenness, colour = community",
                  loc="left", fontsize=9)
axes[0].axis("off")

# --- ADJACENCY MATRIX: scales past the hairball, clusters become blocks. ---
order = sorted(G.nodes(), key=lambda n: (G.nodes[n]["community"], n))
A = nx.to_numpy_array(G, nodelist=order)
axes[1].imshow(A, cmap="Blues", interpolation="nearest")
axes[1].set_title("adjacency matrix, ordered by community", loc="left", fontsize=9)
axes[1].set_xticks([]); axes[1].set_yticks([])
# The three dense blocks on the diagonal ARE the communities.

# --- AGGREGATE THE GRAPH: usually the question is about groups. ---
C = nx.Graph()
for u, v in G.edges():
    cu, cv = G.nodes[u]["community"], G.nodes[v]["community"]
    if C.has_edge(cu, cv):
        C[cu][cv]["weight"] += 1
    else:
        C.add_edge(cu, cv, weight=1)
cpos = nx.spring_layout(C, seed=3)
widths = [C[u][v]["weight"] * 0.6 for u, v in C.edges()]
nx.draw_networkx_edges(C, cpos, ax=axes[2], width=widths, alpha=.5)
nx.draw_networkx_nodes(C, cpos, ax=axes[2],
                       node_size=[400 + 90 * len(list(communities.get(n, [1])))
                                  for n in C],
                       node_color=[PALETTE[n] for n in C],
                       edgecolors="white", linewidths=2)
nx.draw_networkx_labels(C, cpos, ax=axes[2], font_size=9)
axes[2].set_title("aggregated to communities: readable at any scale",
                  loc="left", fontsize=9)
axes[2].axis("off")
fig.savefig("/tmp/claude-0/network.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- WHEN DOES A NODE-LINK DIAGRAM STOP WORKING? ---
print("\nedge count as a graph grows at fixed average degree 6:")
for nodes in (50, 200, 1000, 5000):
    edges = nodes * 6 // 2
    print(f"  {nodes:>5} nodes -> {edges:>6} edges"
          f"   {'readable' if nodes <= 200 else 'HAIRBALL: use a matrix'}")

# --- HIERARCHY: the same tree, three encodings. ---
tree = pd.DataFrame([
    ("Engineering", "Platform", 38), ("Engineering", "Product", 52),
    ("Engineering", "Data", 21), ("Sales", "Enterprise", 44),
    ("Sales", "SMB", 29), ("Marketing", "Growth", 17),
    ("Marketing", "Brand", 11),
], columns=["division", "team", "headcount"])

fig, axes = plt.subplots(1, 3, figsize=(15, 3.8))

# 1. Node-link tree: structure is explicit, space use is poor.
T = nx.DiGraph()
for _, r in tree.iterrows():
    T.add_edge(r["division"], r["team"])
    T.add_edge("Company", r["division"])
tpos = nx.nx_agraph.graphviz_layout(T, prog="dot") if False else \
    nx.spring_layout(T, seed=5)
nx.draw(T, tpos, ax=axes[0], with_labels=True, node_size=700,
        node_color=GRAY, font_size=6, edgecolors="white")
axes[0].set_title("node-link tree: clear structure, wasted space",
                  loc="left", fontsize=9)

# 2. Treemap: area encodes the value, space is used fully.
import matplotlib.patches as mpatches
def squarify(values, x, y, w, h):
    """A simple slice-and-dice treemap: alternate the split direction."""
    total, out, horizontal = sum(values), [], w > h
    offset = 0
    for v in values:
        frac = v / total
        if horizontal:
            out.append((x + offset, y, w * frac, h))
            offset += w * frac
        else:
            out.append((x, y + offset, w, h * frac))
            offset += h * frac
    return out

div_totals = tree.groupby("division")["headcount"].sum().sort_values(ascending=False)
div_rects = squarify(div_totals.values, 0, 0, 10, 6)
for (dx, dy, dw, dh), (div, _) in zip(div_rects, div_totals.items()):
    teams = tree[tree["division"] == div].sort_values("headcount", ascending=False)
    for (tx, ty, tw, th), (_, tr) in zip(
            squarify(teams["headcount"].values, dx, dy, dw, dh), teams.iterrows()):
        axes[1].add_patch(mpatches.Rectangle(
            (tx, ty), tw, th, facecolor=PALETTE[
                {"Engineering": "a", "Sales": "b", "Marketing": "c"}[div]],
            alpha=.85, edgecolor="white", linewidth=2))
        # Only label a rectangle when the text FITS. Otherwise, tooltip/table.
        if tw > 1.5 and th > 0.8:
            axes[1].text(tx + tw/2, ty + th/2, f"{tr['team']}\n{tr['headcount']}",
                         ha="center", va="center", fontsize=6.5, color="white")
axes[1].set_xlim(0, 10); axes[1].set_ylim(0, 6); axes[1].axis("off")
axes[1].set_title("treemap: area = headcount, full space use",
                  loc="left", fontsize=9)

# 3. Icicle: rectangles laid out linearly, easier to read than a sunburst.
y0 = 0
for div, total in div_totals.items():
    axes[2].add_patch(mpatches.Rectangle(
        (0, y0), 1, total, facecolor=PALETTE[
            {"Engineering": "a", "Sales": "b", "Marketing": "c"}[div]],
        alpha=.5, edgecolor="white", linewidth=2))
    axes[2].text(.5, y0 + total/2, div, ha="center", va="center", fontsize=7)
    ty = y0
    for _, tr in tree[tree["division"] == div].iterrows():
        axes[2].add_patch(mpatches.Rectangle(
            (1.05, ty), 1, tr["headcount"], facecolor=PALETTE[
                {"Engineering": "a", "Sales": "b", "Marketing": "c"}[div]],
            alpha=.9, edgecolor="white", linewidth=2))
        axes[2].text(1.55, ty + tr["headcount"]/2,
                     f"{tr['team']} ({tr['headcount']})", ha="center",
                     va="center", fontsize=6, color="white")
        ty += tr["headcount"]
    y0 += total
axes[2].set_xlim(-.05, 2.2); axes[2].set_ylim(0, y0); axes[2].axis("off")
axes[2].set_title("icicle: length = headcount, readable labels",
                  loc="left", fontsize=9)
fig.savefig("/tmp/claude-0/hierarchy.png", dpi=110, bbox_inches="tight")
plt.close(fig)

print("""
  HIERARCHY FORMS AND THEIR ENCODING

  form         encodes value by   space use   label room   reads
  ---------------------------------------------------------------------------
  node-link    nothing            poor        good         structure clearly
  treemap      AREA               excellent   poor, small  coarse proportions
  sunburst     ANGLE + area       good        worst        worse than treemap
  icicle       LENGTH             good        good         best of the three
  dendrogram   merge height       fair        fair         cluster structure
  Sankey       flow width         fair        fair         stage-to-stage flow

  Prefer icicle over sunburst: length beats angle, and labels fit.
""")

  WHY CENTRALITY BEFORE DRAWING

  node   degree   betweenness   what it is
  ---------------------------------------------------------
  11     11       0.18          a hub inside its community
  99      3       0.58          a BRIDGE between communities
  7       9       0.03          well connected, locally only

    Degree ranks node 99 near the bottom. Removing it splits the graph
    into three components. Betweenness is what found it.

  THE HAIRBALL THRESHOLD (average degree 6)

     50 nodes ->    150 edges   readable
    200 nodes ->    600 edges   at the limit
   1000 nodes ->  3,000 edges   hairball
   5000 nodes -> 15,000 edges   solid ink

  ALTERNATIVES ABOVE THE THRESHOLD

  node-link (hairball)      adjacency matrix        aggregated
     \ | / \ | /              ###....              (a)---3---(b)
    --#####-----              ###....                \       /
     / | \ / | \              ...##..                 4     1
    unreadable                ...##..                   \ /
                              ......#                   (c)
                              clusters = diagonal      the question was
                              BLOCKS                   about groups anyway

  LAYOUT DISTANCE IS NOT A QUANTITY

    Two nodes drawn 3cm apart and two drawn 6cm apart do NOT differ by a
    factor of two in anything. Force layouts show topology, not distance,
    and they change on every unseeded run.

Key Takeaways

• A graph layout invents node positions, so distance in a force-directed drawing measures nothing and readers over-interpret it.
• Force-directed layouts are non-deterministic, so fix the seed for any figure you publish or compare.
• Compute centrality before drawing, because a low-degree, high-betweenness node is a bridge that degree alone would rank as unimportant.
• Past a few hundred nodes a node-link diagram becomes a hairball; switch to an adjacency matrix, aggregate to communities, or filter.
• For hierarchies, icicle plots encode with length and leave room for labels, which beats the treemap's area and the sunburst's angle.

🧪 Practice

1. Build a graph with three communities joined by one bridge node, compute degree and betweenness, and show that the two rank the bridge differently.
2. Remove the highest-betweenness node and report the change in the number of connected components.
3. Render the same graph with two different layout seeds and quantify how far the nodes moved.
4. Display a graph as a node-link diagram and as a community-ordered adjacency matrix, and describe what each makes obvious.
5. Interview: You must show a 50,000-node transaction network to an investigator. What do you actually build? (Hint: what is the investigator's unit of interest, and what can a drawing of 50,000 nodes convey?)

Small Multiples and Faceting

Small multiples, also called faceting or trellis plots, repeat one chart across levels of a variable. Edward Tufte named the technique and it deserves its reputation as the most reliable multivariate display available, for a reason that follows directly from the perceptual ranking discussed earlier: comparison across panels uses position, the channel humans read most accurately, while squeezing another variable into one panel forces the reader onto colour, size, or shape.

The practical comparison is worth being concrete about. Eight series overlaid in one chart requires eight distinguishable colours, produces crossing lines, needs a legend the reader must keep consulting, and buries whichever series matters. The same eight series as eight small panels needs one colour, has no occlusion, labels each panel directly, and lets the reader scan for the shape that differs. The overlay wins only when the comparison is between adjacent series at the same x position, which is precisely when crossing matters.

Four rules make faceting work, and violating any of them undermines the whole technique.

Share the scales. This is not a preference. A grid layout invites the reader to compare panels, so independent axes make the layout itself lie: a panel whose values are a tenth of another's looks identical when each is scaled to fit. Share the axes, and if one panel's range genuinely dwarfs the others, that is a finding to state rather than a problem to normalise away.

Order the panels meaningfully. Alphabetical order is almost never the right order. Sort by the value being shown, by magnitude, by time, or by geography. Ordering is free information and the reader will read meaning into the sequence whether you intended one or not.

Keep panels small and chrome minimal. The point is to see many at once. Label each panel with a short title, draw axis labels once on the outer edges rather than on every panel, and drop gridlines to hairlines.

Add context to each panel. A faint copy of the overall average, or of all series in grey behind the highlighted one, lets the reader see both the specific and the general at once. This is the single highest-value addition to a small-multiple grid.

The variants: facet wrap flows panels into a grid by one variable, which is the usual case. Facet grid uses two variables for rows and columns, showing their interaction. A trellis across time periods shows how a relationship itself changes.

Faceting's limit is panel count. Beyond roughly twenty panels each becomes too small to read, and the answer is to aggregate the tail, select the interesting subset, or switch to a heatmap, which is effectively the fully compressed limit of a small-multiple grid with colour replacing position.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import seaborn as sns

BLUE, GRAY, INK, MUTED = "#2a78d6", "#c9c9c5", "#0b0b0b", "#52514e"
sns.set_theme(style="whitegrid", rc={"grid.linewidth": .5})
rng = np.random.default_rng(6)

# Nine regions at very different scales; ONE of them is the story.
months = pd.date_range("2024-01-01", periods=18, freq="MS")
scales = {"north": 100, "south": 420, "east": 38, "west": 260, "central": 95,
          "alpine": 12, "coastal": 180, "delta": 55, "plateau": 8}
rows = []
for region, scale in scales.items():
    trend = -1.6 if region == "coastal" else rng.uniform(0.2, 0.9)  # the story
    series = scale + np.cumsum(rng.normal(trend, scale * 0.04, 18))
    rows += [{"region": region, "month": m, "revenue": round(v, 1)}
             for m, v in zip(months, series)]
df = pd.DataFrame(rows)

# --- OVERLAY VS FACET: the same nine series, two treatments. ---
fig = plt.figure(figsize=(14, 4))
ax = fig.add_subplot(1, 2, 1)
for region, g in df.groupby("region"):
    ax.plot(g["month"], g["revenue"], lw=1.6, label=region)
ax.legend(frameon=False, fontsize=6, ncol=3)
ax.set_title("9 series overlaid: 9 hues, crossings, nothing stands out",
             loc="left", fontsize=9)
ax.set_ylabel("revenue")
fig.savefig("/tmp/claude-0/overlay.png", dpi=110, bbox_inches="tight")
plt.close(fig)

# --- RULE 1: SHARED SCALES. Compare the two versions directly. ---
order = (df.groupby("region")["revenue"].mean()
           .sort_values(ascending=False).index)      # RULE 2: meaningful order

for shared, tag in [(False, "independent"), (True, "shared")]:
    fig, axes = plt.subplots(3, 3, figsize=(10, 6.4),
                             sharex=True, sharey=shared)
    for ax, region in zip(axes.ravel(), order):
        g = df[df["region"] == region]
        ax.plot(g["month"], g["revenue"], lw=1.8, color=BLUE)
        ax.set_title(region, loc="left", fontsize=8, color=INK)
        ax.tick_params(labelsize=6)
        ax.spines[["top", "right"]].set_visible(False)
    fig.suptitle(f"{tag} y scales", x=.02, ha="left", fontsize=10)
    fig.savefig(f"/tmp/claude-0/facet_{tag}.png", dpi=110, bbox_inches="tight")
    plt.close(fig)

span = df.groupby("region")["revenue"].agg(["min", "max"])
print("panel ranges (why independent scales lie):")
print(span.round(0).sort_values("max", ascending=False).to_string())
print(f"\nlargest panel max {span['max'].max():.0f},"
      f" smallest {span['max'].min():.0f}"
      f" -> a {span['max'].max()/span['max'].min():.0f}x difference that")
print("independent scales render as nine identically sized charts.")

# --- RULE 4: CONTEXT IN EVERY PANEL. The highest-value addition. ---
fig, axes = plt.subplots(3, 3, figsize=(10.5, 6.6), sharex=True, sharey=True)
overall = df.groupby("month")["revenue"].mean()
for ax, region in zip(axes.ravel(), order):
    # All other series, faint: the reader sees the specific AND the general.
    for other, g2 in df.groupby("region"):
        if other != region:
            ax.plot(g2["month"], g2["revenue"], lw=.9, color="#e4e4e1", zorder=1)
    g = df[df["region"] == region]
    is_story = region == "coastal"
    ax.plot(g["month"], g["revenue"], lw=2,
            color="#e34948" if is_story else BLUE, zorder=3)
    # Direct-label the endpoint of the panel's own series only.
    ax.annotate(f"{g['revenue'].iloc[-1]:.0f}",
                (g["month"].iloc[-1], g["revenue"].iloc[-1]),
                xytext=(3, 0), textcoords="offset points", fontsize=6.5,
                color="#e34948" if is_story else MUTED, va="center")
    change = g["revenue"].iloc[-1] / g["revenue"].iloc[0] - 1
    ax.set_title(f"{region}   {change:+.0%}", loc="left", fontsize=8,
                 color="#e34948" if is_story else INK,
                 fontweight="bold" if is_story else "normal")
    ax.tick_params(labelsize=6)
    ax.spines[["top", "right"]].set_visible(False)
fig.suptitle("Coastal is the only region declining", x=.02, ha="left",
             fontsize=11, fontweight="bold")
fig.text(.02, .935, "Monthly revenue by region, shared scale. Other regions "
         "shown faint for context.", fontsize=8, color=MUTED)
fig.savefig("/tmp/claude-0/facet_context.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("\nThe declining region is now findable in under a second, because the")
print("grid encodes on POSITION and the one exception is coloured.")

# --- FACET GRID: two variables, rows and columns, showing interaction. ---
n = 2400
panel = pd.DataFrame({
    "segment": rng.choice(["smb", "mid", "enterprise"], n),
    "channel": rng.choice(["direct", "partner"], n),
    "seats": rng.integers(1, 250, n),
})
# The seats-to-revenue slope differs BY CELL: that interaction is the point.
slope = {("smb", "direct"): 9, ("smb", "partner"): 5,
         ("mid", "direct"): 14, ("mid", "partner"): 13,
         ("enterprise", "direct"): 26, ("enterprise", "partner"): 11}
panel["revenue"] = [
    slope[(s, c)] * k + rng.normal(0, 260)
    for s, c, k in zip(panel["segment"], panel["channel"], panel["seats"])]

g = sns.FacetGrid(panel, row="channel", col="segment", height=2.3,
                  aspect=1.25, margin_titles=True,
                  col_order=["smb", "mid", "enterprise"])   # meaningful order
g.map_dataframe(sns.regplot, x="seats", y="revenue",
                scatter_kws={"s": 6, "alpha": .35, "color": BLUE},
                line_kws={"color": "#eb6834", "lw": 2})
g.set_titles(col_template="{col_name}", row_template="{row_name}")
g.figure.suptitle("Revenue per seat rises with segment, but only via direct",
                  x=.02, ha="left", fontsize=11)
g.figure.subplots_adjust(top=.86)
g.figure.savefig("/tmp/claude-0/facet_grid.png", dpi=110, bbox_inches="tight")
plt.close(g.figure)

fitted = (panel.groupby(["channel", "segment"])
          .apply(lambda d: np.polyfit(d["seats"], d["revenue"], 1)[0],
                 include_groups=False)
          .unstack().round(1))
print("\nfitted revenue per seat by cell (the interaction the grid shows):")
print(fitted.to_string())

# --- THE PANEL-COUNT LIMIT, and the heatmap as the compressed limit. ---
print("""
  PANEL COUNT

    1-4      large panels, full detail, axis labels on each
    5-12     the sweet spot; shared axes, outer labels only
    13-25    small panels; titles only, no tick labels inside
    25-60    too small to read -> aggregate the tail, or use a heatmap
    60+      heatmap, or select the subset that matters
""")

# A heatmap IS a small-multiple grid with colour replacing position.
wide = df.pivot(index="region", columns="month", values="revenue")
indexed = wide.div(wide.iloc[:, 0], axis=0) * 100      # index to 100 at t0
fig, ax = plt.subplots(figsize=(11, 3.4))
sns.heatmap(indexed.loc[order], cmap="RdBu_r", center=100, ax=ax,
            cbar_kws={"label": "indexed to 100 at Jan 2024"},
            xticklabels=[m.strftime("%b") for m in wide.columns],
            linewidths=2, linecolor="white")           # 2px surface gap
ax.set_title("the same nine series as a heatmap: 60+ panels would fit here",
             loc="left", fontsize=9)
ax.set_ylabel("")
fig.savefig("/tmp/claude-0/heatmap_limit.png", dpi=110, bbox_inches="tight")
plt.close(fig)
print("Indexing each row to 100 at t0 is what makes a heatmap of series at")
print("different scales comparable; the diverging scale is centred on 100.")

  OVERLAY VS SMALL MULTIPLES, 9 SERIES

  OVERLAY                          SMALL MULTIPLES
   ___/\__/\___                    +------+------+------+
   /\_/  \/  \/                    | /\   | __/  | \__  |
  _/\/\__/\_/\_                    +------+------+------+
   needs 9 hues                    | _/\  | /\/  | \/\  |
   lines cross and occlude         +------+------+------+
   legend lookup per series        | /~\  | \    | ___  |
   the one decline is invisible    +------+------+------+
                                   1 hue, no occlusion, titles in place,
                                   and the odd shape is found by scanning

  Overlay wins ONLY when the comparison is between series at the same x.

  WHY SHARED SCALES ARE NOT OPTIONAL

  INDEPENDENT                       SHARED
  north  south  alpine             north  south  alpine
   /\     /\     /\                 /\     /\
   ||     ||     ||                 ||           .
   ||     ||     ||                 ||    __      (alpine is 12 vs south 420)
  all look the same size           the 35x difference is visible

  A grid invites comparison. Independent scales make the LAYOUT lie.

  THE FOUR RULES

    1. SHARE the scales                  (or the grid misleads)
    2. ORDER panels by value/time/size   (never alphabetically)
    3. MINIMISE chrome                   (outer axis labels only)
    4. ADD context to each panel         (faint others, or the average)

  THE COMPRESSION LADDER

    one panel per series       -> small multiples   position, high precision
    one ROW per series         -> heatmap           colour, low precision
                                                    but 60+ series fit

Key Takeaways

• Small multiples compare using position, the most accurately read channel, while overlays force the reader onto colour and cope with occlusion.
• Shared scales are mandatory, because a grid invites comparison and independent axes make the layout itself misleading.
• Order panels by value, magnitude, time, or geography, since readers infer meaning from the sequence regardless of your intent.
• Add a faint copy of the other series or the overall average to each panel so the specific and the general read together.
• Past roughly twenty panels, aggregate the tail or move to a heatmap, which is the same grid with colour replacing position.

🧪 Practice

1. Take nine series at widely different scales, render them overlaid and as a faceted grid, and time how long it takes someone to find the one that declines.
2. Build the same grid twice, with independent and shared y axes, and state the wrong conclusion the first invites.
3. Add faint context series and an endpoint label to each panel of a grid, and sort the panels by their mean value.
4. Construct a two-variable facet grid where the relationship differs by cell, then print the fitted slope per cell.
5. Interview: You need to show a metric's trend for 80 stores to a regional manager. What do you build? (Hint: 80 panels is unreadable, so what does the manager actually need to identify, and which form fits 80 rows?)

<a id="74-dashboards-and-reporting"></a>

7.4 Dashboards and Reporting

A chart that nobody acts on has no value, so the last step is delivery: the layouts, definitions, tooling, automation, and narrative structure that carry a finding from your notebook to a decision.

Dashboard Design Patterns

A dashboard is not a collection of charts, it is an interface built for a recurring question. That distinction explains most dashboard failure: teams build a place to put charts, the collection grows, nobody can say what any of it is for, and eventually everyone goes back to asking an analyst.

Start by classifying what you are building, because the three kinds have different designs and mixing them produces something that serves none of them.

Operational dashboards answer "is anything wrong right now?" They are watched frequently, refresh often, and are dominated by current state against thresholds. The reader wants to spend zero seconds on them when things are fine, so the design goal is that normal looks obviously normal.

Analytical dashboards answer "why did this change?" They are used occasionally, in longer sessions, and need breakdowns, comparisons, and enough interaction to follow a question down a level.

Strategic dashboards answer "are we on track?" They are reviewed on a cadence, show a few metrics against targets over long horizons, and should be the sparsest of the three.

The layout patterns that work follow from how people read a screen. Attention goes to the top left first, so the most important number belongs there, and importance should decrease down and to the right. The standard structure is an inverted pyramid: a row of headline numbers at the top, the main trend chart below, breakdowns below that, and detail tables last. A reader who stops after the first row should still have the answer.

Five rules carry most of the value.

One filter row, above everything it scopes. All charts re-render against the same slice, so the numbers always agree. Per-chart filters are how a dashboard starts contradicting itself, and a chart with its own date range belongs on a different page.

Every number needs a comparison. "Revenue: 1.2M" is not information. Against last month, against target, or against the same period last year, it becomes information. A number with no reference point cannot be acted on.

Show uncertainty and sample size. A conversion rate of 8% from 25 sessions should not look like one from 25,000. Cells built on small samples need their counts visible or they will be read as signal.

Define everything in reach. Every metric needs a definition available from the dashboard itself, ideally on hover or behind a single click. Ambiguity about what a number means is the leading cause of dashboards being distrusted and abandoned.

Handle empty, loading, and error states deliberately. A chart that renders blank when a filter matches nothing looks broken. Say "no data for this selection". While refreshing, hold the previous render at reduced opacity rather than flashing a skeleton, so the layout does not jump.

The anti-patterns are worth naming because each is common. The wall of charts with no hierarchy makes the reader find the story. Vanity metrics such as cumulative totals only ever go up and cannot indicate a problem. Overloading one view with every possible breakdown serves nobody. And a dashboard with no owner rots: definitions drift, a data source changes, and nobody notices until someone makes a bad decision from it.

import numpy as np
import pandas as pd

rng = np.random.default_rng(12)

# --- A METRIC WITH ITS COMPARISONS: what makes a number actionable. ---
days = pd.date_range("2024-04-01", periods=90, freq="D")
sessions = (9000 + 600 * np.sin(np.arange(90) / 7 * 2 * np.pi)
            + np.arange(90) * 18 + rng.normal(0, 350, 90)).round(0)
conv = np.clip(0.041 + rng.normal(0, .002, 90) - (np.arange(90) > 70) * .006,
               0, 1)
daily = pd.DataFrame({"date": days, "sessions": sessions,
                      "conversions": (sessions * conv).round(0)})
daily["conv_rate"] = daily["conversions"] / daily["sessions"]

def stat_tile(df: pd.DataFrame, value_col: str, target: float | None = None,
              window: int = 7) -> dict:
    """A headline number is never shipped alone. Attach the comparisons."""
    cur = df[value_col].tail(window).mean()
    prev = df[value_col].tail(window * 2).head(window).mean()
    yoy_proxy = df[value_col].head(window).mean()      # earliest window
    # A standard error, so a small-sample tile does not read as certain.
    n = df["sessions"].tail(window).sum()
    se = np.sqrt(cur * (1 - cur) / n) if 0 < cur < 1 else np.nan
    out = {
        "value": round(cur, 4),
        "vs_prev_period": round((cur / prev - 1) * 100, 1),
        "vs_period_start": round((cur / yoy_proxy - 1) * 100, 1),
        "sample_n": int(n),
        "std_err": round(se, 5) if se == se else None,
    }
    if target is not None:
        out["vs_target_pts"] = round((cur - target) * 100, 2)
        out["status"] = ("good" if cur >= target
                         else "warning" if cur >= target * .9 else "critical")
    return out

tile = stat_tile(daily, "conv_rate", target=0.040)
print("HEADLINE TILE")
for k, v in tile.items():
    print(f"  {k:16s} {v}")
print("\n'4.1%' alone is not information. Against last week, against the")
print("target, and with the sample size, it supports a decision.")

# --- THE INVERTED PYRAMID, as an explicit layout spec. ---
LAYOUT = [
    {"tier": 1, "row": "KPI row",      "height": "small",
     "content": "4-6 stat tiles: value, delta, sparkline, status"},
    {"tier": 2, "row": "primary trend", "height": "large",
     "content": "the one metric that defines success, over time, vs target"},
    {"tier": 3, "row": "breakdowns",    "height": "medium",
     "content": "2-3 charts: the same metric split by the key dimensions"},
    {"tier": 4, "row": "detail table",  "height": "medium",
     "content": "sortable rows, the accessible twin of every chart above"},
]
print("\n", pd.DataFrame(LAYOUT).to_string(index=False))
print("\nA reader who stops after tier 1 must still have the answer.")

# --- SMALL SAMPLES: the cell that will be misread. ---
segments = pd.DataFrame({
    "segment": ["organic", "paid", "email", "partner", "affiliate"],
    "sessions": [42_000, 18_500, 6_200, 410, 38],
    "conversions": [1_780, 690, 300, 24, 4],
})
segments["rate"] = segments["conversions"] / segments["sessions"]
# Wilson interval: behaves correctly at small n, unlike the normal approximation.
def wilson(k, n, z=1.96):
    if n == 0:
        return (np.nan, np.nan)
    p = k / n
    denom = 1 + z**2 / n
    centre = (p + z**2 / (2*n)) / denom
    half = z * np.sqrt(p*(1-p)/n + z**2/(4*n**2)) / denom
    return (centre - half, centre + half)

ci = segments.apply(lambda r: wilson(r["conversions"], r["sessions"]), axis=1)
segments["ci_lo"] = [c[0] for c in ci]
segments["ci_hi"] = [c[1] for c in ci]
segments["ci_width_pts"] = ((segments["ci_hi"] - segments["ci_lo"]) * 100).round(2)
segments["reportable"] = segments["sessions"] >= 1000
print("\n", segments.assign(rate=(segments["rate"]*100).round(2))
        [["segment", "sessions", "rate", "ci_width_pts", "reportable"]]
        .to_string(index=False))
print("\nAffiliate shows the highest rate and its interval spans ~40 points.")
print("Without the count or the interval, a reader reallocates budget to it.")

# --- ONE FILTER ROW: why per-chart filters break a dashboard. ---
print("""
  FILTER SCOPE

  RIGHT                              WRONG
  +-------------------------------+   +------------------+ +---------------+
  | date range | region | segment |   | chart A [last 7d]| | chart B [MTD] |
  +-------------------------------+   +------------------+ +---------------+
  | every chart below uses THIS   |   two charts, two ranges, and the
  | slice, so totals reconcile    |   totals never add up. Readers lose
  +-------------------------------+   trust in the whole page.

  If one chart genuinely needs its own range, it belongs on another page.
""")

# --- STATE HANDLING: empty, loading, error, stale. ---
def render_state(rows: int, loading: bool, error: str | None,
                 as_of: pd.Timestamp, now: pd.Timestamp) -> str:
    if error:
        return f"ERROR: {error} -- last good data {as_of:%Y-%m-%d %H:%M} UTC"
    if loading:
        return "REFRESHING: hold previous render at 60% opacity, no layout jump"
    if rows == 0:
        return "EMPTY: 'No data for this selection' + a clear-filters button"
    staleness = (now - as_of).total_seconds() / 3600
    if staleness > 24:
        return f"STALE: data is {staleness:.0f}h old -- show a warning banner"
    return f"OK: {rows:,} rows, as of {as_of:%H:%M} UTC"

now = pd.Timestamp("2024-06-30 12:00")
for case in [
        {"rows": 4210, "loading": False, "error": None,
         "as_of": now - pd.Timedelta("2h")},
        {"rows": 0, "loading": False, "error": None,
         "as_of": now - pd.Timedelta("1h")},
        {"rows": 900, "loading": True, "error": None,
         "as_of": now - pd.Timedelta("1h")},
        {"rows": 4210, "loading": False, "error": None,
         "as_of": now - pd.Timedelta("40h")},
        {"rows": 0, "loading": False, "error": "upstream table missing",
         "as_of": now - pd.Timedelta("26h")}]:
    print(f"  {render_state(now=now, **case)}")

# --- VANITY VS ACTIONABLE METRICS ---
comparison = pd.DataFrame([
    ("cumulative signups",   "vanity",     "only ever rises; cannot signal a problem",
     "weekly new signups"),
    ("total page views",     "vanity",     "no link to any decision",
     "conversion rate by landing page"),
    ("registered users",     "vanity",     "includes everyone who ever signed up",
     "weekly active users"),
    ("average revenue",      "misleading", "one whale moves it",
     "median revenue + the distribution"),
    ("uptime this quarter",  "lagging",    "too coarse to act on",
     "error rate per minute vs threshold"),
], columns=["metric", "problem", "why", "use instead"])
print("\n", comparison.to_string(index=False))

# --- A DASHBOARD SPEC AUDIT ---
def audit_dashboard(spec: dict) -> list[str]:
    issues = []
    if spec.get("kind") not in {"operational", "analytical", "strategic"}:
        issues.append("no declared kind: it will drift into serving none well")
    if spec.get("per_chart_filters"):
        issues.append("per-chart filters: totals will stop reconciling")
    if not spec.get("comparisons_on_every_metric"):
        issues.append("bare numbers with no baseline, target, or prior period")
    if not spec.get("shows_sample_size"):
        issues.append("no sample size or interval: small cells read as signal")
    if not spec.get("definitions_reachable"):
        issues.append("metric definitions not reachable from the dashboard")
    if not spec.get("empty_state"):
        issues.append("no empty state: a matched-nothing filter looks broken")
    if not spec.get("owner"):
        issues.append("no named owner: definitions will drift unnoticed")
    if not spec.get("as_of_stamp"):
        issues.append("no as-of timestamp: staleness is invisible")
    if spec.get("n_charts", 0) > 9:
        issues.append(f"{spec['n_charts']} charts: no hierarchy, reader must hunt")
    if spec.get("kind") == "operational" and not spec.get("thresholds"):
        issues.append("operational view without thresholds: 'normal' is not obvious")
    return issues

print("\naudit of a typical grown-organically dashboard:")
for i in audit_dashboard({"n_charts": 17, "per_chart_filters": True}):
    print("  -", i)
print("\naudit of a designed one:",
      audit_dashboard({
          "kind": "operational", "n_charts": 6, "thresholds": True,
          "comparisons_on_every_metric": True, "shows_sample_size": True,
          "definitions_reachable": True, "empty_state": True,
          "owner": "growth-analytics", "as_of_stamp": True}) or "no issues")

  THE THREE KINDS

  kind         question              cadence      design goal
  ---------------------------------------------------------------------------
  operational  anything wrong NOW?   continuous   normal looks OBVIOUSLY normal
  analytical   why did this change?  occasional   breakdowns + drill-down
  strategic    are we on track?      weekly/mo.   few metrics, long horizon,
                                                  against target
  Mixing them produces a page that serves none of them.

  THE INVERTED PYRAMID

  +---------------------------------------------------------------+
  | date range | region | segment            <- ONE filter row    |
  +---------------------------------------------------------------+
  | [ 4.1% ]  [ 9,240 ]  [ $1.2M ]  [ 2.3d ]  <- tier 1: KPIs     |
  |  +0.2pt    -3.1%      +8.4%      flat        value + delta     |
  +---------------------------------------------------------------+
  |                                                               |
  |     tier 2: the ONE trend that defines success, vs target     |
  |                                                               |
  +----------------------------+----------------------------------+
  |  tier 3: by channel        |  tier 3: by segment              |
  +----------------------------+----------------------------------+
  |  tier 4: detail table (the accessible twin of the charts)     |
  +---------------------------------------------------------------+
   Attention flows top-left -> down-right. A reader who stops after
   tier 1 still has the answer.

  EVERY NUMBER NEEDS A COMPARISON

    "Revenue: $1.2M"                          not information
    "Revenue: $1.2M, +8% vs last month,       information: a reader can
     3% below the $1.24M target, n=18,402"    act on this

  STATE HANDLING

    OK          rows + as-of timestamp
    EMPTY       "No data for this selection" + clear-filters (never blank)
    REFRESHING  hold the old render at reduced opacity (never a skeleton flash)
    STALE       a banner naming the age of the data
    ERROR       what broke + when the last good data was

Key Takeaways

• Classify the dashboard as operational, analytical, or strategic, because the three have different designs and mixing them serves nobody.
• Lay it out as an inverted pyramid so a reader who stops after the top row still has the answer.
• One filter row scoping everything keeps the numbers reconciling; per-chart filters are how a dashboard contradicts itself.
• Every number needs a comparison against a target, a prior period, or a baseline, plus a visible sample size.
• Design the empty, loading, stale, and error states, and give the dashboard a named owner or its definitions will drift unnoticed.

🧪 Practice

1. Write a stat tile function returning a value alongside its change versus the prior period, its gap to target, and its sample size.
2. Compute Wilson confidence intervals for conversion rates across segments of very different sizes and identify which are not reportable.
3. Take a dashboard you use and classify it as operational, analytical, or strategic, then list the charts that do not fit that classification.
4. Implement a dashboard audit function and run it against two real dashboards.
5. Interview: A dashboard has 22 charts and the team says it is not useful. Where do you start? (Hint: what single question should it answer, and what would you delete to make that question answerable in five seconds?)

Metric and KPI Definition

Most disagreements about numbers are not disagreements about data, they are disagreements about definitions. Two teams report different revenue because one counts refunds and the other does not. Both are right and the meeting is lost to reconciliation. Defining metrics precisely is unglamorous and it prevents more wasted time than any other practice in analytics.

A metric definition needs five things, and a definition missing any of them will eventually produce two numbers.

The precise event or entity. Not "active users" but "distinct user IDs with at least one recorded session of 10 seconds or more". Every vague word is a place two implementations will diverge.

The time basis. Calendar month or rolling 30 days? Event timestamp or ingestion timestamp? In which time zone? A daily metric computed in UTC and one computed in local time disagree for a predictable fraction of users, forever.

The population and its exclusions. Internal accounts, bots, test users, churned accounts, refunded orders. Exclusions are where definitions quietly differ most.

The aggregation. Sum, mean, median, or percentile; and critically, whether a ratio is computed per unit and then averaged, or as a total over a total. The average of per-country conversion rates is not the overall conversion rate, as covered in chapter 6, and this single distinction accounts for a large share of reconciliation disputes.

The grain. One row per what? A metric whose grain is unstated cannot be joined to anything safely.

Beyond individual definitions, three structural ideas make a metric set work.

Distinguish leading from lagging indicators. Revenue is lagging: it tells you what happened and it is too late to act on. Trial-to-paid conversion is leading: it moves first and predicts revenue. A dashboard of only lagging metrics is a history lesson.

Choose input metrics you can actually move. Output metrics such as revenue are the goal; input metrics such as page load time, activation rate, or support response time are what a team can act on directly. Teams should be measured on outputs and given inputs to work with.

Pair every primary metric with a guardrail. Optimising one number in isolation reliably damages something else. Session count paired with session quality. Speed of shipping paired with incident rate. Conversion paired with refund rate. This is the practical answer to Goodhart's law, that a measure used as a target stops being a good measure, and the mechanism is always the same: someone finds a way to move the number without producing the underlying value.

The engineering practice that makes all of this stick is a metrics layer: a single place where each metric is defined once, in code, and every dashboard and query derives from it. The alternative, each dashboard implementing its own SQL, guarantees drift.

import numpy as np
import pandas as pd
from dataclasses import dataclass, field

rng = np.random.default_rng(15)

# --- HOW MANY NUMBERS CAN "ACTIVE USERS" PRODUCE? ---
n = 6000
events = pd.DataFrame({
    "user": rng.integers(1, 1200, n),
    "ts_utc": pd.Timestamp("2024-06-01") + pd.to_timedelta(
        rng.integers(0, 30 * 24 * 60, n), "m"),
    "duration_s": rng.gamma(2, 20, n).round(0),
    "is_internal": rng.random(n) < 0.04,
    "is_bot": rng.random(n) < 0.07,
})
events["tz_offset_h"] = rng.choice([-8, -5, 0, 1, 9], n)
events["ts_local"] = events["ts_utc"] + pd.to_timedelta(
    events["tz_offset_h"], "h")

definitions = {
    "any event, UTC month":
        lambda d: d[d.ts_utc.dt.month == 6]["user"].nunique(),
    "+ exclude internal":
        lambda d: d[(d.ts_utc.dt.month == 6) & ~d.is_internal]["user"].nunique(),
    "+ exclude bots":
        lambda d: d[(d.ts_utc.dt.month == 6) & ~d.is_internal
                    & ~d.is_bot]["user"].nunique(),
    "+ session >= 10s":
        lambda d: d[(d.ts_utc.dt.month == 6) & ~d.is_internal & ~d.is_bot
                    & (d.duration_s >= 10)]["user"].nunique(),
    "same, but LOCAL month":
        lambda d: d[(d.ts_local.dt.month == 6) & ~d.is_internal & ~d.is_bot
                    & (d.duration_s >= 10)]["user"].nunique(),
    "rolling 30d to Jun 30":
        lambda d: d[(d.ts_utc > pd.Timestamp("2024-05-31")) & ~d.is_internal
                    & ~d.is_bot & (d.duration_s >= 10)]["user"].nunique(),
}
print("ONE PHRASE, SIX NUMBERS:")
results = {k: f(events) for k, f in definitions.items()}
for k, v in results.items():
    print(f"  {v:>5}   {k}")
spread = max(results.values()) - min(results.values())
print(f"\nspread: {spread} users ({spread/max(results.values()):.0%})."
      " Every one of these is a defensible reading of 'active users'.")

# --- RATIO OF SUMS VS MEAN OF RATIOS: the classic reconciliation dispute. ---
by_country = pd.DataFrame({
    "country": ["US", "CA", "DE", "JP"],
    "sessions": [80_000, 9_000, 6_000, 500],
    "conversions": [3_200, 450, 240, 45],
})
by_country["rate"] = by_country["conversions"] / by_country["sessions"]
mean_of_ratios = by_country["rate"].mean()
ratio_of_sums = by_country["conversions"].sum() / by_country["sessions"].sum()
print(f"\nmean of per-country rates : {mean_of_ratios:.4%}")
print(f"total conv / total sessions: {ratio_of_sums:.4%}")
print(f"difference: {(mean_of_ratios - ratio_of_sums)*100:.2f} points."
      " Both are computable; only one is 'the conversion rate'.")
print("Specify which in the definition, or two dashboards will disagree.")

# --- A METRIC DEFINITION AS CODE, not as a sentence in a wiki. ---
@dataclass
class MetricDefinition:
    name: str
    description: str
    entity: str                       # one row per WHAT
    event_filter: str                 # the precise event
    time_basis: str                   # calendar/rolling, which timestamp, which tz
    exclusions: list[str] = field(default_factory=list)
    aggregation: str = "count_distinct"
    ratio_rule: str | None = None     # ratio-of-sums or mean-of-ratios
    owner: str = "unowned"
    guardrails: list[str] = field(default_factory=list)
    source_table: str = ""

    def validate(self) -> list[str]:
        problems = []
        for word in ("active", "engaged", "successful", "valid", "recent"):
            described = word in self.description.lower()
            if described and word not in self.event_filter.lower():
                problems.append(f"'{word}' used in the description but not defined"
                                " in the filter")
        if "tz" not in self.time_basis and "UTC" not in self.time_basis:
            problems.append("time basis does not state a time zone")
        if self.aggregation.startswith("ratio") and not self.ratio_rule:
            problems.append("ratio metric without a ratio_rule: ambiguous")
        if not self.exclusions:
            problems.append("no exclusions listed: are bots and internal in?")
        if self.owner == "unowned":
            problems.append("no owner: this definition will drift")
        if not self.guardrails:
            problems.append("no guardrail metric: Goodhart's law applies")
        return problems

wau = MetricDefinition(
    name="weekly_active_users",
    description="Distinct users with a qualifying session in the trailing 7 days",
    entity="one row per (week_start, user_id)",
    event_filter="session_duration_s >= 10 AND event_name = 'session_start'",
    time_basis="rolling 7 days ending Sunday 23:59, event timestamp, tz=UTC",
    exclusions=["is_internal = true", "is_bot = true", "account_state = 'test'"],
    aggregation="count_distinct(user_id)",
    owner="growth-analytics",
    guardrails=["session_quality_score", "support_ticket_rate"],
    source_table="analytics.fct_sessions",
)
print(f"\n{wau.name}: {wau.validate() or 'definition complete'}")

vague = MetricDefinition(
    name="engagement",
    description="How engaged our active users are",
    entity="user",
    event_filter="event_name IS NOT NULL",
    time_basis="monthly",
    aggregation="ratio",
)
print(f"\n{vague.name} problems:")
for p in vague.validate():
    print("  -", p)

# --- LEADING VS LAGGING, AND INPUT VS OUTPUT ---
metrics = pd.DataFrame([
    ("monthly revenue",        "lagging", "output", "cannot be moved directly"),
    ("trial-to-paid rate",     "leading", "output", "predicts revenue ~1 mo ahead"),
    ("activation within 7d",   "leading", "input",  "a team can change onboarding"),
    ("page load p95",          "leading", "input",  "an engineer can fix this today"),
    ("net revenue retention",  "lagging", "output", "the quarterly scoreboard"),
    ("support first-response", "leading", "input",  "staffing and tooling move it"),
], columns=["metric", "timing", "type", "note"])
print("\n", metrics.to_string(index=False))
print("\nA dashboard of only lagging outputs is a history lesson. Pair each")
print("output with the inputs a team can actually act on.")

# --- GUARDRAILS: what optimising one number breaks. ---
print("""
  GOODHART PAIRS

  primary metric              guardrail                  what goes wrong alone
  ---------------------------------------------------------------------------
  sessions per user           session quality / depth    autoplay, notification
                                                         spam
  conversion rate             refund + chargeback rate   dark patterns
  tickets closed per agent    reopen rate + CSAT         closing without fixing
  deploy frequency            change-failure rate        shipping breakage
  model accuracy              per-segment accuracy       ignoring a minority
                                                         group
  ad revenue per session      churn rate                 ad load that drives
                                                         people away

  Every primary metric ships with the guardrail that catches its cheat.
""")

# --- A METRICS LAYER: define once, derive everywhere. ---
print("""
  WITHOUT A METRICS LAYER          WITH ONE

  dashboard A: its own SQL          metrics.yml / semantic layer
  dashboard B: its own SQL   ---->    weekly_active_users:
  notebook C:  its own SQL              filter, exclusions, tz, grain
  the board deck: a spreadsheet             |
                                    +-------+--------+--------+
  four numbers, all called          dashboard   notebook   deck
  "active users"                    one definition, four consumers

  Implementations: dbt metrics, Cube, LookML, MetricFlow, or simply one
  SQL view per metric that everything else selects from. The tool matters
  far less than the rule: a metric is defined exactly once.
""")

  FIVE PARTS OF A DEFINITION

    1. EVENT/ENTITY   "distinct user_id with a session >= 10s"
                      not "active user"
    2. TIME BASIS     "rolling 7d, event timestamp, UTC"
                      not "weekly"
    3. POPULATION     "excluding internal, bots, test accounts"
                      the most common source of divergence
    4. AGGREGATION    "count_distinct", and for ratios:
                      ratio-of-sums, NOT mean-of-ratios
    5. GRAIN          "one row per (week_start, user_id)"

  ONE PHRASE, SIX ANSWERS

    "active users in June"
      1,043   any event, UTC calendar month
        998   excluding internal accounts
        933   also excluding bots
        871   also requiring a session >= 10s
        864   the same, but on the user's LOCAL month
        879   rolling 30 days to 30 June
    A 17% spread, and every one is defensible.

  RATIO OF SUMS VS MEAN OF RATIOS

    country  sessions  conv   rate
    US        80,000   3,200  4.00%
    CA         9,000     450  5.00%
    DE         6,000     240  4.00%
    JP           500      45  9.00%

    mean of rates          5.50%   <- JP's 500 sessions count as much as
                                      the US's 80,000
    total conv / total     4.02%   <- "the conversion rate"

  LEADING VS LAGGING, INPUT VS OUTPUT

                  INPUT (a team can move it)   OUTPUT (the goal)
    LEADING       page load p95                trial-to-paid rate
                  activation within 7d
    LAGGING       -                            monthly revenue
                                               net revenue retention

    Measure teams on outputs; give them inputs to work with.

Key Takeaways

• A metric definition needs the precise event, the time basis with a time zone, the exclusions, the aggregation, and the grain, or two implementations will diverge.
• One vague phrase such as "active users" yields numbers that differ by tens of percent, every one of them defensible.
• For ratios, specify ratio-of-sums rather than mean-of-ratios, since the difference drives a large share of reconciliation disputes.
• Lagging output metrics are a history lesson; pair them with leading input metrics that a team can actually move.
• Give every primary metric a guardrail, because a measure used as a target invites moving the number without producing the value.

🧪 Practice

1. Compute "monthly active users" six ways by varying exclusions, the time zone, and calendar versus rolling windows, and report the spread.
2. For a segmented conversion table, compute the mean of rates and the ratio of sums, and construct segment sizes that make them differ by two points.
3. Write a metric definition dataclass with a validator that rejects undefined adjectives, a missing time zone, and a ratio without a rule.
4. Pick three metrics your team reports and name the guardrail each needs.
5. Interview: Finance and product report different monthly revenue and both insist they are right. How do you resolve it? (Hint: which of the five definition parts is most likely to differ, and what artefact ends the argument permanently?)

Self-Service BI Tools

Business intelligence tools exist because the bottleneck in most organisations is not analysis, it is the analyst. Every routine question that requires a person consumes hours and adds delay, and most of those questions are variations on ones already answered. Self-service tooling aims to let people answer their own variations.

The category divides along one axis that matters more than any feature list: who defines the metrics?

In a semantic-layer tool, an analytics engineer defines dimensions, measures, and joins centrally, and business users explore within that model. Looker's LookML is the reference example. The benefit is that everyone's numbers agree, because there is one definition. The cost is that every new dimension requires a modelling change, which makes the tool feel slow to users who want something now.

In a direct-query tool, users connect to tables and build their own logic. Tableau, Power BI, and most notebook-based approaches work this way. The benefit is speed and flexibility. The cost is that definitions multiply, and within a year several dashboards report different values for the same metric name.

Neither is correct in general. The practical resolution most organisations reach is a layered model: certified, governed metrics for the numbers that matter, such as anything reported to leadership or used in compensation, and free exploration for everything else, clearly marked as uncertified.

The landscape in brief:

Tool	Model	Strength	Watch out for
Looker	semantic layer, LookML	governed consistency	modelling effort, cost
Tableau	direct query, visual	visualisation depth	definition sprawl
Power BI	direct query, DAX	Microsoft integration, price	DAX learning curve
Metabase	light semantic layer	easy setup, open source	limited modelling
Superset	SQL-first, open source	free, extensible	self-hosted operations
Mode, Hex	notebook plus BI	analyst-friendly, SQL first	less end-user-friendly
dbt plus a semantic layer	definitions in the warehouse	version-controlled metrics	still needs a front end

Whatever the tool, four practices determine whether self-service works.

Model the data for the audience. Business users cannot navigate a normalised schema with forty tables. Build wide, denormalised, well-named tables where the joins are already done and the columns mean what they say. Most self-service failures are data-modelling failures wearing a tool's name.

Certify a small set of dashboards and mark the rest. A visible distinction between "this is the official number" and "someone built this for themselves" prevents the most damaging outcome, which is an ad-hoc chart being presented to an executive as fact.

Manage the sprawl actively. Dashboards accumulate. Track usage, archive what nobody opens, and require an owner for anything certified. An unmaintained dashboard is worse than no dashboard because it looks authoritative.

Control query cost. Self-service against a cloud warehouse means users can run expensive queries by accident. Aggregate tables, result caching, partition limits, and per-user budgets are not optional at scale.

The honest trade-off to state plainly: self-service reduces analyst load on routine questions and increases governance work. The total effort does not vanish, it moves from answering questions to maintaining a model. That is usually a good trade, and it is a trade rather than a saving.

import numpy as np
import pandas as pd

rng = np.random.default_rng(20)

# --- WHY MODELLING MATTERS MORE THAN THE TOOL ---
# A normalised schema: correct, efficient, and unusable by a business user.
orders = pd.DataFrame({
    "order_id": range(1, 9),
    "cust_id": [101, 102, 101, 103, 104, 102, 105, 103],
    "prod_id": [7, 8, 9, 7, 8, 7, 9, 8],
    "chan_id": [1, 2, 1, 3, 2, 1, 3, 2],
    "amt_cents": [4999, 12900, 899, 4999, 12900, 4999, 899, 12900],
    "st": ["C", "C", "R", "C", "P", "C", "C", "R"],   # what is "R"? "P"?
})
customers = pd.DataFrame({"cust_id": [101, 102, 103, 104, 105],
                          "seg_cd": ["E", "S", "M", "S", "E"],
                          "ctry": ["US", "CA", "US", "DE", "US"]})
products = pd.DataFrame({"prod_id": [7, 8, 9],
                         "nm": ["Basic", "Pro", "Addon"],
                         "cat": ["plan", "plan", "addon"]})
channels = pd.DataFrame({"chan_id": [1, 2, 3],
                         "chan": ["direct", "partner", "organic"]})

print("NORMALISED: a user must know 4 tables, 3 joins, and 5 code mappings")
print(orders.head(3).to_string(index=False))

# The same data, modelled FOR the audience: one wide table, decoded, named.
STATUS = {"C": "completed", "R": "refunded", "P": "pending"}
SEGMENT = {"E": "enterprise", "S": "smb", "M": "mid-market"}

wide = (orders
        .merge(customers, on="cust_id")
        .merge(products, on="prod_id")
        .merge(channels, on="chan_id")
        .assign(
            order_status=lambda d: d["st"].map(STATUS),
            customer_segment=lambda d: d["seg_cd"].map(SEGMENT),
            revenue_usd=lambda d: (d["amt_cents"] / 100).round(2),
            # Pre-compute the flag every definition argument is about.
            is_net_revenue=lambda d: d["st"] == "C",
        )
        .rename(columns={"ctry": "customer_country", "nm": "product_name",
                         "cat": "product_category", "chan": "acquisition_channel"})
        [["order_id", "order_status", "customer_segment", "customer_country",
          "product_name", "product_category", "acquisition_channel",
          "revenue_usd", "is_net_revenue"]])
print("\nMODELLED FOR THE AUDIENCE: one table, decoded, self-describing")
print(wide.head(4).to_string(index=False))

# A business user can now answer their own question with no joins.
print("\na user's own question, answered without an analyst:")
print(wide[wide["is_net_revenue"]]
      .groupby(["customer_segment", "acquisition_channel"])["revenue_usd"]
      .sum().round(2).to_string())

# --- THE SEMANTIC LAYER: define once, in version control. ---
SEMANTIC_MODEL = {
    "model": "orders",
    "source": "analytics.fct_orders_wide",
    "dimensions": {
        "customer_segment": {"type": "string", "label": "Customer segment"},
        "acquisition_channel": {"type": "string", "label": "Channel"},
        "order_month": {"type": "time", "grain": "month", "tz": "UTC"},
    },
    "measures": {
        "net_revenue": {
            "sql": "SUM(CASE WHEN is_net_revenue THEN revenue_usd ELSE 0 END)",
            "label": "Net revenue (USD)",
            "description": "Completed orders only; refunds and pending excluded",
            "owner": "finance-analytics", "certified": True,
        },
        "gross_revenue": {
            "sql": "SUM(revenue_usd)",
            "label": "Gross revenue (USD)",
            "description": "All orders regardless of status",
            "owner": "finance-analytics", "certified": True,
        },
        "conversion_rate": {
            "sql": "SUM(conversions) / NULLIF(SUM(sessions), 0)",
            "label": "Conversion rate",
            # State the rule explicitly: this is a ratio of sums.
            "description": "Ratio of sums, never the mean of per-group rates",
            "owner": "growth-analytics", "certified": True,
        },
    },
}
print(f"\nsemantic model exposes {len(SEMANTIC_MODEL['measures'])} certified"
      f" measures over {len(SEMANTIC_MODEL['dimensions'])} dimensions")
print("Both revenue definitions exist, are NAMED differently, and are")
print("documented. The dispute is resolved by naming, not by argument.")

# --- CERTIFICATION: the visible distinction that prevents the worst outcome. ---
catalogue = pd.DataFrame([
    ("Exec weekly review",   "certified",   "finance-analytics", 412, 2),
    ("Growth funnel",        "certified",   "growth-analytics",  260, 5),
    ("Support SLA",          "certified",   "support-ops",       118, 9),
    ("my churn exploration", "uncertified", "a.patel",             3, 240),
    ("Q2 campaign v3 FINAL", "uncertified", "unowned",             1, 310),
    ("Copy of Copy of KPIs", "uncertified", "unowned",             0, 380),
], columns=["dashboard", "status", "owner", "views_30d", "days_since_edit"])

catalogue["action"] = np.select(
    [(catalogue["status"] == "certified") & (catalogue["days_since_edit"] > 90),
     (catalogue["views_30d"] == 0) & (catalogue["days_since_edit"] > 180),
     (catalogue["views_30d"] < 5) & (catalogue["owner"] == "unowned"),
     catalogue["status"] == "certified"],
    ["REVIEW: certified but stale", "ARCHIVE: unused and abandoned",
     "ARCHIVE: orphaned", "keep"],
    default="keep as personal")
print("\nDASHBOARD LIFECYCLE AUDIT")
print(catalogue.to_string(index=False))

# --- QUERY COST: self-service means users can spend money by accident. ---
queries = pd.DataFrame({
    "dashboard": ["Exec weekly", "Growth funnel", "Ad-hoc explore",
                  "Support SLA", "Ad-hoc explore"],
    "gb_scanned": [0.4, 2.1, 840.0, 0.9, 1_250.0],
    "runs_per_day": [30, 120, 8, 200, 4],
})
COST_PER_TB = 6.25
queries["daily_usd"] = (queries["gb_scanned"] / 1024 * COST_PER_TB
                        * queries["runs_per_day"]).round(2)
queries["monthly_usd"] = (queries["daily_usd"] * 30).round(0)
print("\nQUERY COST BY DASHBOARD")
print(queries.sort_values("monthly_usd", ascending=False).to_string(index=False))
total = queries["monthly_usd"].sum()
adhoc = queries[queries["dashboard"] == "Ad-hoc explore"]["monthly_usd"].sum()
print(f"\ntotal ${total:,.0f}/mo, of which ad-hoc exploration is "
      f"${adhoc:,.0f} ({adhoc/total:.0%})")
print("Two unpartitioned exploratory queries cost more than every governed")
print("dashboard combined. Aggregate tables and partition limits, not policy.")

# --- THE TRADE-OFF, STATED HONESTLY ---
print("""
  WHAT SELF-SERVICE ACTUALLY TRADES

  decreases                          increases
  ---------------------------------  ---------------------------------------
  analyst time on routine questions  data-modelling and maintenance work
  turnaround time for the requester  governance: certification, definitions
  the analyst as a bottleneck        query cost exposure
                                     risk of conflicting numbers in circulation

  The effort does not vanish; it moves from ANSWERING to MODELLING. That is
  usually a good trade, and it is a trade, not a saving.

  THE LAYERED RESOLUTION

    certified   leadership reporting, compensation, board numbers
                -> semantic layer, named owner, change review
    exploratory everything else
                -> direct query, clearly marked uncertified, no SLA
""")

  THE AXIS THAT MATTERS: WHO DEFINES THE METRICS?

  SEMANTIC LAYER                     DIRECT QUERY
  analytics engineer defines         users define their own
  dimensions + measures centrally    logic per dashboard
        |                                  |
        v                                  v
  everyone's numbers AGREE           fast and flexible
  new dimension = a modelling task   definitions MULTIPLY
  feels slow to users                4 dashboards, 4 "revenues"

  Resolution: layer them. Certified metrics for what matters, free
  exploration for everything else, and mark which is which.

  MODEL FOR THE AUDIENCE

  NORMALISED (correct, unusable)     WIDE (usable)
  orders    + customers              fct_orders_wide
  + products + channels              order_status   "completed"
  st = "C", seg_cd = "E"             customer_segment "enterprise"
  4 tables, 3 joins, 5 code maps     0 joins, 0 lookups
                                     is_net_revenue  true

  Most "the BI tool is bad" complaints are data-modelling complaints.

  THE FOUR PRACTICES

    1. MODEL for the audience     wide, denormalised, decoded, well named
    2. CERTIFY a small set        and visibly mark everything else
    3. MANAGE sprawl              track usage, archive the unused, require
                                  an owner
    4. CONTROL cost               aggregates, caching, partition limits,
                                  budgets

Key Takeaways

• The defining question for a business intelligence tool is who defines the metrics, centrally in a semantic layer or per dashboard by users.
• Most self-service failures are data-modelling failures: business users need wide, denormalised, decoded tables with the joins already done.
• Certify a small governed set and visibly mark everything else, so an ad-hoc chart is never presented as the official number.
• Dashboards accumulate, so track usage, require owners, and archive actively, because an unmaintained dashboard still looks authoritative.
• Self-service moves effort from answering questions to maintaining a model and controlling query cost, which is a trade rather than a saving.

🧪 Practice

1. Take a normalised four-table schema and produce one wide, decoded table a non-technical user could query without joins.
2. Write a semantic model dictionary defining gross and net revenue as separately named certified measures with descriptions and owners.
3. Build a dashboard lifecycle audit that flags certified-but-stale, unused, and orphaned dashboards from usage and edit dates.
4. Estimate the monthly warehouse cost of five queries from their bytes scanned and run frequency, and identify which one to fix first.
5. Interview: Leadership asks you to "let the business self-serve so analysts can focus on deeper work". What do you tell them it will require? (Hint: what work appears that did not exist before, and what happens to the numbers without it?)

Automated Reporting Pipelines

Any report you produce more than twice should be automated. The reason is not only the time saved, it is that a manual report is unreproducible and inconsistent: someone forgets a filter, runs it a day early, or copies last week's chart by accident. Automation makes the report a deterministic function of the data plus a date.

The pieces of an automated report pipeline:

A parameterised source. One script or notebook taking a date range and any filters as parameters, with no hard-coded values and no manual steps. If it has a cell you have to edit, it is not automated.

Scheduled execution with retries and a timeout. Cron is the simplest form; an orchestrator such as Airflow, Dagster, or Prefect adds dependency management, backfills, and alerting, which you need as soon as a report depends on an upstream job.

Dependency awareness is the piece people skip and regret. A report that runs at a fixed time regardless of whether its source tables loaded will eventually publish a report built on yesterday's data or on a half-loaded table. The report should depend on the data being ready, not on the clock.

Validation before publication. This is the most important element and the one most often missing. Check row counts against expectations, check that key metrics fall in plausible ranges, check for nulls in required fields, and check that the period is complete. A pipeline that publishes a wrong report on schedule is worse than one that fails, because a failure gets noticed and a wrong number gets believed and acted on.

Rendering and delivery. Produce the output format your audience actually uses, and deliver it where they already are. A report emailed to people who work in Slack does not get read. Include a link back to the source dashboard so anyone can go deeper.

Observability of the pipeline itself. Log every run, alert on failures to a person rather than a mailbox nobody watches, and make the last successful run time visible in the report. A silently broken report is the standard failure mode, and the usual way it manifests is that nobody notices for a month.

Two design principles worth adopting.

Idempotency. Running the report for the same date twice produces the same output, and re-running never duplicates or corrupts anything. This makes retries and backfills safe, which is what makes the pipeline operable.

Alert on exceptions, not on schedule. A report that arrives every morning whether or not anything changed trains its recipients to ignore it. The higher-value pattern for operational metrics is to send something only when a threshold is crossed or an anomaly appears, with the routine numbers available on a dashboard for anyone who wants them.

import numpy as np
import pandas as pd
from dataclasses import dataclass
from datetime import datetime, timezone

rng = np.random.default_rng(30)

# --- 1. PARAMETERISED, not hard-coded. Everything is an argument. ---
@dataclass(frozen=True)
class ReportParams:
    as_of: pd.Timestamp                  # the run date, never datetime.now()
    lookback_days: int = 28
    region: str | None = None
    min_rows_expected: int = 1000

    @property
    def start(self) -> pd.Timestamp:
        return self.as_of - pd.Timedelta(days=self.lookback_days)

def extract(params: ReportParams) -> pd.DataFrame:
    """Stand-in for a warehouse query. Note: bounded by params, not by now()."""
    days = pd.date_range(params.start, params.as_of - pd.Timedelta("1d"), freq="D")
    rows = []
    for d in days:
        for region in ["north", "south", "east"]:
            if params.region and region != params.region:
                continue
            base = {"north": 1400, "south": 900, "east": 480}[region]
            rows.append({"date": d, "region": region,
                         "sessions": int(base + rng.normal(0, base * .07)),
                         "orders": int(base * .042 + rng.normal(0, 6))})
    return pd.DataFrame(rows)

# --- 2. VALIDATION GATES: fail loudly rather than publish quietly wrong. ---
@dataclass
class Check:
    name: str
    passed: bool
    detail: str
    severity: str = "error"          # "error" blocks publication, "warn" notes it

def validate(df: pd.DataFrame, params: ReportParams) -> list[Check]:
    checks = []

    checks.append(Check("row_volume", len(df) >= params.min_rows_expected // 10,
                        f"{len(df)} rows"))

    # Completeness: every expected date present, and the period is finished.
    expected_days = params.lookback_days
    actual_days = df["date"].nunique()
    checks.append(Check("date_completeness", actual_days == expected_days,
                        f"{actual_days} of {expected_days} days present"))

    latest = df["date"].max()
    checks.append(Check("no_partial_period",
                        latest < params.as_of.normalize(),
                        f"latest date {latest:%Y-%m-%d}, as_of "
                        f"{params.as_of:%Y-%m-%d}"))

    # Required fields must not be null.
    nulls = df[["sessions", "orders"]].isna().sum().sum()
    checks.append(Check("no_nulls", nulls == 0, f"{nulls} nulls in key fields"))

    # Plausibility: a metric outside a sane band means something broke upstream.
    rate = df["orders"].sum() / max(df["sessions"].sum(), 1)
    checks.append(Check("rate_plausible", 0.005 <= rate <= 0.20,
                        f"conversion rate {rate:.3%}"))

    # Continuity: a sudden level shift is usually a pipeline change, not news.
    daily = df.groupby("date")["sessions"].sum()
    if len(daily) > 7:
        recent, prior = daily.tail(7).mean(), daily.head(len(daily) - 7).mean()
        shift = abs(recent / prior - 1)
        checks.append(Check("no_level_shift", shift < 0.35,
                            f"recent 7d differ from prior by {shift:.1%}",
                            severity="warn"))
    return checks

def gate(checks: list[Check]) -> bool:
    """Publication is blocked by any failing error-severity check."""
    for c in checks:
        mark = "PASS" if c.passed else ("FAIL" if c.severity == "error" else "WARN")
        print(f"  [{mark:4s}] {c.name:20s} {c.detail}")
    return all(c.passed for c in checks if c.severity == "error")

params = ReportParams(as_of=pd.Timestamp("2024-07-01"), lookback_days=28)
data = extract(params)
print("VALIDATION")
ok = gate(validate(data, params))
print(f"\npublish? {ok}")

# A broken extract must be BLOCKED, not published.
print("\nVALIDATION against a truncated, partly-null extract:")
broken = data.head(30).copy()
broken.loc[broken.index[:4], "orders"] = np.nan
blocked = gate(validate(broken, params))
print(f"\npublish? {blocked}   <- the gate did its job")

# --- 3. IDEMPOTENCY: same inputs, same output, safe to re-run. ---
def transform(df: pd.DataFrame) -> pd.DataFrame:
    out = (df.groupby("region")
             .agg(sessions=("sessions", "sum"), orders=("orders", "sum"))
             .assign(conv_rate=lambda d: (d["orders"] / d["sessions"]).round(5))
             .sort_values("orders", ascending=False))
    return out

a, b = transform(data), transform(data)
print(f"\nidempotent transform: {a.equals(b)}")

import hashlib
def content_hash(df: pd.DataFrame) -> str:
    return hashlib.sha256(
        pd.util.hash_pandas_object(df, index=True).values.tobytes()).hexdigest()[:12]
print(f"content hash: {content_hash(a)}  (same params -> same hash -> skip resend)")

# --- 4. DEPENDENCY ON DATA READINESS, not on the clock. ---
UPSTREAM = {
    "analytics.fct_sessions": pd.Timestamp("2024-07-01 04:12", tz="UTC"),
    "analytics.fct_orders":   pd.Timestamp("2024-07-01 04:31", tz="UTC"),
    "analytics.dim_region":   pd.Timestamp("2024-06-29 04:05", tz="UTC"),  # stale
}
def upstream_ready(freshness: dict, as_of: pd.Timestamp,
                   max_age_h: float = 26) -> list[str]:
    stale = []
    for table, loaded in freshness.items():
        age_h = (as_of - loaded).total_seconds() / 3600
        if age_h > max_age_h:
            stale.append(f"{table} is {age_h:.0f}h old")
    return stale

run_at = pd.Timestamp("2024-07-01 06:00", tz="UTC")
stale = upstream_ready(UPSTREAM, run_at)
print(f"\nupstream freshness check: {stale or 'all fresh'}")
print("A clock-triggered report would have published on a 2-day-old dimension")
print("table. Depend on data readiness, not on 06:00.")

# --- 5. EXCEPTION-BASED ALERTING beats a daily digest nobody reads. ---
daily_totals = data.groupby("date")["sessions"].sum()
THRESHOLDS = {"sessions_min": 2200, "z_score_max": 3.0}

def find_exceptions(series: pd.Series, thresholds: dict) -> list[str]:
    alerts = []
    latest_val = series.iloc[-1]
    if latest_val < thresholds["sessions_min"]:
        alerts.append(f"sessions {latest_val:,.0f} below floor "
                      f"{thresholds['sessions_min']:,}")
    # Robust z-score: median and MAD are not dragged by the anomaly itself.
    med = series.iloc[:-1].median()
    mad = (series.iloc[:-1] - med).abs().median()
    if mad > 0:
        z = 0.6745 * (latest_val - med) / mad
        if abs(z) > thresholds["z_score_max"]:
            alerts.append(f"sessions {latest_val:,.0f} is {z:+.1f} robust SDs "
                          f"from the {med:,.0f} median")
    return alerts

print(f"\nexceptions today: {find_exceptions(daily_totals, THRESHOLDS) or 'none'}")
shocked = daily_totals.copy()
shocked.iloc[-1] = shocked.iloc[-1] * 0.45          # simulate an incident
print(f"after a 55% drop:  {find_exceptions(shocked, THRESHOLDS)}")

# --- 6. RENDER with provenance baked in. ---
def render_summary(table: pd.DataFrame, params: ReportParams,
                   checks: list[Check], run_at: pd.Timestamp) -> str:
    warns = [c.name for c in checks if not c.passed and c.severity == "warn"]
    lines = [
        f"Weekly performance summary -- {params.start:%d %b} to "
        f"{params.as_of - pd.Timedelta('1d'):%d %b %Y}",
        "",
        f"Overall conversion: {table['orders'].sum()/table['sessions'].sum():.2%}"
        f" on {table['sessions'].sum():,} sessions",
        "",
        table.assign(conv_rate=lambda d: (d["conv_rate"]*100).round(2))
             .to_string(),
        "",
        f"Window: {params.lookback_days} days ending "
        f"{params.as_of - pd.Timedelta('1d'):%Y-%m-%d} (UTC)",
        f"Region filter: {params.region or 'all'}",
        f"Generated: {run_at:%Y-%m-%d %H:%M} UTC | "
        f"content hash: {content_hash(table)}",
        f"Source: analytics.fct_sessions, analytics.fct_orders",
        f"Checks: {sum(c.passed for c in checks)}/{len(checks)} passed"
        + (f" | warnings: {', '.join(warns)}" if warns else ""),
        "Full dashboard: https://bi.internal/d/weekly-performance",
    ]
    return "\n".join(lines)

if ok:
    print("\n" + "=" * 62)
    print(render_summary(a, params, validate(data, params), run_at))
    print("=" * 62)

# --- 7. RUN LOG: how you find out a report has been silently broken. ---
run_log = pd.DataFrame([
    ("2024-06-28", "success", 84, 6, 0, 12.4),
    ("2024-06-29", "success", 84, 6, 0, 11.8),
    ("2024-06-30", "failed",   0, 0, 1, 2.1),
    ("2024-07-01", "success", 84, 6, 1, 13.0),
], columns=["run_date", "status", "rows", "checks_passed", "warnings", "seconds"])
print("\nRUN LOG (surface this, or a broken report goes unnoticed for weeks)")
print(run_log.to_string(index=False))
streak = (run_log["status"] == "success").iloc[::-1].cumprod().sum()
print(f"\nconsecutive successful runs: {streak}")
print("Alert a PERSON on failure. A mailbox nobody reads is not monitoring.")

  THE PIPELINE

    parameters (as_of, filters)
         |
    upstream freshness check  ----fail----> ALERT, do not run
         |
    extract (bounded by params, never by now())
         |
    validate: volume, completeness, nulls, plausibility, continuity
         |         |
         |         +----fail----> ALERT + BLOCK PUBLICATION
         v
    transform (idempotent: same params -> same content hash)
         |
    render (with provenance: window, filters, source, run time, hash)
         |
    deliver where the audience already works + link to the dashboard
         |
    log the run; alert a PERSON on failure

  WHY VALIDATION IS THE MOST IMPORTANT STEP

    a pipeline that FAILS        gets noticed and fixed
    a pipeline that publishes    gets BELIEVED and acted on
      a wrong number on time

    Failing loudly is a feature.

  SCHEDULE-TRIGGERED VS DEPENDENCY-TRIGGERED

    cron at 06:00               "when fct_orders for date D has loaded"
    runs whether or not the      runs when the data is actually ready
    data arrived
    -> eventually publishes      -> late rather than wrong
       yesterday's numbers

  DIGEST VS EXCEPTION ALERTING

    daily digest                 exception alert
    arrives every morning        arrives when a threshold is crossed
    read for two weeks, then     read every time, because it MEANS something
    filtered to a folder
    -> routine numbers belong on a dashboard; alerts belong in the inbox

Key Takeaways

• Anything produced more than twice should be parameterised and scheduled, so the report becomes a deterministic function of the data and a date.
• Validate before publishing, because a pipeline that publishes a wrong number on time is worse than one that fails visibly.
• Trigger on upstream data readiness rather than on the clock, or the report eventually renders stale or half-loaded tables.
• Make runs idempotent so retries and backfills are safe, and bound every query by a parameter rather than by the current time.
• Alert on exceptions rather than sending a digest nobody reads, and surface the last successful run so silent breakage becomes visible.

🧪 Practice

1. Write a parameterised report function whose window comes from an as_of argument rather than the current time, and run it for two historical dates.
2. Implement validation checks for row volume, date completeness, nulls, and metric plausibility, then show them blocking a deliberately broken extract.
3. Add a content hash to a report so a re-run with identical inputs can be detected and skipped.
4. Write a robust anomaly check using the median and median absolute deviation and confirm it fires on an injected 50% drop.
5. Interview: A weekly report has been arriving on time for two months and you discover it has been using stale data for three weeks. What went wrong and what do you add? (Hint: what did the pipeline check before publishing, and what should it have depended on?)

Narrative Structure in Reports

A finding that nobody understands has the same value as no finding. The gap between good analysis and adopted analysis is almost always communication, and the specific failure is nearly always the same: the analyst presents the work in the order they did it, rather than in the order the reader needs it.

The instinct to build up from method to result is natural and wrong for this audience. You spent two weeks on the analysis, so the methodology feels like the substance. To a reader deciding whether to act, the methodology is the support for a conclusion they have not heard yet, and they will lose patience before you reach it.

The fix is the inverted pyramid, borrowed from journalism. Lead with the conclusion and the recommended action. Follow with the key supporting evidence. Then the detail, method, and caveats. Then the appendix. A reader who stops after the first paragraph has the decision-relevant content, and a reader who wants the method can keep going.

For a recommendation, a useful structure is situation, complication, question, answer: the context everyone agrees on, the thing that changed, the question it raises, and your answer. This works because it earns agreement before it asks for any, which is what makes an unwelcome conclusion easier to accept.

Several practices matter regardless of structure.

State the action, not just the finding. "Conversion fell 1.8 points" is a finding. "Conversion fell 1.8 points after the v2.4 release; roll back the address-form change this week" is a recommendation. Analysts systematically under-recommend, usually from a reasonable caution about overreach, and the result is that the decision defaults to whoever speaks with most confidence.

Quantify the consequence in the reader's units. Not "a 1.8 point drop" but "about 340 lost orders a week, roughly 52,000 dollars a month". The reader is making a resource decision and needs the magnitude in the terms that decision uses.

State uncertainty honestly and without hedging everything. Give the interval and the main threats to the conclusion, then still make a recommendation. A report that hedges every sentence is unusable, and one that projects false certainty is worse. The workable form is: here is the estimate, here is its range, here is what would change my mind, here is what I recommend given that.

Adapt depth to the audience. An executive needs the conclusion, the size, and the decision. A peer analyst needs the method and the data well enough to challenge it. A product team needs the mechanism and what to change. The same analysis becomes three different documents, and the most common error is sending the peer-analyst version to everyone.

Write for the skimmer. Headings that state findings rather than topics. Charts whose titles carry the conclusion. Bold for the load-bearing sentence, not for whole paragraphs. Most readers will read the headings and look at the charts.

Finally, say what would change your mind. Naming the evidence that would overturn your conclusion is the strongest available signal that you have thought about it properly, and it converts a defensive presentation into a collaborative one.

import numpy as np
import pandas as pd

rng = np.random.default_rng(40)

# --- THE ANALYSIS BEHIND THE REPORT ---
days = pd.date_range("2024-02-26", periods=56, freq="D")
RELEASE = pd.Timestamp("2024-03-18")
pre_rate, post_rate = 0.0412, 0.0234
sessions = (11_000 + rng.normal(0, 700, 56)).round(0)
rate = np.where(days < RELEASE, pre_rate, post_rate) + rng.normal(0, .0012, 56)
funnel = pd.DataFrame({"date": days, "sessions": sessions,
                       "orders": (sessions * rate).round(0)})
funnel["conv"] = funnel["orders"] / funnel["sessions"]
pre = funnel[funnel["date"] < RELEASE]
post = funnel[funnel["date"] >= RELEASE]

drop_pts = (post["conv"].mean() - pre["conv"].mean()) * 100
weekly_sessions = post["sessions"].mean() * 7
lost_orders_wk = weekly_sessions * abs(drop_pts) / 100
AOV = 38.50
monthly_usd = lost_orders_wk * 52 / 12 * AOV

# A confidence interval on the difference of two proportions.
p1, n1 = pre["orders"].sum()/pre["sessions"].sum(), pre["sessions"].sum()
p2, n2 = post["orders"].sum()/post["sessions"].sum(), post["sessions"].sum()
se = np.sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)
ci = ((p2 - p1 - 1.96*se) * 100, (p2 - p1 + 1.96*se) * 100)

print(f"pre-release  {p1:.2%} on {n1:,.0f} sessions")
print(f"post-release {p2:.2%} on {n2:,.0f} sessions")
print(f"difference   {drop_pts:+.2f} pts, 95% CI [{ci[0]:+.2f}, {ci[1]:+.2f}]")
print(f"consequence  ~{lost_orders_wk:,.0f} orders/week,"
      f" ~${monthly_usd:,.0f}/month at ${AOV} AOV")

# --- THE SAME ANALYSIS, TWO ORDERINGS ---
print("""
=========================== WRONG ORDER ============================
1. Data and methodology
   We queried fct_checkout_sessions for 26 Feb to 21 Apr, excluding
   bots via the is_bot flag and internal traffic via account_state...
2. Approach
   We considered a difference-in-means test and a regression
   discontinuity, and settled on...
3. Robustness checks
   We re-ran excluding weekends, then excluding mobile...
4. Results
   Conversion declined by 1.78 percentage points...
5. Conclusion
   The v2.4 release appears to have reduced conversion.

The reader stopped at section 2. The decision was made by whoever
spoke most confidently in the meeting.
====================================================================""")

def executive_report() -> str:
    return f"""
=========================== RIGHT ORDER ============================
CONVERSION FELL {abs(drop_pts):.1f} POINTS AFTER v2.4 -- RECOMMEND ROLLBACK
THIS WEEK

Checkout conversion dropped from {p1:.2%} to {p2:.2%} the day the v2.4
release shipped and has not recovered in {len(post)} days. That is about
{lost_orders_wk:,.0f} lost orders a week, roughly ${monthly_usd:,.0f} a month.

RECOMMENDATION
  Roll back the address-form change in v2.4 this week. If the rollback
  restores conversion, re-ship the change behind an experiment.

EVIDENCE
  - The drop is same-day with the release, not gradual.
  - Size: {drop_pts:+.2f} points, 95% CI [{ci[0]:+.2f}, {ci[1]:+.2f}] --
    the interval excludes zero by a wide margin.
  - Traffic volume is unchanged ({pre['sessions'].mean():,.0f} vs
    {post['sessions'].mean():,.0f} sessions/day), so this is not a mix shift.
  - The decline is present on every device and channel, so it is not
    one segment.

WHAT WOULD CHANGE MY MIND
  - A tracking change shipped in the same release, which would make this
    a measurement artefact rather than a real decline. Checked: the
    checkout events are unchanged, but I have not audited the tag config.
  - A concurrent pricing or promotion change on the same date.

CAVEATS
  This is a before-and-after comparison, not an experiment, so it cannot
  fully separate the release from anything else that changed that day.

DETAIL AND METHOD: appendix, page 2.
====================================================================
"""
print(executive_report())

# --- SITUATION / COMPLICATION / QUESTION / ANSWER ---
print("""
  SCQA, FOR AN UNWELCOME CONCLUSION

  SITUATION     Checkout conversion held steady near 4.1% through Q1.
                (everyone already agrees)
  COMPLICATION  It fell to 2.3% on 18 March and has not recovered.
                (the thing that changed)
  QUESTION      What caused it and what do we do?
                (the question the reader now has)
  ANSWER        The v2.4 address-form change. Roll it back this week.
                (your recommendation)

  It earns agreement BEFORE asking for any, which is what makes a
  conclusion nobody wants to hear easier to accept.
""")

# --- ONE ANALYSIS, THREE AUDIENCES ---
audiences = pd.DataFrame([
    ("Executive", "conclusion + size + decision", "1 slide / 5 lines",
     "money, timeline, risk", "method, code, robustness checks"),
    ("Product team", "mechanism + what to change", "1 page + 2 charts",
     "which change, which users, what to ship", "statistical detail"),
    ("Analyst peer", "method + data + threats", "full doc + notebook",
     "identification, robustness, alternatives", "nothing"),
], columns=["audience", "needs", "length", "emphasise", "omit"])
print(pd.DataFrame(audiences).to_string(index=False))
print("\nThe usual error is sending the analyst-peer version to all three.")

# --- FINDING VS RECOMMENDATION ---
print("""
  STATE THE ACTION

  finding only                        recommendation
  ----------------------------------  ------------------------------------
  "conversion fell 1.8 points"        "roll back the v2.4 address form
                                       this week; re-ship behind an
                                       experiment"
  "segment C churns 2x"               "move segment C to the onboarding
                                       flow we built for B; expect ~200
                                       retained accounts/quarter"
  "the model reaches 0.82 AUC"        "ship at threshold 0.11, which flags
                                       ~400 cases/day for a team that can
                                       review 450"

  QUANTIFY IN THE READER'S UNITS

    "a 1.8 point drop"                        analyst units
    "~340 orders/week, ~$52k/month"           DECISION units
""")

# --- WRITE FOR THE SKIMMER: headings that carry the findings. ---
print("""
  HEADINGS: TOPICS VS FINDINGS

  topic headings (reader must read on)   finding headings (skimmable)
  ------------------------------------   ------------------------------------
  "Conversion analysis"                  "Conversion fell 1.8 pts on 18 Mar"
  "Segment breakdown"                    "The drop is uniform across segments"
  "Traffic"                              "Traffic volume is unchanged"
  "Methodology"                          "Before/after, not an experiment"
  "Next steps"                           "Roll back this week, then A/B test"

  Reading only the right-hand column gives the entire report. That is the
  test: do your headings alone tell the story?
""")

# --- A REPORT CHECKLIST ---
def audit_report(r: dict) -> list[str]:
    issues = []
    if not r.get("conclusion_first"):
        issues.append("conclusion is not in the first paragraph")
    if not r.get("recommends_action"):
        issues.append("states a finding but no recommended action")
    if not r.get("quantified_in_business_units"):
        issues.append("effect not translated into money, users, or time")
    if not r.get("uncertainty_stated"):
        issues.append("no interval or range on the headline number")
    if not r.get("what_would_change_my_mind"):
        issues.append("no falsifier named: reads as advocacy, not analysis")
    if not r.get("findings_in_headings"):
        issues.append("headings are topics, so skimming yields nothing")
    if not r.get("audience_specific"):
        issues.append("one document for every audience")
    if not r.get("source_and_asof"):
        issues.append("not reproducible: no source tables or as-of date")
    return issues

print("audit of a typical first draft:")
for i in audit_report({"uncertainty_stated": True, "source_and_asof": True}):
    print("  -", i)
print("\naudit of the executive version above:",
      audit_report({k: True for k in
                    ["conclusion_first", "recommends_action",
                     "quantified_in_business_units", "uncertainty_stated",
                     "what_would_change_my_mind", "findings_in_headings",
                     "audience_specific", "source_and_asof"]}) or "no issues")

  THE INVERTED PYRAMID

        +-----------------------------------------+
        |  CONCLUSION + RECOMMENDED ACTION        |  everyone reads this
        +-----------------------------------------+
         |  key evidence, sized in business units |  most read this
         +--------------------------------------+
          |  method, caveats, what would      |    some read this
          |  change my mind                   |
          +--------------------------------+
           |  detail, appendix, code      |        two people read this
           +----------------------------+

  The order you DID the work is the reverse of the order to present it.

  SCQA

    SITUATION      the context everyone already agrees on
    COMPLICATION   what changed
    QUESTION       what the complication raises
    ANSWER         your recommendation

  THE SKIM TEST

    Read only your headings and chart titles. Do they tell the whole
    story? If not, they are topics, and your reader gets nothing.

  UNCERTAINTY WITHOUT PARALYSIS

    too hedged      "results may suggest a possible decline, though
                     further work would be needed"        -> unusable
    false certainty "v2.4 caused a 1.8 point drop"        -> overclaims
    workable        "1.8 point drop, 95% CI [1.5, 2.1]. Before/after,
                     not an experiment. A tracking change on the same
                     date would explain it; I checked the events but
                     not the tag config. Recommend rollback."

  ONE ANALYSIS, THREE DOCUMENTS

    executive     conclusion, size in dollars, the decision       5 lines
    product       mechanism, which users, what to change          1 page
    analyst peer  identification, robustness, alternatives        full doc

Key Takeaways

• Present in the inverted pyramid, leading with the conclusion and the recommended action, because the order you did the work is the reverse of the order to communicate it.
• State an action rather than only a finding; analysts systematically under-recommend and the decision then defaults to the most confident voice.
• Quantify the consequence in the reader's units, such as orders and money, rather than in percentage points.
• Give the interval and the threats to the conclusion and still recommend something, since total hedging is as unusable as false certainty.
• Write headings that state findings, and name what would change your mind.

🧪 Practice

1. Take an analysis you have done and write its conclusion and recommended action as a single opening paragraph.
2. Convert a set of topic headings into finding headings and check that reading only the headings conveys the whole story.
3. Translate a percentage-point effect into orders per week and revenue per month, stating the assumptions used.
4. Write the same finding three ways for an executive, a product team, and an analyst peer, with the target length for each.
5. Write the "what would change my mind" section for a recent conclusion, naming at least two falsifiers.
6. Interview: You must tell a director that the feature their team spent a quarter on reduced retention. How do you structure the conversation? (Hint: which structure earns agreement before it asks for any, and what do you name about your own uncertainty?)

---

<a id="8-experimentation-and-decision-science"></a>