15 - Exploratory Data Analysis (EDA) with Python

They hand you a data file and ask for a model. Before choosing the algorithm, you need to know what you are holding: how many rows, which columns are numbers and which are labels, where values are missing, what the distributions look like, which variables move together and whether some value is so far from the others that it drives every calculation. Exploratory data analysis (EDA) is this work. The NIST statistics handbook describes it as an approach, not a list of techniques, that relies mostly on graphics to bring out the structure of the data, spot outliers and test underlying assumptions, postponing the choice of model; its founding work is Tukey's 1977 book. In this article, the last in the Probability and Statistics series, you learn to do it in Python with pandas and matplotlib on a real dataset, and you see with numbers why summary statistics alone are never enough. It is the step that comes before the tests, intervals and models of the previous articles.

Why summarising is not enough: Anscombe's quartet

In 1973 Francis Anscombe published four small datasets, of 11 (x, y) pairs each, built to have almost the same summary statistics and completely different shapes. The listing rebuilds them from the published values and computes for each the mean and variance of x and y, Pearson's correlation and the least-squares line. Then it looks at what those numbers hide: in the second set it measures how much better a parabola describes the data than the line; in the third it compares the least-squares slope with the robust Theil-Sen one (the median of the slopes of all pairs of points); in the fourth it counts the distinct values of x.

import numpy as np
from scipy import stats

# Anscombe's quartet (1973): four datasets of 11 points
x123 = np.array([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5], float)
x4 = np.array([8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8], float)
data = {
    "I":   (x123, np.array([8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68])),
    "II":  (x123, np.array([9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74])),
    "III": (x123, np.array([7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73])),
    "IV":  (x4,   np.array([6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89])),
}
print("set   mean x  mean y  var x  var y   r     line y = a + b x")
for name, (x, y) in data.items():
    b, a = np.polyfit(x, y, 1)
    print(f"{name:5s} {x.mean():6.2f} {y.mean():8.2f} {x.var(ddof=1):6.2f} {y.var(ddof=1):6.2f} "
          f"{np.corrcoef(x, y)[0, 1]:6.3f}   y = {a:.2f} + {b:.3f} x")

# What the eye sees and the statistics do not
x, y = data["II"]
quad = np.polyfit(x, y, 2)
r2_lin = np.corrcoef(x, y)[0, 1] ** 2
r2_quad = 1 - np.sum((y - np.polyval(quad, x)) ** 2) / np.sum((y - y.mean()) ** 2)
print(f"II : R2 line {r2_lin:.3f}  R2 parabola {r2_quad:.3f}  (the relationship is curved)")
x, y = data["III"]
ts = stats.theilslopes(y, x)
b_ols, a_ols = np.polyfit(x, y, 1)
print(f"III: least-squares slope {b_ols:.3f}  Theil-Sen slope {ts.slope:.3f}  (one point drags the line)")
x, y = data["IV"]
print(f"IV : distinct values of x = {sorted(set(x.tolist()))}  without the point x=19 all the x are equal: slope undefined")

# Output:
# set   mean x  mean y  var x  var y   r     line y = a + b x
# I       9.00     7.50  11.00   4.13  0.816   y = 3.00 + 0.500 x
# II      9.00     7.50  11.00   4.13  0.816   y = 3.00 + 0.500 x
# III     9.00     7.50  11.00   4.12  0.816   y = 3.00 + 0.500 x
# IV      9.00     7.50  11.00   4.12  0.817   y = 3.00 + 0.500 x
# II : R2 line 0.666  R2 parabola 1.000  (the relationship is curved)
# III: least-squares slope 0.500  Theil-Sen slope 0.346  (one point drags the line)
# IV : distinct values of x = [8.0, 19.0]  without the point x=19 all the x are equal: slope undefined

The four sets have means 9.00 and 7.50, variances of x equal to 11.00 and of y between 4.12 and 4.13, correlation 0.816-0.817 and the same line, y = 3.00 + 0.500 x. Yet the first is a linear cloud with noise, in the second the relationship is a perfect curve (R² of the line 0.666, of the parabola 1.000), in the third a single outlier drags the slope to 0.500 while the robust estimate gives 0.346, in the fourth all the x equal 8 except one at 19, and without that point the slope would not exist. No number in that table would have made you suspect any of this; a scatter plot would. This is why EDA begins by looking at the data before summarising them.

A first look at a real dataset

As an example I use the wine dataset, included in scikit-learn: 178 wines from one Italian region, 13 chemical measurements for each and a class among three producers (the description is in the documentation of load_wine). The first move is always the same: the shape of the table, the column types and the missing values. Then a numerical summary of the variables, with mean and median side by side, because a gap between the two is already a hint of skewness. Finally Tukey's fence: any value that falls more than 1.5 times the interquartile range (IQR = Q3 − Q1) below the first quartile or above the third is considered suspect, the same rule by which matplotlib draws the whiskers of boxplots. Skewness is positive if the right tail is longer.

import numpy as np
import pandas as pd
from sklearn.datasets import load_wine

wine = load_wine(as_frame=True)
df = wine.frame.rename(columns={"target": "class"})
print("shape:", df.shape, "  missing values:", int(df.isna().sum().sum()))
print("types:", df.dtypes.astype(str).value_counts().to_dict())
print("rows per class:", df["class"].value_counts().sort_index().to_dict())

cols = ["alcohol", "malic_acid", "ash", "magnesium", "proline"]
summary = df[cols].describe().loc[["mean", "std", "min", "50%", "max"]].round(2)
print(summary.to_string())

# Shape of the distributions: skewness and values outside Tukey's "fences" (Q1 - 1.5 IQR, Q3 + 1.5 IQR)
q1, q3 = df[cols].quantile(.25), df[cols].quantile(.75)
iqr = q3 - q1
outside = ((df[cols] < q1 - 1.5 * iqr) | (df[cols] > q3 + 1.5 * iqr)).sum()
print(pd.DataFrame({"skewness": df[cols].skew().round(2), "outside fence": outside}).to_string())

# Output:
# shape: (178, 14)   missing values: 0
# types: {'float64': 13, 'int64': 1}
# rows per class: {0: 59, 1: 71, 2: 48}
#       alcohol  malic_acid   ash  magnesium  proline
# mean    13.00        2.34  2.37      99.74   746.89
# std      0.81        1.12  0.27      14.28   314.91
# min     11.03        0.74  1.36      70.00   278.00
# 50%     13.05        1.87  2.36      98.00   673.50
# max     14.83        5.80  3.23     162.00  1680.00
#             skewness  outside fence
# alcohol        -0.05              0
# malic_acid      1.04              3
# ash            -0.18              3
# magnesium       1.10              4
# proline         0.77              0

The table has 178 rows and 14 columns (13 real numeric ones plus the integer class), no missing values and classes of 59, 71 and 48 wines: the counts match the official description of the dataset. For proline, the mean (746.89) exceeds the median (673.50) by more than 70 units, and the skewness is 0.77, a sign of a right tail; for alcohol, mean and median almost coincide (13.00 and 13.05) and the skewness is −0.05. Tukey's rule flags 3 values outside the fence for malic acid, 3 for ash, 4 for magnesium and none for alcohol and proline. Flagging does not mean removing: a value outside the fence can be a typing error, a measurement in a different unit or a real and rare case, and the three hypotheses call for opposite actions. Before removing it you check where it comes from.

Relationships between variables, and what they do not say

The next step looks at the variables in pairs. Pearson's correlation measures linear tendency and Spearman's (article 13) the monotone one: if they diverge a lot, the relationship is not a straight line or there is some outlier. The listing lists the three most correlated pairs, shows the means of some variables for each class and computes the one-way ANOVA F ratio for proline, that is, the ratio between the variability between classes and that within classes, by hand and with scipy.stats.f_oneway (see ANOVA and multiple comparisons).

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.datasets import load_wine

wine = load_wine(as_frame=True)
df = wine.frame.rename(columns={"target": "class"})

# Relationships: the most correlated pairs of variables (Pearson and Spearman)
feat = df.drop(columns="class")
pearson = feat.corr().where(np.triu(np.ones((13, 13), bool), 1)).stack()
top = pearson.reindex(pearson.abs().sort_values(ascending=False).index).head(3)
for (a, b), r in top.items():
    print(f"{a} ~ {b}: Pearson {r:.3f}  Spearman {feat[a].corr(feat[b], method='spearman'):.3f}")

# Does a variable separate the classes? mean per class and ratio of between-class to within-class variance
means = df.groupby("class")[["alcohol", "proline", "ash"]].mean().round(2)
print(means.to_string())
between = df.groupby("class")["proline"].agg(["mean", "count"])
F = (between["count"] * (between["mean"] - df["proline"].mean()) ** 2).sum() / 2 / (df.groupby("class")["proline"].var() * (between["count"] - 1)).sum() * (len(df) - 3)
groups = [g["proline"].to_numpy() for _, g in df.groupby("class")]
print(f"proline: Fisher F by hand = {F:.1f}   scipy.stats.f_oneway = {stats.f_oneway(*groups).statistic:.1f}")

# Output:
# total_phenols ~ flavanoids: Pearson 0.865  Spearman 0.879
# flavanoids ~ od280/od315_of_diluted_wines: Pearson 0.787  Spearman 0.742
# total_phenols ~ od280/od315_of_diluted_wines: Pearson 0.700  Spearman 0.687
#        alcohol  proline   ash
# class                        
# 0        13.74  1115.71  2.46
# 1        12.28   519.51  2.24
# 2        13.15   629.90  2.44
# proline: Fisher F by hand = 207.9   scipy.stats.f_oneway = 207.9

The three most correlated pairs are all phenols: total_phenols and flavanoids (Pearson 0.865, Spearman 0.879), then flavanoids and od280/od315 (0.787 and 0.742) and total_phenols with od280/od315 (0.700 and 0.687). They are measurements of the same chemical family, hence not independent: a model that uses them all together will have to take this into account. proline varies a lot with the class: mean 1115.71 in class 0 against 519.51 and 629.90 in the others, with F = 207.9 both by hand and with SciPy. This says that the variable separates the classes well; it does not say that it is their cause, and it does not say which class of wine has «more proline» for a reason that the dataset explains.

Simpson's paradox: when the total lies

A trap of EDA is aggregating too early. In Simpson's paradox (described in 1951) a trend that appears in every subgroup reverses, or disappears, in the aggregated data, because the groups do not have the same composition. The following example is constructed, and the numbers are invented: two variants of a page, A and B, and two types of visitor, experts and beginners, the latter much less likely to convert.

import pandas as pd

# INVENTED data for the example: conversions of two variants, split by visitor type
rows = [("A", "expert",   90,  100), ("A", "beginner", 160, 400),
        ("B", "expert",  800, 1000), ("B", "beginner",  10, 100)]
t = pd.DataFrame(rows, columns=["variant", "group", "conversions", "visits"])
t["rate"] = (t["conversions"] / t["visits"]).round(3)
print(t.to_string(index=False))

total = t.groupby("variant")[["conversions", "visits"]].sum()
total["rate"] = (total["conversions"] / total["visits"]).round(3)
print(total.to_string())

# Why: the two variants do not see the same mix of visitors
mix = t.pivot(index="variant", columns="group", values="visits")
print((mix.div(mix.sum(axis=1), axis=0) * 100).round(1).to_string())

# Output:
# variant    group  conversions  visits  rate
#       A   expert           90     100   0.9
#       A beginner          160     400   0.4
#       B   expert          800    1000   0.8
#       B beginner           10     100   0.1
#          conversions  visits   rate
# variant                            
# A                250     500  0.500
# B                810    1100  0.736
# group    beginner  expert
# variant                  
# A            80.0    20.0
# B             9.1    90.9

Variant A converts better than B in both groups (90% against 80% among experts, 40% against 10% among beginners) but worse overall (50.0% against 73.6%). The explanation is in the last table: 80% of A's visitors are beginners, against 9.1% for B. The aggregated comparison measures the variant and the visitor mix together, and does not let you separate them. Which reading is right depends on how the visitors were assigned: if at random, the two mixes should look alike and one this different points to a flaw in the assignment; if not, the comparison by group is more reliable. The way to notice is to explore the data split by the variables that explain its composition. No test after this step fixes a wrong aggregation (see also Test and learn: controlled experiments and A/B tests).

Charts: which one for which question

Numbers summarise, charts show the shape. The choice depends on the question you are asking of the data, and the table follows this criterion:

The listing produces three panels on the wine dataset: the histogram of proline, the boxplots of the same variable by class and the scatter plot of flavanoids against total_phenols coloured by class. It uses matplotlib's Agg backend, which draws to a file without opening windows: this is what you need on a server or in a script. The image is not reproduced in the article; the program's output confirms that the three panels exist and that the last one contains the 178 points.

import matplotlib
matplotlib.use("Agg")                       # no window: save to file
import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import load_wine

df = load_wine(as_frame=True).frame
fig, ax = plt.subplots(1, 3, figsize=(12, 3.6))
ax[0].hist(df["proline"], bins=20, edgecolor="white")          # shape: one hump, tail to the right
ax[0].set(title="proline: distribution", xlabel="proline", ylabel="count")
df.boxplot(column="proline", by="target", ax=ax[1])             # comparison between groups
ax[1].set(title="proline by class", xlabel="class", ylabel="proline")
for c, g in df.groupby("target"):                               # relationship, coloured by group
    ax[2].scatter(g["flavanoids"], g["total_phenols"], s=12, label=f"class {c}")
ax[2].set(title="phenols and flavanoids", xlabel="flavanoids", ylabel="total_phenols"); ax[2].legend()
fig.suptitle("")
fig.tight_layout()
fig.savefig("eda_wine.png", dpi=80)
print("panels:", len(fig.axes), " points in the last chart:", sum(len(c.get_offsets()) for c in ax[2].collections))

# Output:
# panels: 3  points in the last chart: 178

The line df.boxplot(column="proline", by="target", ax=ax[1]) does the comparison between groups without a loop; ax[2].scatter is inside a loop over the classes to assign a colour to each and an entry to the legend. If you remove matplotlib.use("Agg") on a machine without a screen, the program may fail when creating the figure. Opening the image you should find what the numbers have already said: a right tail in proline (skewness 0.77), class 0 with higher values on average (mean 1115.71) and, in the scatter plot, points that roughly follow a line (correlation 0.865). Charts add no facts, they make them easier to grasp and to show to those who do not read tables.

A worked exercise and a test to do on your own

A company measures delivery times in days: 2, 3, 3, 4, 4, 4, 5, 5, 6, 21. The mean is 57/10 = 5.7; the median, the average of the two central values (4 and 4), is 4.0. The first quartile is 3.25 and the third is 5.0 (linear interpolation, as in NumPy), so IQR = 1.75 and the upper fence is 5.0 + 1.5 · 1.75 = 7.625: the 21 is outside, the others are not. Without the 21 the mean falls to 36/9 = 4.0 and the median stays 4.0. A single value moves the mean by 1.7 days and the median by zero: the median is robust, the mean is not. The 21 must be checked (a genuinely very long delivery, or an error?) before choosing between mean and median to describe the typical delay. The numbers were checked with np.percentile.

What you can do now

You can start from shape, types and missing values, describe a variable with mean and median side by side, and flag values outside the fence without removing them by reflex. You can read two different correlations together, remembering that neither says the cause, you can recognise a Simpson's paradox by splitting the data by group, and you can choose the chart based on the question. With the EDA done, the rest of the series becomes the right tool at the right moment: estimators to summarise, hypothesis tests to decide, regression to model, the bootstrap to measure uncertainty. For pandas and matplotlib in detail see Pandas: data tables without surprises and Visualising data with Matplotlib.

The «Probability and Statistics» series is complete. The fifteen articles, in order: probability spaces and Kolmogorov's axioms, discrete variables, continuous variables, multivariate distributions and covariance, the central limit theorem, estimators, hypothesis tests, ANOVA, chi-square, linear regression, Bayesian estimation, Markov chains, non-parametric statistics, bootstrap and this one.