Reading distributions with seaborn

Looking at the shape of a column is a different way of seeing it than reading its mean. The mean of life expectancy gives you one number. The shape might be one tight cluster of countries, or it might be two separate humps, a richer world and a poorer one, that happen to average out to the same place. I find I keep wanting to see which it is. So this post is me working through the plots that let you look.

It is a tour of seven seaborn plots. Each one reads the same kind of column from the World Development Indicators dataset: one value per country for some indicator. Each is a different lens on that column. The histogram asks whether the spread is one hump or two. The scatter is for two numbers per country that might move together. The box and violin compare the shape once you split by region. The heatmap collapses a whole matrix into a grid of color. The rug goes all the way back down to the raw ticks. Seven views of the same data, and each one foregrounds a different feature of the shape.

The code below was re-run on a current seaborn (0.12 and later), and the API notes reflect that.

Every figure in this post reads the same kind of column, one value per country, through a different lens. The plot you reach for is a choice about which feature of the shape you want to make visible.

One table, one prep pattern

The whole post runs off a single sqlite file. The data behind it is the World Bank’s World Development Indicators, packaged into one database.sqlite. You connect once, pull three tables, and you are set.

%pylab inline

import pandas as pd
import seaborn as sns
import sqlite3


db_path = './data/world-development-indicators/database.sqlite'

conn = sqlite3.connect(db_path)
db = conn.cursor()
db.execute("SELECT name FROM sqlite_master WHERE type='table';")
print(db.fetchall())

data_countries = pd.read_sql_query('select * from Country',conn)
data_series = pd.read_sql_query('select * from Series',conn)
data_indicators = pd.read_sql_query('select * from Indicators',conn)

data_indicators is a long table. There is one row per CountryName x IndicatorName x Year, with the number itself in Value. That long shape is the substrate for everything below.

The same boilerplate runs before every plot: keep only rows whose country has a real region, take each country's most recent year, then group down to one value per country. Name it once and the per-section data prep is just a change of indicator string.

Most of the data-prep blocks below are the same four moves. You filter data_indicators down to the indicator strings you want. You drop aggregate rows by keeping only countries that carry a real Region. You sort and take each country’s last available year. You group to one row per country. Where a section needs two indicators side by side, it pivots them into columns at the end. I kept the blocks verbatim so each one is runnable on its own. Once you have seen the histogram prep you have seen the pattern, and you can skim the rest for the one line that changes: selected_indicators.

Histogram: one hump or two?

Start with the simplest lens. You have one column of values, life expectancy at birth, and a histogram asking whether the spread is one cluster or several.

selected_indicators = ['Life expectancy at birth, female (years)',
                       'Life expectancy at birth, male (years)',
                       'Life expectancy at birth, total (years)']
countries = data_countries.CountryCode[data_countries.Region!=''].unique()
condition = data_indicators.IndicatorName.isin(selected_indicators)

data_plot = data_indicators.loc[condition,:]
condition = data_plot.CountryCode.isin(countries)
data_plot = data_plot.loc[condition,:]
data_plot.sort_values(['CountryName','IndicatorName','Year'], inplace=True)

data_plot = data_plot.groupby(['CountryName','IndicatorName'], as_index=False).last()
data_plot[['feature','type']] = data_plot['IndicatorName'].str.split(', ',expand=True)
data_plot.reset_index(inplace=True, drop=True)

nbins = 15
sns.set(style="white",
        palette="pastel",
        color_codes=True,
        rc={'figure.figsize':(12,8),
            'figure.dpi':500})

# sns.distplot is removed in seaborn >= 0.12; use sns.histplot with kde=True
sns.histplot(data_plot.Value[data_plot.type=='female (years)'], bins=nbins, kde=True)
sns.histplot(data_plot.Value[data_plot.type=='male (years)'], bins=nbins, kde=True)
sns.histplot(data_plot.Value[data_plot.type=='total (years)'], bins=nbins, kde=True)
plt.legend(['Female', 'Male', 'Total'], bbox_to_anchor=(1.12,1.04))
plt.xlim((25,100))
plt.grid(color='black',linestyle='-.',linewidth=0.25)
plt.title('Life expectancy at birth ( In years )')

This draws three overlaid histograms, female, male and total, each with its KDE curve laid on top via kde=True. The nbins = 15 and plt.xlim((25,100)) are the two knobs that matter. The bin count decides how coarse the shape reads, and the x-limit frames the same 25-to-100 range for all three so they overlay on one axis. One thing to watch with the overlay is occlusion. The code sets no alpha and no element='step', so the filled bars can hide each other where they stack, and the KDE curves do most of the work of telling the three apart. One detail is baked into the code: sns.distplot is gone in seaborn 0.12 and later, so this is histplot(kde=True) rather than the old one-call helper.

Three overlaid histograms with KDE curves showing the spread of female, male, and total life expectancy at birth across countries — Life expectancy at birth, female, male and total, overlaid across countries. A histogram is the first question you ask a column: where does the mass sit, and is it one lump or several? The bin count is the one number you have to defend.

Scatter: when one number per country becomes two

A histogram reads one column. The moment you have two numbers per country, the question changes. It stops being where the mass sits and becomes whether the two move together, and that is a scatter. Here it is male versus female unemployment.

selected_indicators = ['Unemployment, female (% of female labor force)',
                       'Unemployment, male (% of male labor force)']

countries = data_countries.CountryCode[data_countries.Region!=''].unique()
condition = data_indicators.IndicatorName.isin(selected_indicators)

data_plot = data_indicators.loc[condition,:]
condition = data_plot.CountryCode.isin(countries)
data_plot = data_plot.loc[condition,:]
data_plot.sort_values(['CountryName','IndicatorName','Year'], inplace=True)

data_plot = data_plot.groupby(['CountryName','IndicatorName'], as_index=False).last()
data_plot[['feature','type']] = data_plot['IndicatorName'].str.split(', ',expand=True)
data_plot.reset_index(inplace=True, drop=True)
data_plot['type'] = data_plot.type.str.replace(' \(% of male labor force\)','')
data_plot['type'] = data_plot.type.str.replace(' \(% of female labor force\)','')
data_plot = data_plot.pivot_table(values='Value',index='CountryName',columns='type')

sns.set(style="white",
        palette="pastel",
        rc={'figure.figsize':(7,5),
            'figure.dpi':500})

sns.lmplot(x = 'female', y = 'male', data = data_plot, fit_reg=False, x_jitter=1.5, y_jitter=1.5)
plt.xlim((0,40))
plt.ylim((0,40))
plt.grid(color='black', linestyle='-.', linewidth=0.25)
plt.title('Unemployment (% of total)',)
plt.savefig('./plots/02.scatter.png',orientation='landscape',dpi=500);

fit_reg=False keeps lmplot from drawing a regression line, because the point here is the cloud, not a fitted slope. The x_jitter=1.5 and y_jitter=1.5 nudge overlapping points apart, so a stack of countries at the same rate does not hide behind one marker. plt.xlim((0,40)) with plt.ylim((0,40)) puts both axes on the same scale, which is what makes the diagonal readable. Points on the line are countries where male and female unemployment match, and the spread off it is the gap between them.

Scatter plot of male versus female unemployment rate per country on matched axes from zero to forty percent — Male against female unemployment, one point per country. Two numbers per country need two axes, and the diagonal does the work: how far a country sits off the 45-degree line is how far its two rates disagree.

Density: a histogram without the bin argument

Go back to one column, but drop the bins. A kernel density estimate is the histogram’s smooth cousin. You no longer have a bin count to defend, but you pay for it with smoothness the data never had. The bin count does not vanish so much as move. It becomes the bandwidth, hidden inside seaborn’s default Scott rule rather than spelled out, so the lens has a knob too, it is just one you usually never touch. Here it is adult mortality, male and female, side by side.

selected_indicators = ['Mortality rate, adult, female (per 1,000 female adults)',
                       'Mortality rate, adult, male (per 1,000 male adults)']

countries = data_countries.CountryCode[data_countries.Region!=''].unique()
condition = data_indicators.IndicatorName.isin(selected_indicators)

data_plot = data_indicators.loc[condition,:]
condition = data_plot.CountryCode.isin(countries)
data_plot = data_plot.loc[condition,:]
data_plot.sort_values(['CountryName','IndicatorName','Year'], inplace=True)

data_plot = data_plot.groupby(['CountryName','IndicatorName'], as_index=False).last()
data_plot[['feature','type']] = data_plot['IndicatorName'].str.split(', adult, ',expand=True)
data_plot.reset_index(inplace=True, drop=True)
data_plot['type'] = data_plot.type.str.replace(' \(per 1,000 female adults\)','')
data_plot['type'] = data_plot.type.str.replace(' \(per 1,000 male adults\)','')
data_plot = data_plot.pivot_table(values='Value',index='CountryName',columns='type')

sns.set(style="white",
        palette="pastel",
        color_codes=True,
        rc={
            'figure.figsize':(10,6),
            'figure.dpi':200
           })

# seaborn >= 0.12 deprecates positional Series args; use keyword form
sns.kdeplot(data=data_plot, x='male', color='red')
sns.kdeplot(data=data_plot, x='female', color='blue')
plt.grid(color='black',linestyle='-.', linewidth=0.25)
plt.title('Mortality rate')
plt.ylim((0,0.006))
plt.xlim((-100,700))
plt.savefig('./03.density.png');

This draws two kdeplot curves, male in red and female in blue. The plt.ylim((0,0.006)) and plt.xlim((-100,700)) fix the frame so the two densities are comparable. The negative left edge is a tell worth keeping in mind. The kernel happily smears probability mass below zero, into mortality rates that cannot exist. That is the cost of the lens. It buys you a clean curve by inventing a little continuity at the edges. The seaborn 0.12 note here is about the call signature: positional Series arguments are deprecated, so this passes data= and x= by keyword.

Two kernel density curves comparing adult mortality rate distributions for male and female populations across countries — Adult mortality, male against female, as kernel density curves. The KDE frees you from choosing a bin count, at the price of a smooth tail the raw data never had, including mass smeared left of zero. Reach for it when you want the silhouette, not the exact counts.

Box and violin: now group by region

So far every plot has read one ungrouped column. The next question is whether the shape changes when you split countries into groups. Two lenses answer it. A box plot gives you the five-number summary per group, median, quartiles and whiskers, and throws the rest away. A violin gives you the box plus the kernel shape it discarded.

The box and the violin describe the same distribution. The box keeps five numbers and the violin keeps the whole silhouette, with the dashed box drawn inside to show what it still carries. If a group hides two clusters, the box will look identical to a single-hump group and the violin will not.

First the box. This one is merchandise trade as a share of GDP, split by world region.

selected_indicators = ['Merchandise trade (% of GDP)']

countries = data_countries.CountryCode[data_countries.Region!=''].unique()
condition = data_indicators.IndicatorName.isin(selected_indicators)

data_plot = data_indicators.loc[condition,:]
condition = data_plot.CountryCode.isin(countries)
data_plot = data_plot.loc[condition,:]
data_plot.sort_values(['CountryName','IndicatorName','Year'], inplace=True)

data_plot = data_plot.groupby(['CountryName','IndicatorName'], as_index=False).last()
data_plot.reset_index(inplace=True, drop=True)
data_plot['Region'] = data_plot.merge(right=data_countries,on='CountryCode',how='left')['Region']

columns_order = sorted(data_plot.Region.unique())

sns.set(style="white",
        palette="pastel",
        color_codes=True,
        rc={
            'figure.figsize':(10,6),'figure.dpi':200
           })

sns.boxplot(x='Region',
            y='Value',
            palette='autumn',
            order=columns_order,
            width=0.4,
            fliersize=3,
            data=data_plot);
plt.grid(color='black',linestyle='-.', linewidth=0.25)
plt.xticks(rotation=30)
plt.title('Merchandise trade')
plt.ylabel('% of GDP');
plt.savefig('./04.boxplot.png');

width=0.4 keeps the boxes slim so the region labels breathe, and fliersize=3 shrinks the outlier dots that would otherwise crowd the high-trade regions. You get one box per region, each a five-number summary. The comparison you get for free is which regions trade more relative to their size and which spread wider. The sorted(data_plot.Region.unique()) is deliberate. %pylab inline pulls numpy.sort into the namespace as a bare sort, so I reach for the builtin sorted to be explicit about fixing the column order across every grouped plot below.

Box plot of merchandise trade as a percent of GDP for each world region with outliers shown as small dots — Merchandise trade as a share of GDP, one box per region. The box is the fast way to compare medians and spreads across groups, but it shows you a single hump even when there are two.

Now the violin, which keeps the shape the box drops. This one also splits each region in half, gaseous against liquid CO2 emissions.

selected_indicators = [ 'CO2 emissions from gaseous fuel consumption (% of total)',
                        'CO2 emissions from liquid fuel consumption (% of total)']

countries = data_countries.CountryCode[data_countries.Region!=''].unique()
condition = data_indicators.IndicatorName.isin(selected_indicators)

data_plot = data_indicators.loc[condition,:]
condition = data_plot.CountryCode.isin(countries)
data_plot = data_plot.loc[condition,:]
data_plot.sort_values(['CountryName','IndicatorName','Year'], inplace=True)

data_plot = data_plot.groupby(['CountryName','IndicatorName'], as_index=False).last()
data_plot.reset_index(inplace=True, drop=True)
data_plot['Region'] = data_plot.merge(right=data_countries,on='CountryCode',how='left')['Region']

import matplotlib.patches as mpatches

columns_order = sorted(data_plot.Region.unique())

sns.set(style="white",
        palette="pastel",
        color_codes=True,
        rc={
            'figure.figsize':(16,10),'figure.dpi':250
           })
sns.violinplot(x ='Region',
               y='Value',
               hue='IndicatorName',
               linewidth=0.25,
               inner="quart",
               palette={"CO2 emissions from gaseous fuel consumption (% of total)": "y",
                        "CO2 emissions from liquid fuel consumption (% of total)": "b"},
               data=data_plot,
               split=True)
plt.grid(color='black',linestyle='-.', linewidth=0.25)
plt.title('CO$_2$ emission')
plt.ylabel('% of total')

gas_patch = mpatches.Patch(color='yellow', label='Gaseous',alpha=0.5)
liquid_patch = mpatches.Patch(color='skyblue', label='Liquid')
plt.legend(handles=[gas_patch, liquid_patch], bbox_to_anchor=(0.2, 0.99), fontsize='x-large')
plt.savefig('./plots/05.violinplot.png', dpi=250, bbox_inches='tight');

split=True with inner="quart" is the trick that earns the violin its keep. Gaseous goes on one side, liquid on the other, with quartile lines drawn inside each half so you still get the box-style summary without a second plot. Where a region’s silhouette bulges in two places, you are looking at a cluster a box would have flattened into one.

Split violin plot per region comparing the distribution of gaseous and liquid CO2 emissions as a percent of total — Gaseous against liquid CO2 emissions, one split violin per region. Same five-number summary as a box, drawn as the quartile lines inside, plus the full silhouette around it. The shape is exactly what the box throws away.

Heatmap: a whole matrix at once

The plots so far read a column, or a column split into groups. A heatmap is for when the data is already a matrix, region by trade destination, and you have run out of position axes. Color becomes the only channel left.

selected_indicators_export = [
    'Merchandise exports to developing economies in East Asia & Pacific (% of total merchandise exports)',
    'Merchandise exports to developing economies in Latin America & the Caribbean (% of total merchandise exports)',
    'Merchandise exports to developing economies in Middle East & North Africa (% of total merchandise exports)',
    'Merchandise exports to developing economies in South Asia (% of total merchandise exports)',
    'Merchandise exports to developing economies in Sub-Saharan Africa (% of total merchandise exports)',
    'Merchandise exports to developing economies outside region (% of total merchandise exports)',
    'Merchandise exports to developing economies within region (% of total merchandise exports)',
    'Merchandise exports to economies in the Arab World (% of total merchandise exports)',
    'Merchandise exports to high-income economies (% of total merchandise exports)'
]

selected_indicators_imports = [
    'Merchandise imports from developing economies in East Asia & Pacific (% of total merchandise imports)',
    'Merchandise imports from developing economies in Latin America & the Caribbean (% of total merchandise imports)',
    'Merchandise imports from developing economies in Middle East & North Africa (% of total merchandise imports)',
    'Merchandise imports from developing economies in South Asia (% of total merchandise imports)',
    'Merchandise imports from developing economies in Sub-Saharan Africa (% of total merchandise imports)',
    'Merchandise imports from developing economies outside region (% of total merchandise imports)',
    'Merchandise imports from developing economies within region (% of total merchandise imports)',
    'Merchandise imports from economies in the Arab World (% of total merchandise imports)',
    'Merchandise imports from high-income economies (% of total merchandise imports)'
]

countries = data_countries.CountryCode[data_countries.Region!=''].unique()
condition = data_indicators.IndicatorName.isin(selected_indicators_export)
data_plot = data_indicators.loc[condition,:]
condition = data_plot.CountryCode.isin(countries)
data_plot = data_plot.loc[condition,:]
data_plot.sort_values(['CountryName','IndicatorName','Year'], inplace=True)
data_plot = data_plot.groupby(['CountryName','IndicatorName'], as_index=False).last()
data_plot.reset_index(inplace=True, drop=True)
data_plot['Region'] = data_plot.merge(right=data_countries,on='CountryCode',how='left')['Region']
data_export = data_plot.pivot_table(values='Value',columns='Region',index='IndicatorName')


countries = data_countries.CountryCode[data_countries.Region!=''].unique()
condition = data_indicators.IndicatorName.isin(selected_indicators_imports)
data_plot = data_indicators.loc[condition,:]
condition = data_plot.CountryCode.isin(countries)
data_plot = data_plot.loc[condition,:]
data_plot.sort_values(['CountryName','IndicatorName','Year'], inplace=True)
data_plot = data_plot.groupby(['CountryName','IndicatorName'], as_index=False).last()
data_plot.reset_index(inplace=True, drop=True)
data_plot['Region'] = data_plot.merge(right=data_countries,on='CountryCode',how='left')['Region']
data_import = data_plot.pivot_table(values='Value',columns='Region',index='IndicatorName')

sns.set(style="white",
        color_codes=True,
        rc={
            'figure.figsize':(20,8),
            'figure.dpi':250
           })
fig, (imports, exports) = plt.subplots(1, 2, sharex=True)

im1 = sns.heatmap(data_import.loc[:,xlabels],
                  ax=imports,
                  center=50,
                  cbar=False,
                  cmap="YlGnBu")
imports.set_yticklabels(ylabels)
imports.set_ylabel('')
imports.set_xlabel('')
imports.set_title('Imports :');

im2 = sns.heatmap(data_export.loc[:,xlabels],
                  ax=exports,
                  center=50,
                  yticklabels=False,
                  cbar=False,
                  cmap="YlGnBu")
exports.set_ylabel('')
exports.set_xlabel('')
exports.set_title('Exports :');
fig.subplots_adjust(wspace=0.05, hspace=0)

mappable = im1.get_children()[0]
fig.colorbar(mappable, ax = [imports,exports],orientation = 'vertical')
plt.savefig('./plots/06.heatmap.png', dpi=250, bbox_inches='tight');

Two heatmaps share a colorbar, imports on the left and exports on the right. center=50 anchors the colormap so the eye reads above and below the halfway mark, and wspace=0.05, hspace=0 pulls the panels tight so they read as one figure. One gap is worth flagging. This block references xlabels and ylabels that I defined further up in the original notebook and never carried into the post. They are just the column order and the row labels for the two pivot tables, so the quickest stub is to derive them from what is already in scope:

xlabels = sorted(data_import.columns)
ylabels = [name.split(' (% ')[0] for name in data_import.index]

Side by side heatmaps of merchandise imports and exports by region sharing one color bar centered at fifty — Merchandise imports and exports by region, two matrices sharing a colorbar. When the data is region by destination, position is spent on both axes and color carries the value. A heatmap trades precise reading for the ability to scan a whole grid at a glance.

Rug: the raw data under everything

Close the tour by going all the way back down. Every plot above summarized, binned or smoothed the column. The rug does none of that. It draws every country as a single tick, the raw marks the other six plots were built on top of.

selected_indicators = ['Merchandise trade (% of GDP)']

countries = data_countries.CountryCode[data_countries.Region!=''].unique()
condition = data_indicators.IndicatorName.isin(selected_indicators)

data_plot = data_indicators.loc[condition,:]
condition = data_plot.CountryCode.isin(countries)
data_plot = data_plot.loc[condition,:]
data_plot.sort_values(['CountryName','IndicatorName','Year'], inplace=True)

data_plot = data_plot.groupby(['CountryName','IndicatorName'], as_index=False).last()
data_plot.reset_index(inplace=True, drop=True)
data_plot['Region'] = data_plot.merge(right=data_countries,on='CountryCode',how='left')['Region']

columns_order = sorted(data_plot.Region.unique())

sns.set(style="white",
        palette="pastel",
        color_codes=True,
        rc={
            'figure.figsize':(12,8),'figure.dpi':500
           })

g = sns.FacetGrid(data_plot,
                  col="Region",
                  col_wrap=4,
                  col_order=columns_order,subplot_kws={'ylim':(0,0.02)})
# sns.distplot removed in seaborn >= 0.12; use sns.kdeplot with rug=True
g.map(sns.kdeplot, "Value")
g.map(sns.rugplot, "Value");
plt.savefig('./plots/07.rugplot.png', dpi=500, bbox_inches='tight');

col_wrap=4 lays the regions out four to a row as a small-multiples grid. subplot_kws={'ylim':(0,0.02)} pins every panel to the same y-range so the densities are comparable across regions, and each panel layers a kdeplot curve over a rugplot of the actual ticks. The rug is the reality check on the curve above it. Where the density says smooth bump, the rug shows you whether that bump rests on twenty countries or three. And that reality check is also the limit of the whole tour. None of these seven plots test anything. They form hypotheses, they do not confirm them. A violin that bulges twice is a reason to go check, not a finding. What the rug keeps honest is the eye that the other six lenses train: glance at the raw ticks whenever a smooth curve looks too clean to trust.

Small multiples of merchandise trade per region, each panel a density curve over a rug of individual country ticks — Merchandise trade by region, a density curve over a rug of individual country ticks, faceted four to a row. The rug is every country as one mark, the raw data the previous six plots were summarizing. Glance at it whenever a smooth curve looks too good.

Which one to reach for

The heuristic is short. One column, reach for a histogram or a density. Two numbers per country, a scatter. A column that splits into groups, a box for the fast comparison or a violin when you suspect the box is hiding a second hump. A full matrix, a heatmap.

The deeper point is that the plot you pick decides what you are allowed to notice. A box can never show you bimodality. A density can never tell you it rests on three points. The mean cannot tell you anything at all. So the move is to look through more than one lens before you believe what you see.

Share on

Bluesky Facebook LinkedIn Mastodon X (formerly Twitter)