Scatter plots: show relationships and correlation

#Why it matters

A scatter plot is how you answer "does X move with Y?", the question behind pricing, marketing spend, delivery times and almost every model a data scientist builds. Visualisation appears in 94% of the Data Analyst and 65% of mid and senior Data Scientist job descriptions we read, and "you have a million points, how do you plot them?" is a common interview question. The scatter is also where analysts most often imply a cause they cannot prove.

#The concept

#Which variable goes on which axis

Put the variable you think explains (the input, the thing you could change or that comes first) on the x-axis, and the outcome on the y-axis: distance on x, fare on y; ad spend on x, sales on y; income on x, life expectancy on y. Readers instinctively read "as x increases, y does this". If neither variable drives the other (height and weight), choose the order that matches the question you were asked.

#Overplotting: when points hide each other

With more than a few thousand points, a plain scatter turns into a solid blob. Worse, it misleads: a region with 600,000 points and a region with 600 both look fully inked, so rare outliers look as important as the dense core. Fixes, in order of effort:

A log colour scale on a hexbin is usually essential: the densest bin often holds thousands of times more points than the median bin.

#Trend lines and their honest use

A fitted line summarises the relationship, and it also invites the reader to extrapolate and to assume cause. Use one when the relationship is roughly linear (on the axes you show), say what was fitted ("least squares on metered trips only"), and never extend it beyond the data. If two populations follow different rules (metered fares and flat airport fares), one line through both describes neither.

#Bubble size by area, not radius

A bubble chart encodes a third number as size. People judge circles by area. If you scale the radius with the value, a country with four times the population gets a circle 16 times as big. In matplotlib, the s= argument of scatter is already an area (points squared), so pass values proportional to the number, never its square. Keep sizes as a secondary encoding: area ranks low in accuracy, so bubbles show "much bigger", not "37% bigger".

#Colour as a third variable, with few categories

Colour works for a category with 2 to 5 values (region, plan type). Beyond that, colours become indistinguishable and the legend becomes a lookup table. Better patterns: grey for everything and one highlight colour for the group you are talking about, or small multiples by category.

#Log axes for heavy-tailed data

GDP per person runs from a few hundred to over a hundred thousand dollars. On a linear axis, most countries pile up in the first centimetre. A log x-axis gives every doubling the same distance, which usually matches the relationship too: going from $1,000 to $2,000 per person matters as much for health as going from $20,000 to $40,000.Uselog-logwhenbothvariablesareheavy-tailed(citypopulationagainstnumberofstartups).Labelticksinplainnumbers($1k, $10k) and write "log scale".

#Label a few points, not all

Label the points the reader will ask about: your own country, the extremes, the points mentioned in your text. Labelling all 200 creates an unreadable mat of text.

#Correlation is not causation, shown visually

A tight scatter shows that two variables move together. It does not show why. Richer countries live longer, but wealth also brings clean water, vaccines and education, each of which could be the cause. And the correlation coefficient itself depends on the scale: it measures linear association, so the same data can give 0.6 on linear axes and 0.86 after a log transform. Report the transform with the number.

#Worked example: GDP per person vs life expectancy, India highlighted

The data is the World Bank's World Development Indicators (CC BY 4.0). Download three indicators as CSV from data.worldbank.org: NY.GDP.PCAP.CD (GDP per capita, current US$), SP.DYN.LE00.IN (life expectancy at birth) and SP.POP.TOTL (population). Each zip holds the data file, whose name ends in a version number (hence the glob), and a Metadata_Country file that tells countries apart from aggregates like "World" or "South Asia". The figures below are from the edition last updated 13 July 2026, using 2024 values.

import glob
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

def wdi(code, year="2024"):
    path = glob.glob(f"API_{code}_DS2_en_csv_v2_*.csv")[0]
    return pd.read_csv(path, skiprows=4).set_index("Country Code")[year].rename(code)

meta = pd.read_csv(glob.glob("Metadata_Country_API_SP.DYN.LE00.IN_*.csv")[0]).set_index("Country Code")
names = pd.read_csv(glob.glob("API_SP.POP.TOTL_DS2_en_csv_v2_*.csv")[0], skiprows=4).set_index("Country Code")["Country Name"]
df = pd.concat([wdi("NY.GDP.PCAP.CD"), wdi("SP.DYN.LE00.IN"), wdi("SP.POP.TOTL"), names], axis=1)
df.columns = ["gdp_pc", "life_exp", "pop", "name"]
df = df.join(meta["Region"]).dropna(subset=["Region"])        # aggregates such as WLD have no region
print(len(df), "economies;", df[["gdp_pc", "life_exp"]].dropna().shape[0], "with both values for 2024")
df = df.dropna(subset=["gdp_pc", "life_exp", "pop"])

print("correlation, linear GDP:", round(df["gdp_pc"].corr(df["life_exp"]), 3))
print("correlation, log GDP:   ", round(np.log10(df["gdp_pc"]).corr(df["life_exp"]), 3))
show = df.loc[["IND", "BGD", "CHN", "USA", "NGA"], ["name", "gdp_pc", "life_exp"]].round(1)
show["pop_million"] = (df["pop"] / 1e6).round(1)
print(show.to_string())

fig, ax = plt.subplots(figsize=(9, 5.5))
others = df.drop("IND")
# Bubble AREA proportional to population: matplotlib's s= is already area (points squared)
size = lambda p: p / 1e6 * 0.6
ax.scatter(others["gdp_pc"], others["life_exp"], s=size(others["pop"]),
           color="#9aa7b4", alpha=0.6, edgecolor="white", linewidth=0.5)
ind = df.loc["IND"]
ax.scatter(ind["gdp_pc"], ind["life_exp"], s=size(ind["pop"]), color="#d9822b", edgecolor="#333")
ax.annotate(f"India\n${ind['gdp_pc']:,.0f}, {ind['life_exp']:.1f} years", (ind["gdp_pc"], ind["life_exp"]),
            xytext=(-95, 30), textcoords="offset points", fontsize=8, arrowprops={"arrowstyle": "-"})
for code in ["BGD", "CHN", "USA", "NGA"]:
    r = df.loc[code]
    ax.annotate(r["name"], (r["gdp_pc"], r["life_exp"]), xytext=(6, -10), textcoords="offset points", fontsize=7)
ax.set_xscale("log")
ax.set_xticks([500, 1000, 2000, 5000, 10000, 20000, 50000, 100000],
              ["$500", "$1k", "$2k", "$5k", "$10k", "$20k", "$50k", "$100k"])
ax.set_xlabel("GDP per person, current US$, 2024 (log scale)")
ax.set_ylabel("life expectancy at birth, years, 2024")
ax.set_title("Richer countries live longer, with sharply diminishing returns\n"
             "Bubble area = population", loc="left", fontsize=10)
ax.spines[["top", "right"]].set_visible(False)
fig.text(0.01, 0.01, "Source: World Bank, World Development Indicators (CC BY 4.0)", fontsize=7, color="#666")
fig.tight_layout(rect=(0, 0.03, 1, 1))
fig.savefig("wdi_scatter.png", dpi=150)

217 economies; 200 with both values for 2024
correlation, linear GDP: 0.601
correlation, log GDP:    0.861
                       name   gdp_pc  life_exp  pop_million
Country Code                                               
IND                   India   2592.0      72.2       1450.9
BGD              Bangladesh   2593.4      74.9        173.6
CHN                   China  13293.1      78.0       1409.0
USA           United States  86169.7      78.9        340.0
NGA                 Nigeria   1084.2      54.6        232.7

What the output tells you:

• Aggregates removed. The file has 265 rows. Joining the metadata and dropping rows without a region (regions, income groups and other aggregates) leaves 217 economies, of which 200 have all three values for 2024. Forgetting this step plots "World" and "South Asia" as if they were countries, and South Asia largely double-counts India.

• Log changes the correlation. 0.601 on linear GDP, 0.861 on log GDP. The relationship is close to linear in the log of income, which is why the log axis is the right picture.

• India against its peers. India and Bangladesh have almost identical GDP per person ($2,592 and $2,593) but Bangladesh's life expectancy is 2.7 years higher (74.9 against 72.2). The highlight makes India findable among 200 bubbles; the Bangladesh label gives the comparison that makes the reader think.

• Bubbles by area. size = p / 1e6 * 0.6 is linear in population, and matplotlib treats s as area. India's bubble has about 8.4 times the area of Bangladesh's, matching 1,450.9 against 173.6 million.

• Causation. The chart supports "richer countries tend to live longer". It does not support "raising GDP will raise life expectancy by N years". Nigeria (54.6 years at $1,084) and the countries above it at the same income show how much else is going on.

#Worked example: 3 million taxi trips, distance vs fare

The NYC TLC yellow taxi trips for January 2025, as in the earlier lessons (TLC trip record data).

import duckdb
import matplotlib.pyplot as plt
import numpy as np

t = duckdb.sql("""
    select trip_distance as miles, fare_amount as fare, RatecodeID as rate
    from 'yellow_tripdata_2025-01.parquet'
    where trip_distance > 0 and trip_distance < 60 and fare_amount > 0 and fare_amount < 250
""").df()
print(f"{len(t):,} trips")

fig, axes = plt.subplots(1, 3, figsize=(14, 4), sharex=True, sharey=True)
axes[0].scatter(t.miles, t.fare, s=1, color="#1f4e79")
axes[0].set_title("Plain scatter: one solid blob", loc="left", fontsize=10)
sample = t.sample(20_000, random_state=1)
axes[1].scatter(sample.miles, sample.fare, s=3, alpha=0.15, color="#1f4e79")
axes[1].set_title("20,000-trip sample, alpha 0.15", loc="left", fontsize=10)
hb = axes[2].hexbin(t.miles, t.fare, gridsize=60, bins="log", cmap="viridis", mincnt=1)
fig.colorbar(hb, ax=axes[2], label="trips per hexagon (log scale)")
axes[2].set_title("Hexbin of all trips, log colour", loc="left", fontsize=10)
for ax in axes:
    ax.set_xlabel("trip distance, miles")
axes[0].set_ylabel("fare, USD")
fig.tight_layout()
fig.savefig("taxi_overplotting.png", dpi=150)

counts = np.ma.compressed(hb.get_array())          # only hexagons with at least one trip
print(f"hexagons used: {len(counts):,}; busiest holds {int(counts.max()):,} trips; median hexagon {int(np.median(counts))}")

# A trend line, honestly: fit only the metered trips (rate code 1), and say so
metered = t[t.rate == 1]
slope, intercept = np.polyfit(metered.miles, metered.fare, 1)
print(f"metered trips: {len(metered):,}; fare = {intercept:.2f} + {slope:.2f} x miles")
print("share of all trips on rate code 2 (JFK flat fare):", f"{(t.rate == 2).mean():.1%}")

3,253,245 trips
hexagons used: 1,955; busiest holds 679,176 trips; median hexagon 5
metered trips: 2,684,091; fare = 5.78 + 3.79 x miles
share of all trips on rate code 2 (JFK flat fare): 3.1%

Look at the three panels in taxi_overplotting.png:

• Plain scatter. 3.25 million points become a solid wedge. The few hundred trips with fares of $150 or $200 at short distances look as prominent as the core, because a pixel cannot be more than fully inked.

• Sample with transparency. 20,000 random trips with alpha=0.15 show the main relationship clearly and the horizontal line of $70 JFK flat fares from about 15 to 25 miles. But the rare structures (fixed fares at $150 and $200, short trips with huge fares) almost vanish, because a 0.6% sample keeps only a handful of them.

• Hexbin with log colour. Every trip counts. The busiest hexagon holds 679,176 trips and the median used hexagon holds 5; only a log colour scale can show both. The dense diagonal (metered trips), the $70 line and the scattered errors are all visible, each at its real density.

The trend line is fitted on metered trips only (rate code 1, 2,684,091 trips): fare = 5.78 + 3.79 x miles. Including the flat-fare airport trips would bend the line towards $70 and describe neither group. Say what you fitted in the chart's note.

#Common mistakes

• Outcome on the x-axis. Readers read left to right as cause to effect. Put the driver on x.

• A plain scatter of a million points. It is slow to render and visually wrong. Sample, use transparency, or hexbin with a log colour scale.

• Bubble radius proportional to the value. Areas grow with the square; big values look far bigger than they are. Scale area.

• A rainbow of 12 region colours. Use grey plus one highlight, or facets.

• Linear axes for incomes, populations, prices. Everything squashes into a corner. Use a log axis with plain tick labels.

• One trend line through two populations. Fit separately, or fit only the group you describe.

• Labelling every point. Label the ones your text talks about.

• "X drives Y" from a scatter. Write "X and Y move together"; causal claims need an experiment or a careful causal design.

• Forgetting aggregates in country data. World Bank and similar files include regions and income groups as rows. Filter them out with the metadata.

#Interview angle

#Practice task

1. Run wdi_scatter.py with the latest WDI files you download. Your numbers may differ slightly if the World Bank has revised them; note any change.

2. Colour the bubbles by the Region column from the metadata (7 regions), then redo it with grey plus South Asia highlighted. Write one sentence on which is easier to read.

3. Redraw the GDP chart for 2004 and 2024 as two panels with shared axes, and label India in both.

4. In the taxi data, make a hexbin of trip_distance against trip duration in minutes, on log-log axes. What do the trips at the bottom right (long distance, short time) tell you about data quality?

5. Done when: each chart states the scale, labels only the points you discuss, and its title describes the relationship without claiming a cause.

#Further reading

• World Bank World Development Indicators: indicator definitions and the CC BY 4.0 licence.

• matplotlib hexbin reference: gridsize, bins="log" and mincnt.

• matplotlib scatter reference: why s is marker area in points squared.

• Gapminder: the classic interactive version of this chart, useful to compare your result.