25.1 Exploratory Data Analysis as Inquiry

"I never start with a hypothesis. I start with a histogram. The hypothesis comes later, usually uninvited, around the third violin plot."

A Quartile With Philosophical Leanings

1. Tukey's Revolution: Analysis Before Models

Consider a researcher who fits a linear model to 10,000 patient records and publishes a significant result. Only later does she discover that the outcome variable is bimodal (having two distinct peaks rather than one): two distinct populations hid in a single column, and a five-second histogram would have caught it. In 1977, John Tukey published Exploratory Data Analysis precisely to prevent this kind of blindness. Before Tukey, the workflow was rigid: specify a model, estimate its parameters, test a hypothesis. Tukey argued that this approach was backwards. If you do not know the shape of your data, any model you specify is a guess, and any hypothesis you test is a question you may not have needed to ask.

Tukey's core insight was that data analysis should be detective work, not judicial proceedings. The detective looks for clues without knowing what crime was committed. The judge evaluates evidence for a specific charge. Both are necessary, but the detective must go first. As discussed in Chapter 1, discovery is search through a space of possibilities; EDA is the reconnaissance phase that maps the terrain before you commit to a search direction.

What EDA Is, Why It Matters, and When to Use It

What. Exploratory data analysis is the systematic examination of a dataset's structure, distributions, relationships, and anomalies, conducted before formal modeling or hypothesis testing.

Why. Models encode assumptions (linearity, normality, independence). EDA reveals whether those assumptions hold. More importantly, EDA surfaces unexpected patterns that suggest new questions: the outlier cluster that turns out to be a distinct cell type, the bimodal distribution that reveals two populations mixed together, the correlation that vanishes once you condition on a third variable.

How. Through a cycle of question, visualization, computation, and refinement. Each pass sharpens the questions. The first pass asks broad questions ("What does the distribution look like?"). Later passes ask targeted ones ("Why do samples from batch 3 cluster separately?").

When. Always. Every dataset deserves an exploratory pass, whether it contains 50 rows or 50 million. The tools change with scale; the intellectual discipline does not. In short: the five-second histogram you skip today is the retracted paper you publish tomorrow.

2. The EDA Inquiry Cycle

Effective EDA follows a structured cycle that mirrors the scientific method introduced in Chapter 2. The cycle has four phases, and you repeat it until you run out of questions or time. Figure 25.1 illustrates this iterative loop.

Phase 1: Orient. What are the variables? What are their types (continuous, categorical, ordinal)? How many observations? How much missingness? This is the "look before you leap" phase.

Phase 2: Describe. Compute summary statistics and visualize univariate distributions. Histograms, kernel density estimates (smooth curves that approximate a variable's probability distribution), box plots. The goal is to understand each variable in isolation.

Phase 3: Relate. Explore bivariate and multivariate relationships. Scatter plots, correlation matrices (tables showing the Pearson correlation coefficient, where r = 1 means perfect positive linear association, r = -1 means perfect negative, and r = 0 means no linear relationship), pair plots. The goal is to find associations, clusters, and outliers.

Phase 4: Refine. Ask "why?" about every surprise. Subset the data, condition on confounders (variables that influence both the predictor and the outcome, potentially creating spurious associations), compute group-wise statistics. Each answer generates new questions for the next cycle. Figure 25.1.1 illustrates EDA Inquiry Cycle.

Skipping this cycle has real consequences. In one reported case, a widely cited genomics study was retracted after reviewers discovered that batch effects in the sequencing pipeline had created artificial clusters, a problem that a single stratified box plot during Phase 2 (Describe) would have revealed before any model was trained.

The following example applies this cycle to the Wisconsin Breast Cancer dataset, a classic benchmark with 30 numeric features computed from digitized images of cell nuclei.

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

# Phase 1: Orient
data = load_breast_cancer()
df = pd.DataFrame(data.data, columns=data.feature_names)
df["diagnosis"] = pd.Categorical(
    data.target_names[data.target]  # 'malignant' or 'benign'
)

print(f"Shape: {df.shape}")
print(f"Features: {len(data.feature_names)} numeric")
print(f"Classes: {df['diagnosis'].value_counts().to_dict()}")
print(f"\nMissing values: {df.isnull().sum().sum()}")

# Phase 2: Describe (summary of key features)
summary_cols = ["mean radius", "mean texture", "mean smoothness",
                "mean compactness", "mean concavity"]
print(f"\n{df[summary_cols].describe().round(3)}")

# Phase 3: Relate (correlation among mean features)
mean_features = [c for c in df.columns if c.startswith("mean")]
corr_matrix = df[mean_features].corr()
high_corr = (corr_matrix.abs() > 0.8) & (corr_matrix != 1.0)
pairs = [(corr_matrix.index[i], corr_matrix.columns[j],
          corr_matrix.iloc[i, j])
         for i in range(len(corr_matrix))
         for j in range(i+1, len(corr_matrix))
         if high_corr.iloc[i, j]]

print(f"\nHighly correlated feature pairs (|r| > 0.8):")
for f1, f2, r in sorted(pairs, key=lambda x: -abs(x[2])):
    print(f"  {f1} <-> {f2}: r = {r:.3f}")

Shape: (569, 31)
Features: 30 numeric
Classes: {'benign': 357, 'malignant': 212}

Missing values: 0

Highly correlated feature pairs (|r| > 0.8):
  mean radius <-> mean perimeter: r = 0.998
  mean radius <-> mean area: r = 0.987
  mean perimeter <-> mean area: r = 0.987
  mean compactness <-> mean concavity: r = 0.883
  mean concavity <-> mean concave points: r = 0.921

Phase 4 (Refine) asks: do these correlations hold equally within the malignant and benign subgroups? A correlation that exists only in one class is more scientifically interesting than one that is universal.

# Phase 4: Refine - check if correlations differ by diagnosis
for label in ["malignant", "benign"]:
    subset = df[df["diagnosis"] == label]
    r = subset["mean compactness"].corr(subset["mean concavity"])
    print(f"  Compactness-concavity correlation ({label}): r = {r:.3f}")

  Compactness-concavity correlation (malignant): r = 0.737
  Compactness-concavity correlation (benign): r = 0.675

3. Modern Tools for Exploratory Wrangling

Tukey worked with pencil and paper. We have better options. The Python ecosystem offers three distinct tools for data wrangling during EDA, each optimized for a different scale and workflow.

3.1 pandas: The Default Choice

pandas remains the most widely used DataFrame library in Python. Its API is mature, its ecosystem is vast (every plotting library speaks pandas), and for datasets under a few million rows, it is fast enough. The code in Listings 25.1 and 25.2 used pandas, and for most exploratory work, it is the right tool.

The main limitation is memory. pandas loads the entire dataset into RAM as a contiguous array, so a 10 GB CSV requires roughly 10 GB of free memory (often more, due to intermediate copies). For datasets that approach or exceed available RAM, you need a different approach.

3.2 Polars: Speed Without Complexity

Polars is a DataFrame library written in Rust that provides a pandas-like API with dramatically better performance. Three design decisions drive the speedup: lazy evaluation (Polars builds a query plan describing the operations before executing any of them), columnar memory layout (cache-friendly for analytical queries), and automatic parallelism (all CPU cores work without explicit threading).

import polars as pl

# Read the same data into Polars
df_pl = pl.DataFrame({
    name: data.data[:, i] for i, name in enumerate(data.feature_names)
})
df_pl = df_pl.with_columns(
    pl.Series("diagnosis", data.target_names[data.target])
)

# Lazy evaluation: build the query plan, then execute
result = (
    df_pl.lazy()
    .group_by("diagnosis")
    .agg([
        pl.col("mean radius").mean().alias("avg_radius"),
        pl.col("mean texture").mean().alias("avg_texture"),
        pl.col("mean concavity").mean().alias("avg_concavity"),
        pl.col("mean radius").std().alias("std_radius"),
        pl.len().alias("count"),
    ])
    .sort("diagnosis")
    .collect()  # execute the plan
)

print(result)

shape: (2, 6)
┌────────────┬────────────┬─────────────┬───────────────┬────────────┬───────┐
│ diagnosis ┆ avg_radius ┆ avg_texture ┆ avg_concavity ┆ std_radius ┆ count │
│ ---       ┆ ---        ┆ ---         ┆ ---           ┆ ---        ┆ ---   │
│ str       ┆ f64        ┆ f64         ┆ f64           ┆ f64        ┆ u32   │
╞═══════════╪════════════╪═════════════╪═══════════════╪════════════╪═══════╡
│ benign    ┆ 12.147     ┆ 17.915      ┆ 0.046         ┆ 1.781      │ 357   │
│ malignant ┆ 17.463     ┆ 21.605      ┆ 0.161         ┆ 3.204      │ 212   │
└───────────┴────────────┴─────────────┴───────────────┴────────────┴───────┘

3.3 DuckDB: SQL for Exploratory Analysis

DuckDB is an embedded analytical database engine. Unlike PostgreSQL or MySQL, it runs inside your Python process with no server, no configuration, and no network overhead. You query CSV files, Parquet files, and pandas DataFrames directly with SQL. For EDA, this is powerful because many exploratory questions are naturally expressed as SQL queries, especially those involving grouping, filtering, windowing, and joining.

import duckdb

# DuckDB can query pandas DataFrames directly - no import needed
result = duckdb.sql("""
    SELECT
        diagnosis,
        COUNT(*) AS n,
        AVG("mean radius") AS avg_radius,
        PERCENTILE_CONT(0.5) WITHIN GROUP
            (ORDER BY "mean concavity") AS median_concavity,
        AVG("mean radius") / AVG("mean smoothness") AS radius_smoothness_ratio
    FROM df
    GROUP BY diagnosis
    ORDER BY diagnosis
""")
print(result.df())  # convert result to pandas for display

    diagnosis    n  avg_radius  median_concavity  radius_smoothness_ratio
0     benign  357   12.146524          0.036360               126.069
1  malignant  212   17.462830          0.151700               176.988

Library Shortcut: DuckDB for Large-File EDA

For Parquet files too large to fit in memory, DuckDB can query them directly without loading: duckdb.sql("SELECT * FROM 'measurements.parquet' WHERE temperature > 300"). This single line replaces a typical multi-step pipeline of chunked pandas reads, filters, and concatenations. DuckDB handles predicate pushdown, column pruning, and parallel scanning internally.

4. The Grammar of Graphics for Discovery

Wrangling tools prepare data for inspection, but the patterns that matter most often reveal themselves visually.

Visualization is the language of EDA. But not all visualizations are equal. A well-chosen plot answers a specific question; a poorly chosen plot hides the answer behind irrelevant decoration. The grammar of graphics (a formal system that decomposes any chart into independent, combinable layers such as data, aesthetics, and geometry), formalized by Leland Wilkinson and implemented in libraries like ggplot2 and Plotly, provides a systematic framework for mapping data to visual properties.

The grammar decomposes every chart into layers: data (the DataFrame), aesthetics (which columns map to x, y, color, size), geometry (points, lines, bars, densities), facets (subplots by category), and statistics (transformations like binning or smoothing). By thinking in these terms, you can systematically explore relationships without falling into the trap of "trying random charts until one looks interesting."

Formally, the grammar treats each visualization as a pipeline. Raw data flows through an aesthetic mapping that assigns columns to visual channels (position, color, size), then through a geometric renderer that chooses marks (points, bars, lines), and optionally through statistical transforms (binning, smoothing) and faceting (subplots by category). This composability replaces trial-and-error chart selection with principled design: each visual decision maps to a specific data property. Reach for a grammar-based library when exploring multivariate relationships; fall back to plain matplotlib for a quick single-variable histogram.

import plotly.express as px

# Scatter matrix of key features, colored by diagnosis
fig = px.scatter_matrix(
    df,
    dimensions=["mean radius", "mean texture",
                 "mean concavity", "mean symmetry"],
    color="diagnosis",
    opacity=0.5,
    title="Pairwise Feature Relationships by Diagnosis",
    labels={col: col.replace("mean ", "") for col in df.columns},
)
fig.update_traces(diagonal_visible=False, marker=dict(size=3))
fig.update_layout(width=800, height=800)
fig.show()

The scatter matrix in Listing 25.5 reveals at a glance which feature pairs separate the two classes. Mean radius and mean concavity show clear separation; mean texture and mean symmetry are nearly useless alone. Notice that malignant tumors average 0.161 concavity versus 0.046 for benign, a 3.5x difference, yet this dramatic gap is invisible in any single summary statistic table and only leaps out when color encodes diagnosis in a scatter plot. This observation, which emerges in seconds from a four-line plot, would take pages of numerical summaries to convey.

5. Systematic Anomaly Screening

While the grammar of graphics equips you to see patterns, it is equally important to spot the observations that do not fit any pattern at all.

Every dataset contains surprises. Some are errors (a temperature reading of negative 1000 degrees). Some are genuine anomalies (a cell with gene expression levels ten standard deviations above the mean). EDA must distinguish between the two, because errors corrupt your models while genuine anomalies may be the most interesting observations in the dataset.

A systematic screening protocol begins with univariate outlier detection, progresses to multivariate methods, and finishes with domain-specific sanity checks.

from scipy import stats

# Univariate screening: flag features with extreme values
def screen_outliers(df, features, z_threshold=4.0):
    """Flag observations with |z-score| > threshold in any feature."""
    results = {}
    for feat in features:
        z_scores = np.abs(stats.zscore(df[feat].values))
        outlier_idx = np.where(z_scores > z_threshold)[0]
        if len(outlier_idx) > 0:
            results[feat] = {
                "count": len(outlier_idx),
                "max_z": z_scores[outlier_idx].max(),
                "indices": outlier_idx.tolist()
            }
    return results

mean_features = [c for c in df.columns if c.startswith("mean")]
outliers = screen_outliers(df, mean_features)

print("Features with extreme observations (|z| > 4.0):")
for feat, info in sorted(outliers.items(), key=lambda x: -x[1]["max_z"]):
    print(f"  {feat}: {info['count']} outliers, max |z| = {info['max_z']:.1f}")

Features with extreme observations (|z| > 4.0):
  mean area: 5 outliers, max |z| = 5.2
  mean perimeter: 3 outliers, max |z| = 4.8
  mean radius: 2 outliers, max |z| = 4.6
  mean concave points: 1 outliers, max |z| = 4.1

Univariate screening catches obvious outliers, but some anomalies are multivariate: a data point that is unremarkable on any single feature but bizarre when you consider the combination. We will address multivariate anomaly detection systematically in Chapter 30. For now, the scatter matrix from Listing 25.5 provides a visual check: points that fall outside the main clusters in multiple pairwise views deserve investigation.

6. From Visual to Computational Exploration

Visual exploration has limits. You can visualize two dimensions at a time (three with color), but real datasets have dozens or hundreds of features. Beyond five or six dimensions, scatter matrices become unwieldy, and the number of pairwise relationships grows quadratically. This is where computational methods, specifically clustering and dimensionality reduction, take over from visual inspection.

The transition from visual to computational EDA marks a shift in how we think about patterns. A histogram shows you a distribution; a clustering algorithm finds the distributions automatically. A scatter plot shows you a projection; a dimensionality reduction algorithm finds the best projection. These methods do not replace visual EDA; they extend it into spaces your eyes cannot reach.

Section 25.2 develops the mathematical foundations of clustering and dimensionality reduction, building the tools that let us explore datasets with hundreds or thousands of features as naturally as Tukey explored datasets with five.